Changeset 804 in Sophya for trunk/SophyaLib/TArray/tmatrix.cc


Ignore:
Timestamp:
Apr 3, 2000, 7:36:01 PM (25 years ago)
Author:
ansari
Message:

Amelioation / debugging de la classe TArray<T> - TVector et TMatrix

heritent maintenant de TArray<T> - Classe RCMatrix rendu prive au fichier
sopemtx.cc - linfit.cc integre a sopemtx.cc

Reza 03/04/2000

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/tmatrix.cc

    r772 r804  
    1 // $Id: tmatrix.cc,v 1.2 2000-03-10 15:13:22 ansari Exp $
     1// $Id: tmatrix.cc,v 1.3 2000-04-03 17:35:58 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    44#include <stdio.h>
    55#include <stdlib.h>
    6 #include <iostream.h>
    7 #include <values.h>
    8 #include <complex>
    96#include "pexceptions.h"
    107#include "tmatrix.h"
    11 #include "objfio.h"
     8
     9
     10
    1211
    1312////////////////////////////////////////////////////////////////
    1413//**** Createur, Destructeur
    15 
    1614template <class T>
    1715TMatrix<T>::TMatrix()
    1816// 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)
     17  : TArray<T>()
     18{
     19}
     20
     21template <class T>
     22TMatrix<T>::TMatrix(uint_4 r,uint_4 c, short mm)
    2523// 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 {
     24  :  TArray<T>()
     25{
     26  if ( (r == 0) || (c == 0) )
     27    throw ParmError("TMatrix<T>::TMatrix(uint_4 r,uint_4 c) NRows or NCols = 0");
     28  ReSize(r, c, mm);
    3629}
    3730
    … …  
    3932TMatrix<T>::TMatrix(const TMatrix<T>& a)
    4033// 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)
     34  : TArray<T>(a)
     35{
     36}
     37
     38template <class T>
     39TMatrix<T>::TMatrix(const TMatrix<T>& a, bool share)
    4740// Constructeur par copie avec possibilite de forcer le partage ou non.
    48 : mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
    49 {
     41: TArray<T>(a, share)
     42{
     43}
     44
     45
     46template <class T>
     47TMatrix<T>::TMatrix(const TArray<T>& a)
     48: TArray<T>(a)
     49{
     50  if (a.NbDimensions() != 2)
     51    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions() != 2 ");
     52}
     53
     54template <class T>
     55TMatrix<T>::TMatrix(const TArray<T>& a, bool share, short mm )
     56: TArray<T>(a, share)
     57{
     58  if (a.NbDimensions() != 2)
     59    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions() != 2 ");
     60  UpdateMemoryMapping(a, mm);
    5061}
    5162
    … …  
    5465// Destructeur
    5566{
    56 }
    57 
    58 ////////////////////////////////////////////////////////////////
    59 // Operations matricielles
    60 template <class T>
    61 TMatrix<T> TMatrix<T>::Transpose(void) const
     67
     68}
     69
     70template <class T>
     71TArray<T>& TMatrix<T>::Set(const TArray<T>& a)
     72{
     73  if (a.NbDimensions() != 2)
     74    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) a.NbDimensions() != 2 ");
     75  return(TArray<T>::Set(a));
     76}
     77
     78template <class T>
     79void TMatrix<T>::ReSize(uint_4 r, uint_4 c, short mm)
     80{
     81  if(r==0||c==0)
     82    throw(SzMismatchError("TMatrix::ReSize r or c==0 "));
     83  uint_4 size[BASEARRAY_MAXNDIMS];
     84  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
     85  size[0] = r;  size[1] = c;
     86  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     87  UpdateMemoryMapping(mm, size[0], size[1]);
     88  TArray<T>::ReSize(2, size, 1);
     89  UpdateMemoryMapping(mm, size[0], size[1]);
     90}
     91
     92template <class T>
     93void TMatrix<T>::Realloc(uint_4 r,uint_4 c, short mm, bool force)
     94{
     95  if(r==0||c==0)
     96    throw(SzMismatchError("TMatrix::Realloc r or c==0 "));
     97  uint_4 size[BASEARRAY_MAXNDIMS];
     98  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
     99  UpdateMemoryMapping(mm, size[0], size[1]);
     100  TArray<T>::Realloc(2, size, 1, force);
     101  UpdateMemoryMapping(mm, size[0], size[1]);
     102}
     103
     104// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
     105template <class T>
     106TMatrix<T> TMatrix<T>::operator () (Range rline, Range rcol) const
     107{
     108  uint_4 nx, ny;
     109  nx = ny = 0;
     110  if (GetMemoryMapping() == CMemoryMapping ) {
     111    TMatrix sm(SubArray(rcol, rline, 0, 0, 0),true,CMemoryMapping);
     112    sm.UpdateMemoryMapping(CMemoryMapping, nx, ny);
     113    sm.SetTemp(true);
     114    return(sm);
     115  }
     116  else {
     117    TMatrix sm(SubArray(rline, rcol, 0, 0, 0),true,CMemoryMapping);
     118    sm.UpdateMemoryMapping(FortranMemoryMapping, nx, ny);
     119    sm.SetTemp(true);
     120    return(sm);
     121  }
     122}
     123
     124////////////////////////////////////////////////////////////////
    62125// 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 }
     126template <class T>
     127TMatrix<T>& TMatrix<T>::Transpose()
     128{
     129  uint_4 rci = macoli_;
     130  macoli_ = marowi_;
     131  marowi_ = rci;
     132  return(*this);
     133}
     134
     135
     136template <class T>
     137TMatrix<T> TMatrix<T>::Transpose(short mm)
     138{
     139  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     140  TMatrix<T> tm(NCols(), NRows(), mm);
     141  for(uint_4 i=0; i<NRows(); i++)
     142    for(uint_4 j=0; j<NCols(); j++)
     143      tm(j,i) = (*this)(i,j);
     144  tm.SetTemp(true);
     145  return tm;
     146}
     147
     148// Rearrangement memoire
     149template <class T>
     150TMatrix<T> TMatrix<T>::Rearrange(short mm)
     151{
     152  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     153  TMatrix<T> tm(NRows(), NCols(), mm);
     154  for(uint_4 i=0; i<NRows(); i++)
     155    for(uint_4 j=0; j<NCols(); j++)
     156      tm(i,j) = (*this)(i,j);
     157  tm.SetTemp(true);
     158  return tm;
     159}
     160
     161template <class T>
     162TMatrix<T>& TMatrix<T>::SetIdentity(IdentityMatrix imx)
     163{
     164  if (ndim_ == 0) {
     165    uint_4 sz = imx.Size();
     166    if (sz < 1) sz = 1;
     167    ReSize(sz, sz);
     168  }
     169  T diag = (T)imx.Diag();
     170  if (NRows() != NCols())
     171    throw SzMismatchError("TMatrix::operator= (IdentityMatrix) NRows() != NCols()") ;
     172  for(uint_4 i=0; i<NRows(); i++) (*this)(i,i) = diag;
     173
     174  return (*this);
     175}
     176
     177
    72178
    73179////////////////////////////////////////////////////////////////
    74180//**** 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"
     181template <class T>
     182void TMatrix<T>::Print(ostream& os, int_4 maxprt, bool si) const
     183{
     184  os << "... TMatrix<T>(NRows=" << NRows() << " x NCols=" << NCols()
     185     << ")::Print() (MemOrg=" << GetMemoryMapping() << ")" << endl;
     186  os << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() << " VectKA=" << VectKA() << endl;
     187  if (maxprt < 0)  maxprt = max_nprt_;
     188  int_4 npr = 0;
     189  Show(os, si);
     190  uint_4 kc,kr; 
     191  for(kr=0; kr<size_[marowi_]; kr++) {
     192    if ( (size_[marowi_] > 1) && (size_[macoli_] > 10) ) cout << "----- Ligne Line= " << kr << endl;
     193    for(kc=0; kc<size_[macoli_]; kc++) {
     194      if(kc > 0) os << ", "; 
     195      os << (*this)(kr, kc);   npr++;
     196      if (npr >= maxprt) {
     197        if (npr < totsize_)  os << "\n     .... " << endl; return;
     198      }
     199    }
     200    os << endl;
     201  }
     202  os << "\n" << endl;
     203}
     204
     205////////////////////////////////////////////////////////////////
     206//****  Multiplication matricielle  *****
     207////////////////////////////////////////////////////////////////
     208
     209template <class T>
     210TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
     211{
     212  if (NCols() != b.NRows())
     213    throw(SzMismatchError("TMatrix<T>::Multiply(b) NCols() != b.NRows() ") );
     214  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     215  TMatrix<T> rm(NRows(), b.NCols(), mm);
     216
     217  const T * pea;
     218  const T * peb;
     219  T sum;
     220  uint_4 r,c,k;
     221  uint_4 stepa = Step(ColsKA());
     222  uint_4 stepb = b.Step(RowsKA());
     223  // Calcul de C=rm = A*B   (A=*this)
     224  for(r=0; r<rm.NRows(); r++)      // Boucle sur les lignes de A
     225    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
     226      sum = 0;
     227      pea = &((*this)(r,0));       // 1er element de la ligne r de A
     228      peb = &(b(0,c));                // 1er element de la colonne c de B
     229      for(k=0; k<NCols(); k++)  sum += pea[k*stepa]*peb[k*stepb];
     230      rm(r,c) = sum;
     231    }
     232
     233  rm.SetTemp(true);
     234  return rm;
     235}
    256236
    257237
    258238///////////////////////////////////////////////////////////////
    259239#ifdef __CXX_PRAGMA_TEMPLATES__
    260 #pragma define_template TMatrix<uint_1>
    261240#pragma define_template TMatrix<uint_2>
    262 #pragma define_template TMatrix<int_2>
    263241#pragma define_template TMatrix<int_4>
    264242#pragma define_template TMatrix<int_8>
    265 #pragma define_template TMatrix<uint_4>
    266 #pragma define_template TMatrix<uint_8>
    267243#pragma define_template TMatrix<r_4>
    268 #pragma define_template TMatrix<r_8>
     244#pragma define_template TMatrix<r_8> 
    269245#pragma define_template TMatrix< complex<r_4> >
    270246#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> >
    283247#endif
    284248
    285249#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    286 template class TMatrix<uint_1>;
    287250template class TMatrix<uint_2>;
    288 template class TMatrix<int_2>;
    289251template class TMatrix<int_4>;
    290252template class TMatrix<int_8>;
    291 template class TMatrix<uint_4>;
    292 template class TMatrix<uint_8>;
    293253template class TMatrix<r_4>;
    294254template class TMatrix<r_8>;
    295255template class TMatrix< complex<r_4> >;
    296256template 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> >;
    309257#endif
Note: See TracChangeset for help on using the changeset viewer.