| 1 | // This may look like C code, but it is really -*- C++ -*-
|
|---|
| 2 | // C.Magneville 04/99
|
|---|
| 3 | #ifndef TMatrix_SEEN
|
|---|
| 4 | #define TMatrix_SEEN
|
|---|
| 5 |
|
|---|
| 6 | #include "machdefs.h"
|
|---|
| 7 | #include <stdio.h>
|
|---|
| 8 | #include <iostream.h>
|
|---|
| 9 | #include <complex>
|
|---|
| 10 | #include "ppersist.h"
|
|---|
| 11 | #include "anydataobj.h"
|
|---|
| 12 | #include "ndatablock.h"
|
|---|
| 13 |
|
|---|
| 14 | namespace SOPHYA {
|
|---|
| 15 |
|
|---|
| 16 | template <class T> class TVector;
|
|---|
| 17 | template <class T> class TMatrixRC;
|
|---|
| 18 |
|
|---|
| 19 | template <class T>
|
|---|
| 20 | class TMatrix : public AnyDataObj {
|
|---|
| 21 | friend class TMatrixRC<T>;
|
|---|
| 22 | friend class TVector<T>;
|
|---|
| 23 | public:
|
|---|
| 24 |
|
|---|
| 25 | // Creation / destruction
|
|---|
| 26 | TMatrix();
|
|---|
| 27 | TMatrix(uint_4 r,uint_4 c);
|
|---|
| 28 | TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br=NULL);
|
|---|
| 29 | TMatrix(const TMatrix<T>& a);
|
|---|
| 30 | TMatrix(const TMatrix<T>& a,bool share);
|
|---|
| 31 | virtual ~TMatrix();
|
|---|
| 32 |
|
|---|
| 33 | // Temporaire?
|
|---|
| 34 | inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
|
|---|
| 35 | inline void SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
|
|---|
| 36 |
|
|---|
| 37 | // Gestion taille/Remplissage
|
|---|
| 38 | inline void Clone(const TMatrix<T>& a) // Clone: copie des donnees de "a"
|
|---|
| 39 | {mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;}
|
|---|
| 40 | inline void Reset(T v=0) {mNDBlock.Reset(v);}
|
|---|
| 41 | inline void ReSize(uint_4 r,uint_4 c) // Reallocation de place
|
|---|
| 42 | {if(r==0||c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
|
|---|
| 43 | mNDBlock.ReSize(r*c); mNr = r; mNc = c;}
|
|---|
| 44 | inline void Realloc(uint_4 r,uint_4 c,bool force=false)
|
|---|
| 45 | {if(r==0||c==0) throw(SzMismatchError("TMatrix::Realloc r ou c==0\n"));
|
|---|
| 46 | mNDBlock.Realloc(r*c,force); mNr = r; mNc = c;}
|
|---|
| 47 |
|
|---|
| 48 | // Informations pointeur/data
|
|---|
| 49 | inline uint_4 NRows() const {return mNr;}
|
|---|
| 50 | inline uint_4 NCols() const {return mNc;}
|
|---|
| 51 | inline uint_4 NCol() const {return mNc;} // back-compat Peida
|
|---|
| 52 | inline T const& operator()(uint_4 r,uint_4 c) const
|
|---|
| 53 | {return *(mNDBlock.Begin()+r*mNc+c);}
|
|---|
| 54 | inline T& operator()(uint_4 r,uint_4 c)
|
|---|
| 55 | {return *(mNDBlock.Begin()+r*mNc+c);}
|
|---|
| 56 | inline T const& operator[](uint_4 ip) const
|
|---|
| 57 | {return *(mNDBlock.Begin()+ip);}
|
|---|
| 58 | inline T& operator[](uint_4 ip)
|
|---|
| 59 | {return *(mNDBlock.Begin()+ip);}
|
|---|
| 60 | inline T* Data() {return mNDBlock.Begin();}
|
|---|
| 61 | inline const T* Data() const {return mNDBlock.Begin();}
|
|---|
| 62 | inline NDataBlock<T>& DataBlock() {return mNDBlock;}
|
|---|
| 63 | inline const NDataBlock<T>& DataBlock() const {return mNDBlock;}
|
|---|
| 64 |
|
|---|
| 65 | // Operations matricielles
|
|---|
| 66 | TMatrix<T> Transpose(void) const;
|
|---|
| 67 |
|
|---|
| 68 | // Operateur d'affectation
|
|---|
| 69 | // A = x (matrice diagonale x*Identite)
|
|---|
| 70 | inline TMatrix<T>& operator = (T x)
|
|---|
| 71 | {if(mNr!=mNc || mNr==0) throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0\n"));
|
|---|
| 72 | for(uint_4 r=0;r<mNr;r++) for(uint_4 c=0;c<mNc;c++) (*this)(r,c)=(r==c)?x:(T)0;
|
|---|
| 73 | return *this;}
|
|---|
| 74 | // A = B : partage les donnees si "a" est temporaire, clone sinon.
|
|---|
| 75 | inline TMatrix<T>& operator = (const TMatrix<T>& a)
|
|---|
| 76 | {if(this == &a) return *this; CloneOrShare(a); return *this;}
|
|---|
| 77 |
|
|---|
| 78 | // Impression
|
|---|
| 79 | void Print(ostream& os,int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const;
|
|---|
| 80 | inline void Print(int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const
|
|---|
| 81 | {Print(cout,lp,i0,ni,j0,nj);}
|
|---|
| 82 |
|
|---|
| 83 | // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) (T) x
|
|---|
| 84 | inline TMatrix<T>& operator += (T b) {mNDBlock += b; return *this;}
|
|---|
| 85 | inline TMatrix<T>& operator -= (T b) {mNDBlock -= b; return *this;}
|
|---|
| 86 | inline TMatrix<T>& operator *= (T b) {mNDBlock *= b; return *this;}
|
|---|
| 87 | inline TMatrix<T>& operator /= (T b) {mNDBlock /= b; return *this;}
|
|---|
| 88 |
|
|---|
| 89 | // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) B
|
|---|
| 90 | inline TMatrix<T>& operator += (const TMatrix<T>& a)
|
|---|
| 91 | {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
|
|---|
| 92 | throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
|
|---|
| 93 | mNDBlock += a.mNDBlock; return *this;}
|
|---|
| 94 | inline TMatrix<T>& operator -= (const TMatrix<T>& a)
|
|---|
| 95 | {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
|
|---|
| 96 | throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
|
|---|
| 97 | mNDBlock -= a.mNDBlock; return *this;}
|
|---|
| 98 | TMatrix<T>& operator *= (const TMatrix<T>& a);
|
|---|
| 99 |
|
|---|
| 100 | // Pour surcharge d'operateurs C = A (+,-,*) B
|
|---|
| 101 | TMatrix<T> Add(const TMatrix<T>& b) const;
|
|---|
| 102 | TMatrix<T> Sub(const TMatrix<T>& b) const;
|
|---|
| 103 | TMatrix<T> Mul(const TMatrix<T>& b) const;
|
|---|
| 104 |
|
|---|
| 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;
|
|---|
| 114 |
|
|---|
| 115 | protected:
|
|---|
| 116 | // partage les donnees si "a" temporaire, clone sinon.
|
|---|
| 117 | inline void CloneOrShare(const TMatrix<T>& a)
|
|---|
| 118 | {mNDBlock.CloneOrShare(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
|
|---|
| 119 | // Share: partage les donnees de "a"
|
|---|
| 120 | inline void Share(const TMatrix<T>& a)
|
|---|
| 121 | {mNDBlock.Share(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
|
|---|
| 122 |
|
|---|
| 123 | uint_4 mNr,mNc;
|
|---|
| 124 | NDataBlock<T> mNDBlock;
|
|---|
| 125 | };
|
|---|
| 126 |
|
|---|
| 127 | ////////////////////////////////////////////////////////////////
|
|---|
| 128 | // Impression
|
|---|
| 129 |
|
|---|
| 130 | template <class T>
|
|---|
| 131 | inline ostream& operator << (ostream& os, const TMatrix<T>& a)
|
|---|
| 132 | {a.Print(os); return(os);}
|
|---|
| 133 |
|
|---|
| 134 | ////////////////////////////////////////////////////////////////
|
|---|
| 135 | // Surcharge d'operateurs A (+,-,*,/) (T) x
|
|---|
| 136 |
|
|---|
| 137 | template <class T> inline TMatrix<T> operator + (const TMatrix<T>& a, T b)
|
|---|
| 138 | {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
|
|---|
| 139 |
|
|---|
| 140 | template <class T> inline TMatrix<T> operator + (T b,const TMatrix<T>& a)
|
|---|
| 141 | {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
|
|---|
| 142 |
|
|---|
| 143 | template <class T> inline TMatrix<T> operator - (const TMatrix<T>& a, T b)
|
|---|
| 144 | {TMatrix<T> result(a); result.SetTemp(true); result -= b; return result;}
|
|---|
| 145 |
|
|---|
| 146 | template <class T> inline TMatrix<T> operator - (T b,const TMatrix<T>& a)
|
|---|
| 147 | {TMatrix<T> result(a); result.SetTemp(true);
|
|---|
| 148 | result.DataBlock() = b-result.DataBlock(); return result;}
|
|---|
| 149 |
|
|---|
| 150 | template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, T b)
|
|---|
| 151 | {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
|
|---|
| 152 |
|
|---|
| 153 | template <class T> inline TMatrix<T> operator * (T b,const TMatrix<T>& a)
|
|---|
| 154 | {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
|
|---|
| 155 |
|
|---|
| 156 | template <class T> inline TMatrix<T> operator / (const TMatrix<T>& a, T b)
|
|---|
| 157 | {TMatrix<T> result(a); result.SetTemp(true); result /= b; return result;}
|
|---|
| 158 |
|
|---|
| 159 | ////////////////////////////////////////////////////////////////
|
|---|
| 160 | // Surcharge d'operateurs C = A (+,-,*,/) B
|
|---|
| 161 |
|
|---|
| 162 | template <class T>
|
|---|
| 163 | inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
|
|---|
| 164 | {return a.Add(b);}
|
|---|
| 165 |
|
|---|
| 166 | template <class T>
|
|---|
| 167 | inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
|
|---|
| 168 | {return a.Sub(b);}
|
|---|
| 169 |
|
|---|
| 170 | template <class T>
|
|---|
| 171 | inline TMatrix<T> operator * (const TMatrix<T>& a,const TMatrix<T>& b)
|
|---|
| 172 | {return a.Mul(b);}
|
|---|
| 173 |
|
|---|
| 174 | ////////////////////////////////////////////////////////////////
|
|---|
| 175 | // Typedef pour simplifier et compatibilite Peida
|
|---|
| 176 | typedef TMatrix<r_8> Matrix;
|
|---|
| 177 |
|
|---|
| 178 | /////////////////////////////////////////////////////////////////////////
|
|---|
| 179 | // Classe pour la gestion de persistance
|
|---|
| 180 | template <class T>
|
|---|
| 181 | class FIO_TMatrix : public PPersist {
|
|---|
| 182 | public:
|
|---|
| 183 | FIO_TMatrix();
|
|---|
| 184 | FIO_TMatrix(string const & filename);
|
|---|
| 185 | FIO_TMatrix(const TMatrix<T> & obj);
|
|---|
| 186 | FIO_TMatrix(TMatrix<T> * obj);
|
|---|
| 187 | virtual ~FIO_TMatrix();
|
|---|
| 188 | virtual AnyDataObj* DataObj();
|
|---|
| 189 | virtual void SetDataObj(AnyDataObj & o);
|
|---|
| 190 | inline operator TMatrix<T>() { return(*dobj); }
|
|---|
| 191 | protected :
|
|---|
| 192 | virtual void ReadSelf(PInPersist&);
|
|---|
| 193 | virtual void WriteSelf(POutPersist&) const;
|
|---|
| 194 | TMatrix<T> * dobj;
|
|---|
| 195 | bool ownobj;
|
|---|
| 196 | };
|
|---|
| 197 |
|
|---|
| 198 | template <class T>
|
|---|
| 199 | inline POutPersist& operator << (POutPersist& os, TMatrix<T> & obj)
|
|---|
| 200 | { FIO_TMatrix<T> fio(&obj); fio.Write(os); return(os); }
|
|---|
| 201 | template <class T>
|
|---|
| 202 | inline PInPersist& operator >> (PInPersist& is, TMatrix<T> & obj)
|
|---|
| 203 | { FIO_TMatrix<T> fio(&obj); fio.Read(is); return(is); }
|
|---|
| 204 |
|
|---|
| 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;
|
|---|
| 314 |
|
|---|
| 315 | } // Fin du namespace
|
|---|
| 316 |
|
|---|
| 317 | #endif
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
