source: Sophya/trunk/SophyaLib/NTools/tmatrix.h@ 501

Last change on this file since 501 was 501, checked in by ansari, 26 years ago

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