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