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