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

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

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