| 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 | class GeneralFit;
|
|---|
| 17 | template <class T> class TVector;
|
|---|
| 18 | template <class T> class TMatrixRC;
|
|---|
| 19 |
|
|---|
| 20 | template <class T>
|
|---|
| 21 | class TMatrix : public AnyDataObj {
|
|---|
| 22 | friend class TMatrixRC<T>;
|
|---|
| 23 | friend class TVector<T>;
|
|---|
| 24 | public:
|
|---|
| 25 |
|
|---|
| 26 | // Creation / destruction
|
|---|
| 27 | TMatrix();
|
|---|
| 28 | TMatrix(uint_4 r,uint_4 c);
|
|---|
| 29 | TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br=NULL);
|
|---|
| 30 | TMatrix(const TMatrix<T>& a);
|
|---|
| 31 | TMatrix(const TMatrix<T>& a,bool share);
|
|---|
| 32 | virtual ~TMatrix();
|
|---|
| 33 |
|
|---|
| 34 | // Temporaire?
|
|---|
| 35 | inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
|
|---|
| 36 | inline void SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
|
|---|
| 37 |
|
|---|
| 38 | // Gestion taille/Remplissage
|
|---|
| 39 | inline void Clone(const TMatrix<T>& a) // Clone: copie des donnees de "a"
|
|---|
| 40 | {mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;}
|
|---|
| 41 | inline void Reset(T v=0) {mNDBlock.Reset(v);}
|
|---|
| 42 | inline void ReSize(uint_4 r,uint_4 c) // Reallocation de place
|
|---|
| 43 | {if(r==0||c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
|
|---|
| 44 | mNDBlock.ReSize(r*c); mNr = r; mNc = c;}
|
|---|
| 45 | inline void Realloc(uint_4 r,uint_4 c,bool force=false)
|
|---|
| 46 | {if(r==0||c==0) throw(SzMismatchError("TMatrix::Realloc r ou c==0\n"));
|
|---|
| 47 | mNDBlock.Realloc(r*c,force); mNr = r; mNc = c;}
|
|---|
| 48 |
|
|---|
| 49 | // Informations pointeur/data
|
|---|
| 50 | inline uint_4 NRows() const {return mNr;}
|
|---|
| 51 | inline uint_4 NCols() const {return mNc;}
|
|---|
| 52 | inline uint_4 NCol() const {return mNc;} // back-compat Peida
|
|---|
| 53 | inline T const& operator()(uint_4 r,uint_4 c) const
|
|---|
| 54 | {return *(mNDBlock.Begin()+r*mNc+c);}
|
|---|
| 55 | inline T& operator()(uint_4 r,uint_4 c)
|
|---|
| 56 | {return *(mNDBlock.Begin()+r*mNc+c);}
|
|---|
| 57 | inline T const& operator[](uint_4 ip) const
|
|---|
| 58 | {return *(mNDBlock.Begin()+ip);}
|
|---|
| 59 | inline T& operator[](uint_4 ip)
|
|---|
| 60 | {return *(mNDBlock.Begin()+ip);}
|
|---|
| 61 | inline T* Data() {return mNDBlock.Begin();}
|
|---|
| 62 | inline const T* Data() const {return mNDBlock.Begin();}
|
|---|
| 63 | inline NDataBlock<T>& DataBlock() {return mNDBlock;}
|
|---|
| 64 | inline const NDataBlock<T>& DataBlock() const {return mNDBlock;}
|
|---|
| 65 |
|
|---|
| 66 | // Operations matricielles
|
|---|
| 67 | TMatrix<T> Transpose(void) const;
|
|---|
| 68 |
|
|---|
| 69 | // Operateur d'affectation
|
|---|
| 70 | // A = x (matrice diagonale x*Identite)
|
|---|
| 71 | inline TMatrix<T>& operator = (T x)
|
|---|
| 72 | {if(mNr!=mNc || mNr==0) throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0\n"));
|
|---|
| 73 | for(uint_4 r=0;r<mNr;r++) for(uint_4 c=0;c<mNc;c++) (*this)(r,c)=(r==c)?x:0;
|
|---|
| 74 | return *this;}
|
|---|
| 75 | // A = B : partage les donnees si "a" est temporaire, clone sinon.
|
|---|
| 76 | inline TMatrix<T>& operator = (const TMatrix<T>& a)
|
|---|
| 77 | {if(this == &a) return *this; CloneOrShare(a); return *this;}
|
|---|
| 78 |
|
|---|
| 79 | // Impression
|
|---|
| 80 | void Print(ostream& os,int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const;
|
|---|
| 81 | inline void Print(int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const
|
|---|
| 82 | {Print(cout,lp,i0,ni,j0,nj);}
|
|---|
| 83 |
|
|---|
| 84 | // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) (T) x
|
|---|
| 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 | inline TMatrix<T>& operator /= (T b) {mNDBlock /= b; return *this;}
|
|---|
| 89 |
|
|---|
| 90 | // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) B
|
|---|
| 91 | inline TMatrix<T>& operator += (const TMatrix<T>& a)
|
|---|
| 92 | {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
|
|---|
| 93 | throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
|
|---|
| 94 | mNDBlock += a.mNDBlock; return *this;}
|
|---|
| 95 | inline TMatrix<T>& operator -= (const TMatrix<T>& a)
|
|---|
| 96 | {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
|
|---|
| 97 | throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
|
|---|
| 98 | mNDBlock -= a.mNDBlock; return *this;}
|
|---|
| 99 | TMatrix<T>& operator *= (const TMatrix<T>& a);
|
|---|
| 100 |
|
|---|
| 101 | // Pour surcharge d'operateurs C = A (+,-,*) B
|
|---|
| 102 | TMatrix<T> Add(const TMatrix<T>& b) const;
|
|---|
| 103 | TMatrix<T> Sub(const TMatrix<T>& b) const;
|
|---|
| 104 | TMatrix<T> Mul(const TMatrix<T>& b) const;
|
|---|
| 105 |
|
|---|
| 106 | // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
|
|---|
| 107 | // operations sur la matrice B
|
|---|
| 108 | TMatrix<T> Inverse() const;
|
|---|
| 109 | static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
|
|---|
| 110 |
|
|---|
| 111 | // Residus et fonction fittees.
|
|---|
| 112 | TMatrix<T> FitResidus(GeneralFit& gfit
|
|---|
| 113 | ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);
|
|---|
| 114 | TMatrix<T> FitFunction(GeneralFit& gfit
|
|---|
| 115 | ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);
|
|---|
| 116 |
|
|---|
| 117 | // Acces aux rangees et colonnes
|
|---|
| 118 | TMatrixRC<T> Row(uint_4 r) const;
|
|---|
| 119 | TMatrixRC<T> Col(uint_4 c) const;
|
|---|
| 120 | TMatrixRC<T> Diag() const;
|
|---|
| 121 |
|
|---|
| 122 | protected:
|
|---|
| 123 | // partage les donnees si "a" temporaire, clone sinon.
|
|---|
| 124 | inline void CloneOrShare(const TMatrix<T>& a)
|
|---|
| 125 | {mNDBlock.CloneOrShare(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
|
|---|
| 126 | // Share: partage les donnees de "a"
|
|---|
| 127 | inline void Share(const TMatrix<T>& a)
|
|---|
| 128 | {mNDBlock.Share(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
|
|---|
| 129 |
|
|---|
| 130 | uint_4 mNr,mNc;
|
|---|
| 131 | NDataBlock<T> mNDBlock;
|
|---|
| 132 | };
|
|---|
| 133 |
|
|---|
| 134 | ////////////////////////////////////////////////////////////////
|
|---|
| 135 | // Impression
|
|---|
| 136 |
|
|---|
| 137 | template <class T>
|
|---|
| 138 | inline ostream& operator << (ostream& os, const TMatrix<T>& a)
|
|---|
| 139 | {a.Print(os); return(os);}
|
|---|
| 140 |
|
|---|
| 141 | ////////////////////////////////////////////////////////////////
|
|---|
| 142 | // Surcharge d'operateurs A (+,-,*,/) (T) x
|
|---|
| 143 |
|
|---|
| 144 | template <class T> inline TMatrix<T> operator + (const TMatrix<T>& a, T b)
|
|---|
| 145 | {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
|
|---|
| 146 |
|
|---|
| 147 | template <class T> inline TMatrix<T> operator + (T b,const TMatrix<T>& a)
|
|---|
| 148 | {TMatrix<T> result(a); result.SetTemp(true); result += b; 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);
|
|---|
| 155 | result.DataBlock() = b-result.DataBlock(); return result;}
|
|---|
| 156 |
|
|---|
| 157 | template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, T b)
|
|---|
| 158 | {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
|
|---|
| 159 |
|
|---|
| 160 | template <class T> inline TMatrix<T> operator * (T b,const TMatrix<T>& a)
|
|---|
| 161 | {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
|
|---|
| 162 |
|
|---|
| 163 | template <class T> inline TMatrix<T> operator / (const TMatrix<T>& a, T b)
|
|---|
| 164 | {TMatrix<T> result(a); result.SetTemp(true); result /= b; return result;}
|
|---|
| 165 |
|
|---|
| 166 | ////////////////////////////////////////////////////////////////
|
|---|
| 167 | // Surcharge d'operateurs C = A (+,-,*,/) B
|
|---|
| 168 |
|
|---|
| 169 | template <class T>
|
|---|
| 170 | inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
|
|---|
| 171 | {return a.Add(b);}
|
|---|
| 172 |
|
|---|
| 173 | template <class T>
|
|---|
| 174 | inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
|
|---|
| 175 | {return a.Sub(b);}
|
|---|
| 176 |
|
|---|
| 177 | template <class T>
|
|---|
| 178 | inline TMatrix<T> operator * (const TMatrix<T>& a,const TMatrix<T>& b)
|
|---|
| 179 | {return a.Mul(b);}
|
|---|
| 180 |
|
|---|
| 181 | ////////////////////////////////////////////////////////////////
|
|---|
| 182 | // Typedef pour simplifier et compatibilite Peida
|
|---|
| 183 | typedef TMatrix<r_8> Matrix;
|
|---|
| 184 |
|
|---|
| 185 | /////////////////////////////////////////////////////////////////////////
|
|---|
| 186 | // Classe pour la gestion de persistance
|
|---|
| 187 | template <class T>
|
|---|
| 188 | class FIO_TMatrix : public PPersist {
|
|---|
| 189 | public:
|
|---|
| 190 | FIO_TMatrix();
|
|---|
| 191 | FIO_TMatrix(string const & filename);
|
|---|
| 192 | FIO_TMatrix(const TMatrix<T> & obj);
|
|---|
| 193 | FIO_TMatrix(TMatrix<T> * obj);
|
|---|
| 194 | virtual ~FIO_TMatrix();
|
|---|
| 195 | virtual AnyDataObj* DataObj();
|
|---|
| 196 | inline operator TMatrix<T>() { return(*dobj); }
|
|---|
| 197 | protected :
|
|---|
| 198 | virtual void ReadSelf(PInPersist&);
|
|---|
| 199 | virtual void WriteSelf(POutPersist&) const;
|
|---|
| 200 | TMatrix<T> * dobj;
|
|---|
| 201 | bool ownobj;
|
|---|
| 202 | };
|
|---|
| 203 |
|
|---|
| 204 | template <class T>
|
|---|
| 205 | inline POutPersist& operator << (POutPersist& os, TMatrix<T> & obj)
|
|---|
| 206 | { FIO_TMatrix<T> fio(&obj); fio.Write(os); return(os); }
|
|---|
| 207 | template <class T>
|
|---|
| 208 | inline PInPersist& operator >> (PInPersist& is, TMatrix<T> & obj)
|
|---|
| 209 | { FIO_TMatrix<T> fio(&obj); fio.Read(is); return(is); }
|
|---|
| 210 |
|
|---|
| 211 | /////////////////////////////////////////////////////////////////////////
|
|---|
| 212 | // Classe de lignes/colonnes de matrices
|
|---|
| 213 | enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
|
|---|
| 214 | template <class T>
|
|---|
| 215 | class TMatrixRC {
|
|---|
| 216 | friend class TVector<T>;
|
|---|
| 217 | friend class TMatrix<T>;
|
|---|
| 218 | public:
|
|---|
| 219 | TMatrixRC();
|
|---|
| 220 |
|
|---|
| 221 | virtual ~TMatrixRC() {}
|
|---|
| 222 |
|
|---|
| 223 | int_4 Next();
|
|---|
| 224 | int_4 Prev();
|
|---|
| 225 | int_4 SetCol(int_4 c);
|
|---|
| 226 | int_4 SetRow(int_4 r);
|
|---|
| 227 | int_4 SetDiag();
|
|---|
| 228 |
|
|---|
| 229 | static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
|
|---|
| 230 | static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
|
|---|
| 231 |
|
|---|
| 232 | TRCKind Kind() const { return kind; }
|
|---|
| 233 | uint_4 NElts() const;
|
|---|
| 234 | T& operator()(uint_4 i);
|
|---|
| 235 | T operator()(uint_4 i) const;
|
|---|
| 236 |
|
|---|
| 237 | TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
|
|---|
| 238 | TVector<T> GetVect() const;
|
|---|
| 239 |
|
|---|
| 240 | TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
|
|---|
| 241 | TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
|
|---|
| 242 |
|
|---|
| 243 | TMatrixRC<T>& operator *= (T x);
|
|---|
| 244 | TMatrixRC<T>& operator /= (T x);
|
|---|
| 245 | TMatrixRC<T>& operator -= (T x);
|
|---|
| 246 | TMatrixRC<T>& operator += (T x);
|
|---|
| 247 |
|
|---|
| 248 | TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
|
|---|
| 249 | TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
|
|---|
| 250 |
|
|---|
| 251 | uint_4 IMaxAbs(uint_4 first=0);
|
|---|
| 252 |
|
|---|
| 253 | static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
|
|---|
| 254 |
|
|---|
| 255 | protected:
|
|---|
| 256 | TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
|
|---|
| 257 | TMatrix<T>* matrix;
|
|---|
| 258 | inline static double Abs_Value(uint_1 v) {return (double) v;}
|
|---|
| 259 | inline static double Abs_Value(uint_2 v) {return (double) v;}
|
|---|
| 260 | inline static double Abs_Value(int_2 v) {return (v>0)? (double) v: (double) -v;}
|
|---|
| 261 | inline static double Abs_Value(int_4 v) {return (v>0)? (double) v: (double) -v;}
|
|---|
| 262 | inline static double Abs_Value(int_8 v) {return (v>0)? (double) v: (double) -v;}
|
|---|
| 263 | inline static double Abs_Value(uint_4 v) {return (double) v;}
|
|---|
| 264 | inline static double Abs_Value(uint_8 v) {return (double) v;}
|
|---|
| 265 | inline static double Abs_Value(r_4 v) {return (double) fabsf(v);}
|
|---|
| 266 | inline static double Abs_Value(r_8 v) {return fabs(v);}
|
|---|
| 267 | inline static double Abs_Value(complex<float> v)
|
|---|
| 268 | {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
|
|---|
| 269 | inline static double Abs_Value(complex<double> v)
|
|---|
| 270 | {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
|
|---|
| 271 |
|
|---|
| 272 | T* data;
|
|---|
| 273 | int_4 index;
|
|---|
| 274 | uint_4 step;
|
|---|
| 275 | TRCKind kind;
|
|---|
| 276 | };
|
|---|
| 277 |
|
|---|
| 278 |
|
|---|
| 279 | template <class T>
|
|---|
| 280 | inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
|
|---|
| 281 | {
|
|---|
| 282 | if ( a.NElts() != b.NElts() )
|
|---|
| 283 | throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
|
|---|
| 284 | if ( a.Kind() != b.Kind() )
|
|---|
| 285 | throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
|
|---|
| 286 | T sum = 0;
|
|---|
| 287 | for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
|
|---|
| 288 | return sum;
|
|---|
| 289 | }
|
|---|
| 290 |
|
|---|
| 291 | template <class T>
|
|---|
| 292 | inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
|
|---|
| 293 | { switch (rckind) { case TmatrixRow : return 1;
|
|---|
| 294 | case TmatrixCol : return m.mNc;
|
|---|
| 295 | case TmatrixDiag : return m.mNc+1; }
|
|---|
| 296 | return 0; }
|
|---|
| 297 |
|
|---|
| 298 | template <class T>
|
|---|
| 299 | inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
|
|---|
| 300 | { switch (rckind) { case TmatrixRow : return const_cast<T *>(m.Data()) + index * m.mNc;
|
|---|
| 301 | case TmatrixCol : return const_cast<T *>(m.Data()) + index;
|
|---|
| 302 | case TmatrixDiag : return const_cast<T *>(m.Data()); }
|
|---|
| 303 | return NULL; }
|
|---|
| 304 |
|
|---|
| 305 | template <class T> inline uint_4 TMatrixRC<T>::NElts() const
|
|---|
| 306 | { if (!matrix) return 0;
|
|---|
| 307 | switch (kind) { case TmatrixRow : return matrix->mNc;
|
|---|
| 308 | case TmatrixCol : return matrix->mNr;
|
|---|
| 309 | case TmatrixDiag : return matrix->mNc; }
|
|---|
| 310 | return 0; }
|
|---|
| 311 |
|
|---|
| 312 | template <class T>
|
|---|
| 313 | inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
|
|---|
| 314 | template <class T>
|
|---|
| 315 | inline T TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
|
|---|
| 316 |
|
|---|
| 317 | ////////////////////////////////////////////////////////////////
|
|---|
| 318 | // Typedef pour simplifier et compatibilite Peida
|
|---|
| 319 | typedef TMatrixRC<r_8> MatrixRC;
|
|---|
| 320 |
|
|---|
| 321 | } // Fin du namespace
|
|---|
| 322 |
|
|---|
| 323 | #endif
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
