Changeset 772 in Sophya for trunk/SophyaLib/TArray
- Timestamp:
- Mar 10, 2000, 4:13:22 PM (26 years ago)
- Location:
- trunk/SophyaLib/TArray
- Files:
-
- 6 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/TArray/tarrinit.cc
r762 r772 5 5 #include "tmatrix.h" 6 6 #include "tvector.h" 7 #include "fioarr.h" 7 8 8 9 … … 17 18 // Enregistrement des classes PPersist du modules TArray 18 19 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> >); 19 31 20 32 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:24ansari Exp $1 // $Id: tmatrix.cc,v 1.2 2000-03-10 15:13:22 ansari Exp $ 2 2 // C.Magneville 04/99 3 3 #include "machdefs.h" … … 161 161 } 162 162 163 ////////////////////////////////////////////////////////////////164 //**** Pour gestion des lignes et des colonnes165 166 template <class T> TMatrixRC<T> TMatrix<T>::Row(uint_4 r) const167 {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) const173 {174 TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixCol, c);175 return rc;176 }177 178 template <class T> TMatrixRC<T> TMatrix<T>::Diag() const179 {180 TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixDiag);181 return rc;182 }183 184 ////////////////////////////////////////////////////////////////185 //**** Pour inversion186 #ifndef M_LN2187 #define M_LN2 0.69314718055994530942188 #endif189 #ifndef LN_MINDOUBLE190 #define LN_MINDOUBLE (M_LN2 * (DMINEXP - 1))191 #endif192 #ifndef LN_MAXDOUBLE193 #define LN_MAXDOUBLE (M_LN2 * DMAXEXP)194 #endif195 196 r_8 TMatrix<r_8>::GausPiv(TMatrix<r_8>& a, TMatrix<r_8>& b)197 // Pivot de Gauss198 // * Attention: egcs impose que cette fonction soit mise dans le .cc199 // 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 constructors248 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 * pivRowa252 b.Row(i).LinComb(r, pivRowb);253 }254 }255 det *= a(n-1, n-1);256 257 // on remonte258 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 constructors262 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() const279 // Inversion280 {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 287 163 288 164 /////////////////////////////////////////////////////////// … … 379 255 #include "tvector.h" 380 256 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() const463 {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 }561 257 562 258 /////////////////////////////////////////////////////////////// … … 585 281 #pragma define_template FIO_TMatrix< complex<r_4> > 586 282 #pragma define_template FIO_TMatrix< complex<r_8> > 587 // Instances gestion lignes/colonnes588 #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> >599 283 #endif 600 284 … … 623 307 template class FIO_TMatrix< complex<r_4> >; 624 308 template class FIO_TMatrix< complex<r_8> >; 625 // Instances gestion lignes/colonnes626 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> >;637 309 #endif -
trunk/SophyaLib/TArray/tmatrix.h
r762 r772 14 14 namespace SOPHYA { 15 15 16 template <class T> class TVector;17 template <class T> class TMatrixRC;18 19 16 template <class T> 20 17 class TMatrix : public AnyDataObj { 21 friend class TMatrixRC<T>;22 friend class TVector<T>;23 18 public: 24 19 … … 103 98 TMatrix<T> Mul(const TMatrix<T>& b) const; 104 99 105 // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes106 // operations sur la matrice B107 TMatrix<T> Inverse() const;108 static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);109 110 // Acces aux rangees et colonnes111 TMatrixRC<T> Row(uint_4 r) const;112 TMatrixRC<T> Col(uint_4 c) const;113 TMatrixRC<T> Diag() const;114 100 115 101 protected: … … 203 189 { FIO_TMatrix<T> fio(&obj); fio.Read(is); return(is); } 204 190 205 /////////////////////////////////////////////////////////////////////////206 // Classe de lignes/colonnes de matrices207 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() const300 { 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 Peida313 typedef TMatrixRC<r_8> MatrixRC;314 191 315 192 } // Fin du namespace -
trunk/SophyaLib/TArray/tvector.h
r762 r772 55 55 {return TVector<T>(a * ((TMatrix<T> const&)(b)));} 56 56 57 // Resolution du systeme A*C = B58 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 B66 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 }74 57 75 58 // Typedef pour simplifier et compatibilite Peida
Note:
See TracChangeset
for help on using the changeset viewer.