Changeset 772 in Sophya for trunk/SophyaLib/TArray


Ignore:
Timestamp:
Mar 10, 2000, 4:13:22 PM (26 years ago)
Author:
ansari
Message:

Separation MatrixRC et TMatrix, etc ... - Creation de TArray<T> Reza 10/3/2000

Location:
trunk/SophyaLib/TArray
Files:
6 added
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaLib/TArray/tarrinit.cc

    r762 r772  
    55#include "tmatrix.h"
    66#include "tvector.h"
     7#include "fioarr.h"
    78
    89
     
    1718//   Enregistrement des classes PPersist du modules TArray
    1819 
     20
     21  PPRegister(FIO_TArray<uint_1>);
     22  PPRegister(FIO_TArray<uint_2>);
     23  PPRegister(FIO_TArray<int_2>);
     24  PPRegister(FIO_TArray<int_4>);
     25  PPRegister(FIO_TArray<int_8>);
     26  PPRegister(FIO_TArray<uint_8>);
     27  PPRegister(FIO_TArray<r_4>);
     28  PPRegister(FIO_TArray<r_8>);
     29  PPRegister(FIO_TArray< complex<r_4> >);
     30  PPRegister(FIO_TArray< complex<r_8> >);
    1931
    2032  PPRegister(FIO_TMatrix<uint_1>);
  • trunk/SophyaLib/TArray/tmatrix.cc

    r762 r772  
    1 // $Id: tmatrix.cc,v 1.1.1.1 2000-03-02 16:48:24 ansari Exp $
     1// $Id: tmatrix.cc,v 1.2 2000-03-10 15:13:22 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
     
    161161}
    162162
    163 ////////////////////////////////////////////////////////////////
    164 //**** Pour gestion des lignes et des colonnes
    165 
    166 template <class T> TMatrixRC<T> TMatrix<T>::Row(uint_4 r) const
    167 {
    168 TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixRow, r);
    169 return rc;
    170 }
    171 
    172 template <class T> TMatrixRC<T> TMatrix<T>::Col(uint_4 c) const
    173 {
    174 TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixCol, c);
    175 return rc;
    176 }
    177 
    178 template <class T> TMatrixRC<T> TMatrix<T>::Diag() const
    179 {
    180 TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixDiag);
    181 return rc;
    182 }
    183 
    184 ////////////////////////////////////////////////////////////////
    185 //**** Pour inversion
    186 #ifndef M_LN2
    187 #define M_LN2 0.69314718055994530942
    188 #endif
    189 #ifndef LN_MINDOUBLE
    190 #define LN_MINDOUBLE  (M_LN2 * (DMINEXP - 1))
    191 #endif
    192 #ifndef LN_MAXDOUBLE
    193 #define LN_MAXDOUBLE  (M_LN2 * DMAXEXP)
    194 #endif
    195 
    196 r_8 TMatrix<r_8>::GausPiv(TMatrix<r_8>& a, TMatrix<r_8>& b)
    197 // Pivot de Gauss
    198 // * Attention: egcs impose que cette fonction soit mise dans le .cc
    199 //              avant ::Inverse() (car Inverse() l'utilise)
    200 // {TMatrix A(a); TMatrix B(b); return (r_8) TMatrix::GausPiv(A,B);}
    201 {
    202 uint_4 n = a.NRows();
    203 if(n!=b.NRows())
    204   throw(SzMismatchError("TMatrix::GausPiv size mismatch\n"));
    205 // On fait une normalisation un peu brutale...
    206 double vmin=MAXDOUBLE;
    207 double vmax=0;
    208 for(uint_4 iii=0; iii<a.NRows(); iii++)
    209   for(uint_4 jjj=0; jjj<a.NCols(); jjj++) {
    210     double v = TMatrixRC<r_8>::Abs_Value(a(iii,jjj));
    211     if(v>vmax) vmax = v;
    212     if(v<vmin && v>0) vmin = v;
    213 }
    214 double nrm = sqrt(vmin*vmax);
    215 if(nrm > 1.e5 || nrm < 1.e-5) {
    216   a /= nrm;
    217   b /= nrm;
    218   //cout << "normalisation matrice " << nrm << endl;
    219 } else nrm=1;
    220 
    221 double det = 1.0;
    222 if(nrm != 1) {
    223   double ld = a.NRows() * log(nrm);
    224   if (ld <= LN_MINDOUBLE || ld >= LN_MAXDOUBLE) {
    225    // cerr << "TMatrix warning, overflow for det" << endl;
    226   } else {
    227     det = exp(ld);
    228   }
    229 }
    230 
    231 TMatrixRC<r_8> pivRowa(a,TmatrixRow);
    232 TMatrixRC<r_8> pivRowb(b,TmatrixRow);
    233 
    234 for(uint_4 k=0; k<n-1; k++) {
    235   uint_4 iPiv = a.Col(k).IMaxAbs(k);
    236   if(iPiv != k) {
    237     TMatrixRC<r_8> aIPiv(a.Row(iPiv));
    238     TMatrixRC<r_8> aK(a.Row(k));
    239     TMatrixRC<r_8>::Swap(aIPiv,aK);
    240     TMatrixRC<r_8> bIPiv(b.Row(iPiv));
    241     TMatrixRC<r_8> bK(b.Row(k));
    242     TMatrixRC<r_8>::Swap(bIPiv,bK);
    243   }
    244   double pivot = a(k,k);
    245   if (fabs(pivot) < 1.e-50) return 0.0;
    246   //det *= pivot;
    247   pivRowa.SetRow(k); // to avoid constructors
    248   pivRowb.SetRow(k);
    249   for (uint_4 i=k+1; i<n; i++) {
    250     double r = -a(i,k)/pivot;
    251     a.Row(i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
    252     b.Row(i).LinComb(r, pivRowb);
    253   }
    254 }
    255 det *= a(n-1, n-1);
    256 
    257 // on remonte
    258 for(uint_4 kk=n-1; kk>0; kk--) {
    259   double pivot = a(kk,kk);
    260   if (fabs(pivot) <= 1.e-50) return 0.0;
    261   pivRowa.SetRow(kk); // to avoid constructors
    262   pivRowb.SetRow(kk);
    263   for(uint_4 jj=0; jj<kk; jj++) {
    264     double r = -a(jj,kk)/pivot;
    265     a.Row(jj).LinComb(r, pivRowa);
    266     b.Row(jj).LinComb(r, pivRowb);
    267   }
    268 }
    269 
    270 for(uint_4 l=0; l<n; l++) {
    271   if (fabs(a(l,l)) <= 1.e-50) return 0.0;
    272   b.Row(l) /= a(l,l);
    273 }
    274 
    275 return det;
    276 }
    277 
    278 TMatrix<r_8> TMatrix<r_8>::Inverse() const
    279 // Inversion
    280 {
    281 TMatrix<r_8> b(mNc,mNr); TMatrix<r_8> a(*this); b = 1.;
    282 if(fabs(TMatrix<r_8>::GausPiv(a,b)) < 1.e-50)
    283   throw(MathExc("TMatrix Inverse() Singular OMatrix"));
    284 return b;
    285 }
    286 
    287163
    288164///////////////////////////////////////////////////////////
     
    379255#include "tvector.h"
    380256
    381 template <class T> TMatrixRC<T>::TMatrixRC()
    382 : matrix(NULL), data(NULL), index(0), step(0)
    383 {}
    384 
    385 template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
    386 : matrix(&m), data(Org(m,rckind,ind)),
    387   index(ind), step(Step(m,rckind)), kind(rckind)
    388 {
    389 if (kind == TmatrixDiag && m.mNc != m.mNr)
    390   throw(SzMismatchError("TMatrixRC::TMatrixRC(...,TmatrixDiag,...) size mismatch\n"));
    391 }
    392 
    393 template <class T> int_4 TMatrixRC<T>::Next()
    394 {
    395 if (!matrix || kind==TmatrixDiag) return -1;
    396 index++;
    397 if(kind == TmatrixRow) {
    398   if(index > (int_4)matrix->mNr) {index = (int_4)matrix->mNr; return -1;}
    399   data += matrix->mNc;
    400 } else {
    401   if (index > (int_4)matrix->mNc) {index = (int_4)matrix->mNc; return -1;}
    402   data++;
    403 }
    404 return index;
    405 }
    406 
    407 template <class T> int_4 TMatrixRC<T>::Prev()
    408 {
    409 if (!matrix || kind == TmatrixDiag) return -1;
    410 index--;
    411 if(index < 0) {index = 0; return -1;}
    412 if(kind == TmatrixRow) data -= matrix->mNc;
    413 else data--;
    414 return index;
    415 }
    416 
    417 template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
    418 {
    419 if(!matrix) return -1;
    420 if(c<0 || c>(int_4)matrix->mNc) return -1;
    421 kind = TmatrixCol;
    422 index = c;
    423 step = Step(*matrix, TmatrixCol);
    424 data = Org(*matrix, TmatrixCol, c);
    425 return c;
    426 }
    427 
    428 template <class T> int_4 TMatrixRC<T>::SetRow(int_4 r)
    429 {
    430 if(!matrix) return -1;
    431 if(r<0 && r>(int_4)matrix->mNr) return -1;
    432 kind = TmatrixRow;
    433 index = r;
    434 step = Step(*matrix, TmatrixRow);
    435 data = Org(*matrix, TmatrixRow, r);
    436 return r;
    437 }
    438 
    439 template <class T> int_4 TMatrixRC<T>::SetDiag()
    440 {
    441 if (!matrix) return -1;
    442 if (matrix->mNc != matrix->mNr)
    443   throw(SzMismatchError("TMatrixRC::SetDiag size mismatch\n"));
    444 kind = TmatrixDiag;
    445 index = 0;
    446 step = Step(*matrix, TmatrixDiag);
    447 data = Org(*matrix, TmatrixDiag);
    448 return 0;
    449 }
    450 
    451 
    452 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator = (const TMatrixRC<T>& rc)
    453 {
    454 matrix = rc.matrix;
    455 data   = rc.data;
    456 index  = rc.index;
    457 step   = rc.step;
    458 kind   = rc.kind;
    459 return *this;
    460 }
    461 
    462 template <class T> TVector<T> TMatrixRC<T>::GetVect() const
    463 {
    464 TVector<T> v(NElts());
    465 for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
    466 return v;
    467 }
    468 
    469 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
    470 {
    471 if ( NElts() != rc.NElts() )
    472   throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
    473 if ( kind != rc.kind )
    474   throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
    475 for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
    476 return *this;
    477 }
    478 
    479 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
    480 {
    481 if( NElts() != rc.NElts() )
    482   throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
    483 if( kind != rc.kind )
    484   throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
    485 for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
    486 return *this;
    487 }
    488 
    489 
    490 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator *= (T x)
    491 {
    492 for(uint_4 i=0; i<NElts(); i++) (*this)(i) *= x;
    493 return *this;
    494 }
    495 
    496 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator /= (T x)
    497 {
    498 for(uint_4 i=0; i<NElts(); i++) (*this)(i) /= x;
    499 return *this;
    500 }
    501 
    502 
    503 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
    504 {
    505 for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
    506 return *this;
    507 }
    508 
    509 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
    510 {
    511 for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
    512 return *this;
    513 }
    514 
    515 template <class T>
    516 TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
    517 {
    518 if ( NElts() != rc.NElts() )
    519   throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
    520 if ( kind != rc.kind )
    521   throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
    522 for(uint_4 i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;
    523 return *this;
    524 }
    525 
    526 template <class T>
    527 TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
    528 {
    529 if ( NElts() != rc.NElts() )
    530   throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
    531 if ( kind != rc.kind )
    532   throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
    533 for(uint_4 i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;
    534 return *this;
    535 }
    536 
    537 template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
    538 {
    539 if (first>NElts())
    540   throw(SzMismatchError("TMatrixRC::IMaxAbs size mismatch\n"));
    541 uint_4 imax = first;
    542 double vmax = Abs_Value((*this)(first));
    543 for(uint_4 i=first+1; i<NElts(); i++) {
    544   double v = Abs_Value((*this)(i));
    545   if(v > vmax) {vmax = v; imax = i;}
    546 }
    547 return imax;
    548 }
    549 
    550 template <class T>
    551 void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
    552 {
    553 if(rc1.NElts() != rc2.NElts())
    554   throw(SzMismatchError("TMatrixRC::Swap size mismatch\n"));
    555 if(rc1.kind != rc2.kind)
    556   throw(SzMismatchError("TMatrixRC::Swap type mismatch\n"));
    557 if(rc1.data == rc2.data) return;
    558 for(uint_4 i=0; i<rc1.NElts(); i++)
    559   {T tmp = rc1(i); rc1(i) = rc2(i); rc2(i) = tmp;}
    560 }
    561257
    562258///////////////////////////////////////////////////////////////
     
    585281#pragma define_template FIO_TMatrix< complex<r_4> >
    586282#pragma define_template FIO_TMatrix< complex<r_8> >
    587 // Instances gestion lignes/colonnes
    588 #pragma define_template TMatrixRC<uint_1>
    589 #pragma define_template TMatrixRC<uint_2>
    590 #pragma define_template TMatrixRC<int_2>
    591 #pragma define_template TMatrixRC<int_4>
    592 #pragma define_template TMatrixRC<int_8>
    593 #pragma define_template TMatrixRC<uint_4>
    594 #pragma define_template TMatrixRC<uint_8>
    595 #pragma define_template TMatrixRC<r_4>
    596 #pragma define_template TMatrixRC<r_8>
    597 #pragma define_template TMatrixRC< complex<r_4> >
    598 #pragma define_template TMatrixRC< complex<r_8> >
    599283#endif
    600284
     
    623307template class FIO_TMatrix< complex<r_4> >;
    624308template class FIO_TMatrix< complex<r_8> >;
    625 // Instances gestion lignes/colonnes
    626 template class TMatrixRC<uint_1>;
    627 template class TMatrixRC<uint_2>;
    628 template class TMatrixRC<int_2>;
    629 template class TMatrixRC<int_4>;
    630 template class TMatrixRC<int_8>;
    631 template class TMatrixRC<uint_4>;
    632 template class TMatrixRC<uint_8>;
    633 template class TMatrixRC<r_4>;
    634 template class TMatrixRC<r_8>;
    635 template class TMatrixRC< complex<r_4> >;
    636 template class TMatrixRC< complex<r_8> >;
    637309#endif
  • trunk/SophyaLib/TArray/tmatrix.h

    r762 r772  
    1414namespace SOPHYA {
    1515
    16 template <class T> class TVector;
    17 template <class T> class TMatrixRC;
    18 
    1916template <class T>
    2017class TMatrix : public AnyDataObj {
    21   friend class TMatrixRC<T>;
    22   friend class TVector<T>;
    2318public:
    2419
     
    10398  TMatrix<T> Mul(const TMatrix<T>& b) const;
    10499
    105   // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
    106   // operations sur la matrice B
    107   TMatrix<T> Inverse() const;
    108   static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
    109 
    110   // Acces aux rangees et colonnes
    111   TMatrixRC<T> Row(uint_4 r) const;
    112   TMatrixRC<T> Col(uint_4 c) const;
    113   TMatrixRC<T> Diag() const;
    114100
    115101protected:
     
    203189{ FIO_TMatrix<T> fio(&obj);  fio.Read(is);  return(is); }
    204190
    205 /////////////////////////////////////////////////////////////////////////
    206 // Classe de lignes/colonnes de matrices
    207 enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
    208 template <class T>
    209 class TMatrixRC {
    210   friend class TVector<T>;
    211   friend class TMatrix<T>;
    212 public:
    213   TMatrixRC();
    214 
    215   virtual ~TMatrixRC() {}
    216 
    217   int_4 Next();
    218   int_4 Prev();
    219   int_4 SetCol(int_4 c);
    220   int_4 SetRow(int_4 r);
    221   int_4 SetDiag();
    222 
    223   static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
    224   static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
    225 
    226   TRCKind Kind() const { return kind; }
    227   uint_4 NElts() const;
    228   T& operator()(uint_4 i);
    229   T  operator()(uint_4 i) const;
    230 
    231   TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
    232   TVector<T> GetVect() const;
    233 
    234   TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
    235   TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
    236 
    237   TMatrixRC<T>& operator *= (T x);
    238   TMatrixRC<T>& operator /= (T x);
    239   TMatrixRC<T>& operator -= (T x);
    240   TMatrixRC<T>& operator += (T x);
    241 
    242   TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
    243   TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
    244 
    245   uint_4 IMaxAbs(uint_4 first=0);
    246 
    247   static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
    248 
    249 protected:
    250   TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
    251   TMatrix<T>* matrix;
    252   inline static double Abs_Value(uint_1 v) {return (double) v;}
    253   inline static double Abs_Value(uint_2 v) {return (double) v;}
    254   inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
    255   inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
    256   inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
    257   inline static double Abs_Value(uint_4 v) {return (double) v;}
    258   inline static double Abs_Value(uint_8 v) {return (double) v;}
    259   inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
    260   inline static double Abs_Value(r_8 v)    {return fabs(v);}
    261   inline static double Abs_Value(complex<float> v)
    262                 {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
    263   inline static double Abs_Value(complex<double> v)
    264                 {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
    265 
    266   T*          data;
    267   int_4       index;
    268   uint_4      step;
    269   TRCKind     kind;
    270 };
    271 
    272 
    273 template <class T>
    274 inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
    275   {
    276   if ( a.NElts() != b.NElts() )
    277     throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
    278   if ( a.Kind() != b.Kind() )
    279     throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
    280   T sum = 0;
    281   for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
    282   return sum;
    283   }
    284 
    285 template <class T>
    286 inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
    287   { switch (rckind) { case TmatrixRow  : return 1;
    288                       case TmatrixCol  : return m.mNc;
    289                       case TmatrixDiag : return m.mNc+1; }
    290     return 0; }
    291 
    292 template <class T>
    293 inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
    294   { switch (rckind) { case TmatrixRow  : return const_cast<T *>(m.Data()) + index * m.mNc;
    295                       case TmatrixCol  : return const_cast<T *>(m.Data()) + index;
    296                       case TmatrixDiag : return const_cast<T *>(m.Data()); }
    297     return NULL; }
    298 
    299 template <class T> inline uint_4 TMatrixRC<T>::NElts() const
    300   { if (!matrix) return 0;
    301     switch (kind) { case TmatrixRow  : return matrix->mNc;
    302                     case TmatrixCol  : return matrix->mNr;
    303                     case TmatrixDiag : return matrix->mNc; }
    304     return 0; }
    305 
    306 template <class T>
    307 inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
    308 template <class T>
    309 inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
    310 
    311 ////////////////////////////////////////////////////////////////
    312 // Typedef pour simplifier et compatibilite Peida
    313 typedef TMatrixRC<r_8> MatrixRC;
    314191
    315192} // Fin du namespace
  • trunk/SophyaLib/TArray/tvector.h

    r762 r772  
    5555   {return TVector<T>(a * ((TMatrix<T> const&)(b)));}
    5656
    57 // Resolution du systeme A*C = B
    58 inline r_8 LinSolveInPlace(TMatrix<r_8>& a, TVector<r_8>& b)
    59 {
    60 if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
    61   throw(SzMismatchError("LinSolveInPlace(TMatrix<r_8>,TVector<r_8>) size mismatch"));
    62 return TMatrix<r_8>::GausPiv(a,b);
    63 }
    64 
    65 // Resolution du systeme A*C = B, avec C retourne dans B
    66 inline r_8 LinSolve(const TMatrix<r_8>& a, const TVector<r_8>& b, TVector<r_8>& c)
    67 {
    68 if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
    69   throw(SzMismatchError("LinSolve(TMatrix<r_8>,TVector<r_8>) size mismatch"));
    70 c = b;
    71 TMatrix<r_8> a1(a);
    72 return TMatrix<r_8>::GausPiv(a1,c);
    73 }
    7457
    7558// Typedef pour simplifier et compatibilite Peida
Note: See TracChangeset for help on using the changeset viewer.