| 1 | // $Id: tmatrix.cc,v 1.4 1999-05-05 05:51:54 ansari Exp $
|
|---|
| 2 | // C.Magneville 04/99
|
|---|
| 3 | #include "machdefs.h"
|
|---|
| 4 | #include <stdio.h>
|
|---|
| 5 | #include <stdlib.h>
|
|---|
| 6 | #include <iostream.h>
|
|---|
| 7 | #include <complex>
|
|---|
| 8 | #include "pexceptions.h"
|
|---|
| 9 | #include "tmatrix.h"
|
|---|
| 10 | #include "objfio.h"
|
|---|
| 11 |
|
|---|
| 12 | using namespace PlanckDPC;
|
|---|
| 13 |
|
|---|
| 14 | ////////////////////////////////////////////////////////////////
|
|---|
| 15 | //**** Createur, Destructeur
|
|---|
| 16 |
|
|---|
| 17 | template <class T>
|
|---|
| 18 | TMatrix<T>::TMatrix()
|
|---|
| 19 | // Constructeur par defaut.
|
|---|
| 20 | : mNr(0), mNc(0), mNDBlock()
|
|---|
| 21 | {
|
|---|
| 22 | }
|
|---|
| 23 |
|
|---|
| 24 | template <class T>
|
|---|
| 25 | TMatrix<T>::TMatrix(uint_4 r,uint_4 c)
|
|---|
| 26 | // Construit une matrice de r lignes et c colonnes.
|
|---|
| 27 | : mNr(r), mNc(c), mNDBlock(r*c)
|
|---|
| 28 | {
|
|---|
| 29 | }
|
|---|
| 30 |
|
|---|
| 31 | template <class T>
|
|---|
| 32 | TMatrix<T>::TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br)
|
|---|
| 33 | // Construit une matrice de r lignes et c colonnes. On fournit
|
|---|
| 34 | // le tableau des valeurs et eventuellement un Bridge.
|
|---|
| 35 | : mNr(r), mNc(c), mNDBlock(r*c,values,br)
|
|---|
| 36 | {
|
|---|
| 37 | }
|
|---|
| 38 |
|
|---|
| 39 | template <class T>
|
|---|
| 40 | TMatrix<T>::TMatrix(const TMatrix<T>& a)
|
|---|
| 41 | // Constructeur par copie (partage si "a" temporaire).
|
|---|
| 42 | : mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
|
|---|
| 43 | {
|
|---|
| 44 | }
|
|---|
| 45 |
|
|---|
| 46 | template <class T>
|
|---|
| 47 | TMatrix<T>::TMatrix(const TMatrix<T>& a,bool share)
|
|---|
| 48 | // Constructeur par copie avec possibilite de forcer le partage ou non.
|
|---|
| 49 | : mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
|
|---|
| 50 | {
|
|---|
| 51 | }
|
|---|
| 52 |
|
|---|
| 53 | template <class T>
|
|---|
| 54 | TMatrix<T>::~TMatrix()
|
|---|
| 55 | // Destructeur
|
|---|
| 56 | {
|
|---|
| 57 | }
|
|---|
| 58 |
|
|---|
| 59 | ////////////////////////////////////////////////////////////////
|
|---|
| 60 | // Operations matricielles
|
|---|
| 61 | template <class T>
|
|---|
| 62 | TMatrix<T> TMatrix<T>::Transpose(void) const
|
|---|
| 63 | // Transposition
|
|---|
| 64 | {
|
|---|
| 65 | TMatrix<T> a; a.Clone(*this); a.SetTemp(true);
|
|---|
| 66 | a.mNr = mNc; a.mNc = mNr;
|
|---|
| 67 | {for(int i=0; i<mNr; i++)
|
|---|
| 68 | for(int j=0; j<mNc; j++) {
|
|---|
| 69 | a(j,i) = (*this)(i,j);
|
|---|
| 70 | }}
|
|---|
| 71 | return a;
|
|---|
| 72 | }
|
|---|
| 73 |
|
|---|
| 74 | ////////////////////////////////////////////////////////////////
|
|---|
| 75 | //**** Impression
|
|---|
| 76 |
|
|---|
| 77 | template <class T>
|
|---|
| 78 | void TMatrix<T>::Print(ostream& os,int lp
|
|---|
| 79 | ,uint_4 i0,uint_4 ni,uint_4 j0,uint_4 nj) const
|
|---|
| 80 | // Impression de la sous-matrice (i:i+ni-1,i:j+nj-1)
|
|---|
| 81 | {
|
|---|
| 82 | os<<"TMatrix::Print("<<mNr<<","<<mNc<<")"<<endl;
|
|---|
| 83 | if(lp>0)
|
|---|
| 84 | {os<<" this="<<this<<endl; mNDBlock.Print(0,0);}
|
|---|
| 85 | if(mNr==0 || mNc==0) return;
|
|---|
| 86 | if(i0>=mNr || j0>=mNc || ni==0 || nj==0) return;
|
|---|
| 87 | uint_4 i1 = i0+ni; if(i1>mNr) i1 = mNr;
|
|---|
| 88 | uint_4 j1 = j0+nj; if(j1>mNc) j1 = mNc;
|
|---|
| 89 | for(uint_4 i=i0;i<i1;i++) {
|
|---|
| 90 | for(uint_4 j=j0;j<j1;j++) cout<<" "<<(*this)(i,j);
|
|---|
| 91 | cout<<endl;
|
|---|
| 92 | }
|
|---|
| 93 | }
|
|---|
| 94 |
|
|---|
| 95 | ////////////////////////////////////////////////////////////////
|
|---|
| 96 | //**** Surcharge de *= (INPLACE): TMatrix *= TMatrix;
|
|---|
| 97 |
|
|---|
| 98 | template <class T>
|
|---|
| 99 | TMatrix<T>& TMatrix<T>::operator *= (const TMatrix<T>& a)
|
|---|
| 100 | // A = A*B -> A(n,m) = A(n,m)*B(m,m)
|
|---|
| 101 | {
|
|---|
| 102 | uint_4 ndata = mNr*mNc;
|
|---|
| 103 | if(ndata==0 || mNc!=a.mNr || a.mNr!=a.mNc)
|
|---|
| 104 | throw(SzMismatchError("TMatrix::operator*=A size mismatch"));
|
|---|
| 105 | // A(i,j) = Sum(k) A(i,k)*B(k,j) ... il faut sauver la ligne "i" de A
|
|---|
| 106 | // Vecteur oi : vecteur ou est sauve la ligne i de la matrice *this
|
|---|
| 107 | // oi,oe = pointeur de debut et de fin du vecteur temporaire
|
|---|
| 108 | // oij = pointeur parcourant le vecteur oi
|
|---|
| 109 | // Matrice *this: i = pointeur du debut de la ligne i
|
|---|
| 110 | // ij = pointeur parcourant la ligne i
|
|---|
| 111 | // Matrice a : aj = pointeur de debut de la colonne j
|
|---|
| 112 | // aji = pointeur parcourant la colonne j
|
|---|
| 113 | T* oi = new T[mNc]; T* oe = oi+mNc;
|
|---|
| 114 | for(T *i=Data(); i<Data()+ndata; i+=mNc) {
|
|---|
| 115 | {for(T *oij=oi,*ij=i; oij<oe;) *oij++ = *ij++;}
|
|---|
| 116 | {for(T *ij=i,*aj=const_cast<T *>(a.Data()); aj<a.Data()+a.mNc; ij++,aj++) {
|
|---|
| 117 | T sum = 0;
|
|---|
| 118 | for(T *oij=oi,*aji=aj; oij<oe; oij++,aji+=a.mNc) sum += *oij * *aji;
|
|---|
| 119 | *ij = sum;
|
|---|
| 120 | }}
|
|---|
| 121 | }
|
|---|
| 122 | delete [] oi;
|
|---|
| 123 | return *this;
|
|---|
| 124 | }
|
|---|
| 125 |
|
|---|
| 126 | ////////////////////////////////////////////////////////////////
|
|---|
| 127 | //**** Pour surcharge d'operateurs C = A (+,-,*,/) B
|
|---|
| 128 |
|
|---|
| 129 | template <class T> TMatrix<T> TMatrix<T>::Add(const TMatrix<T>& b) const
|
|---|
| 130 | {
|
|---|
| 131 | if(mNr!=b.mNr || mNc!=b.mNc)
|
|---|
| 132 | throw(SzMismatchError("NDataBlock operator C=A+B size mismatch\n"));
|
|---|
| 133 | TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
|
|---|
| 134 | result.mNDBlock = mNDBlock+b.mNDBlock;
|
|---|
| 135 | return result;
|
|---|
| 136 | }
|
|---|
| 137 |
|
|---|
| 138 | template <class T> TMatrix<T> TMatrix<T>::Sub(const TMatrix<T>& b) const
|
|---|
| 139 | {
|
|---|
| 140 | if(mNr!=b.mNr || mNc!=b.mNc)
|
|---|
| 141 | throw(SzMismatchError("NDataBlock operator C=A-B size mismatch\n"));
|
|---|
| 142 | TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
|
|---|
| 143 | result.mNDBlock = mNDBlock-b.mNDBlock;
|
|---|
| 144 | return result;
|
|---|
| 145 | }
|
|---|
| 146 |
|
|---|
| 147 | template <class T> TMatrix<T> TMatrix<T>::Mul(const TMatrix<T>& b) const
|
|---|
| 148 | // C = A(this)*B : Cij = Aik Bkj (allocation forcee dans tous les cas)
|
|---|
| 149 | {
|
|---|
| 150 | if(mNr==0 || mNc==0 || b.mNc==0 || mNc!=b.mNr)
|
|---|
| 151 | throw(SzMismatchError("NDataBlock operator C=A*B size mismatch\n"));
|
|---|
| 152 | TMatrix<T> r; r.SetTemp(true); r.ReSize(mNr,b.mNc);
|
|---|
| 153 | T *ai,*aik,*bj,*bkj,*ri,*rij;
|
|---|
| 154 | for(ri=const_cast<T *>(r.Data()),ai=const_cast<T *>(Data());
|
|---|
| 155 | ri<r.Data()+r.mNr*r.mNc;ri+=r.mNc,ai+=mNc) {
|
|---|
| 156 | for(rij=ri,bj=const_cast<T *>(b.Data());rij<ri+r.mNc;rij++,bj++) {
|
|---|
| 157 | *rij = 0;
|
|---|
| 158 | for(aik=ai,bkj=bj;aik<ai+mNc;aik++,bkj+=b.mNc) *rij += *aik * *bkj;
|
|---|
| 159 | }
|
|---|
| 160 | }
|
|---|
| 161 | return r;
|
|---|
| 162 | }
|
|---|
| 163 |
|
|---|
| 164 | ////////////////////////////////////////////////////////////////
|
|---|
| 165 | // -------------------------------------------------------------------------
|
|---|
| 166 | // Les objets delegues pour la gestion de persistance
|
|---|
| 167 | // -------------------------------------------------------------------------
|
|---|
| 168 |
|
|---|
| 169 | template <class T>
|
|---|
| 170 | FIO_TMatrix<T>::FIO_TMatrix()
|
|---|
| 171 | {
|
|---|
| 172 | dobj=new TMatrix<T>;
|
|---|
| 173 | ownobj=true;
|
|---|
| 174 | }
|
|---|
| 175 |
|
|---|
| 176 | template <class T>
|
|---|
| 177 | FIO_TMatrix<T>::FIO_TMatrix(string const & filename)
|
|---|
| 178 | {
|
|---|
| 179 | dobj=new TMatrix<T>;
|
|---|
| 180 | ownobj=true;
|
|---|
| 181 | Read(filename);
|
|---|
| 182 | }
|
|---|
| 183 |
|
|---|
| 184 | template <class T>
|
|---|
| 185 | FIO_TMatrix<T>::FIO_TMatrix(const TMatrix<T> & obj)
|
|---|
| 186 | {
|
|---|
| 187 | dobj = new TMatrix<T>(obj);
|
|---|
| 188 | ownobj=true;
|
|---|
| 189 | }
|
|---|
| 190 |
|
|---|
| 191 | template <class T>
|
|---|
| 192 | FIO_TMatrix<T>::FIO_TMatrix(TMatrix<T> * obj)
|
|---|
| 193 | {
|
|---|
| 194 | dobj = obj;
|
|---|
| 195 | ownobj=false;
|
|---|
| 196 | }
|
|---|
| 197 |
|
|---|
| 198 | template <class T>
|
|---|
| 199 | FIO_TMatrix<T>::~FIO_TMatrix()
|
|---|
| 200 | {
|
|---|
| 201 | if (ownobj && dobj) delete dobj;
|
|---|
| 202 | }
|
|---|
| 203 |
|
|---|
| 204 | template <class T>
|
|---|
| 205 | AnyDataObj* FIO_TMatrix<T>::DataObj()
|
|---|
| 206 | {
|
|---|
| 207 | return(dobj);
|
|---|
| 208 | }
|
|---|
| 209 |
|
|---|
| 210 |
|
|---|
| 211 | template <class T>
|
|---|
| 212 | void FIO_TMatrix<T>::ReadSelf(PInPersist& is)
|
|---|
| 213 | {
|
|---|
| 214 | // On lit les 3 premiers uint_8
|
|---|
| 215 | // 0: Numero de version, 1 : NRows, 2 : NCol
|
|---|
| 216 | uint_4 itab[3];
|
|---|
| 217 | is.Get(itab,3);
|
|---|
| 218 | if (dobj == NULL) dobj = new TMatrix<T>(itab[1],itab[2]);
|
|---|
| 219 | else dobj->ReSize(itab[1],itab[2]);
|
|---|
| 220 | // On lit le NDataBlock
|
|---|
| 221 | FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
|
|---|
| 222 | fio_nd.Read(is);
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 | template <class T>
|
|---|
| 226 | void FIO_TMatrix<T>::WriteSelf(POutPersist& os) const
|
|---|
| 227 | {
|
|---|
| 228 | if (dobj == NULL) return;
|
|---|
| 229 | // On ecrit 3 uint_4 ....
|
|---|
| 230 | // 0: Numero de version, 1 : NRows, 2 : NCol
|
|---|
| 231 | uint_4 itab[3];
|
|---|
| 232 | itab[0] = 1; // Numero de version a 1
|
|---|
| 233 | itab[1] = dobj->NRows();
|
|---|
| 234 | itab[2] = dobj->NCols();
|
|---|
| 235 | os.Put(itab,3);
|
|---|
| 236 | // On ecrit le NDataBlock
|
|---|
| 237 | //cmv encore des problemes avec les consteries
|
|---|
| 238 | FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
|
|---|
| 239 | fio_nd.Write(os);
|
|---|
| 240 | }
|
|---|
| 241 |
|
|---|
| 242 | ///////////////////////////////////////////////////////////////
|
|---|
| 243 | #ifdef __CXX_PRAGMA_TEMPLATES__
|
|---|
| 244 | #pragma define_template TMatrix<uint_1>
|
|---|
| 245 | #pragma define_template TMatrix<uint_2>
|
|---|
| 246 | #pragma define_template TMatrix<int_2>
|
|---|
| 247 | #pragma define_template TMatrix<int_4>
|
|---|
| 248 | #pragma define_template TMatrix<int_8>
|
|---|
| 249 | #pragma define_template TMatrix<uint_4>
|
|---|
| 250 | #pragma define_template TMatrix<uint_8>
|
|---|
| 251 | #pragma define_template TMatrix<r_4>
|
|---|
| 252 | #pragma define_template TMatrix<r_8>
|
|---|
| 253 | #pragma define_template TMatrix< complex<float> >
|
|---|
| 254 | #pragma define_template TMatrix< complex<double> >
|
|---|
| 255 | // Instances des delegues FileIO (PPersist)
|
|---|
| 256 | #pragma define_template FIO_TMatrix<uint_1>
|
|---|
| 257 | #pragma define_template FIO_TMatrix<uint_2>
|
|---|
| 258 | #pragma define_template FIO_TMatrix<int_2>
|
|---|
| 259 | #pragma define_template FIO_TMatrix<int_4>
|
|---|
| 260 | #pragma define_template FIO_TMatrix<int_8>
|
|---|
| 261 | #pragma define_template FIO_TMatrix<uint_4>
|
|---|
| 262 | #pragma define_template FIO_TMatrix<uint_8>
|
|---|
| 263 | #pragma define_template FIO_TMatrix<r_8>
|
|---|
| 264 | #pragma define_template FIO_TMatrix<r_4>
|
|---|
| 265 | #pragma define_template FIO_TMatrix< complex<float> >
|
|---|
| 266 | #pragma define_template FIO_TMatrix< complex<double> >
|
|---|
| 267 | #endif
|
|---|
| 268 |
|
|---|
| 269 | #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
|
|---|
| 270 | template class TMatrix<uint_1>;
|
|---|
| 271 | template class TMatrix<uint_2>;
|
|---|
| 272 | template class TMatrix<int_2>;
|
|---|
| 273 | template class TMatrix<int_4>;
|
|---|
| 274 | template class TMatrix<int_8>;
|
|---|
| 275 | template class TMatrix<uint_4>;
|
|---|
| 276 | template class TMatrix<uint_8>;
|
|---|
| 277 | template class TMatrix<r_4>;
|
|---|
| 278 | template class TMatrix<r_8>;
|
|---|
| 279 | template class TMatrix< complex<float> >;
|
|---|
| 280 | template class TMatrix< complex<double> >;
|
|---|
| 281 | // Instances des delegues FileIO (PPersist)
|
|---|
| 282 | template class FIO_TMatrix<uint_1>;
|
|---|
| 283 | template class FIO_TMatrix<uint_2>;
|
|---|
| 284 | template class FIO_TMatrix<int_2>;
|
|---|
| 285 | template class FIO_TMatrix<int_4>;
|
|---|
| 286 | template class FIO_TMatrix<int_8>;
|
|---|
| 287 | template class FIO_TMatrix<uint_4>;
|
|---|
| 288 | template class FIO_TMatrix<uint_8>;
|
|---|
| 289 | template class FIO_TMatrix<r_8>;
|
|---|
| 290 | template class FIO_TMatrix<r_4>;
|
|---|
| 291 | template class FIO_TMatrix< complex<float> >;
|
|---|
| 292 | template class FIO_TMatrix< complex<double> >;
|
|---|
| 293 | #endif
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
