Changeset 286 in Sophya for trunk/SophyaLib/NTools/tmatrix.cc

Timestamp:

Apr 30, 1999, 1:02:52 PM (26 years ago)

Author:

ansari

Message:

tmatrix suite et fin cmv 30/4/99

File:

• trunk/SophyaLib/NTools/tmatrix.cc (modified) (23 diffs)

trunk/SophyaLib/NTools/tmatrix.cc

@@ -1 @@
-// $Id: tmatrix.cc,v 1.1 1999-04-29 08:57:15 ansari Exp $
+// $Id: tmatrix.cc,v 1.2 1999-04-30 11:02:52 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -18 +18 @@
 TMatrix<T>::TMatrix()
 // Constructeur par defaut.
-: mNr(0), mNc(0), mNDBlock(NULL)
+: mNr(0), mNc(0), mNDBlock()
 {
 }
@@ -25 +25 @@
 TMatrix<T>::TMatrix(uint_4 r,uint_4 c)
 // Construit une matrice de r lignes et c colonnes.
-: mNr(r), mNc(c), mNDBlock(new NDataBlock<T>(r*c))
+: mNr(r), mNc(c), mNDBlock(r*c)
 {
 }
@@ -33 +33 @@
 // Construit une matrice de r lignes et c colonnes. On fournit
 // le tableau des valeurs et eventuellement un Bridge.
-: mNr(r), mNc(c), mNDBlock(new NDataBlock<T>(r*c,values,br))
+: mNr(r), mNc(c), mNDBlock(r*c,values,br)
 {
 }
@@ -39 +39 @@
 template <class T>
 TMatrix<T>::TMatrix(const TMatrix<T>& a)
-// Constructeur par copie.
-: mNr(a.mNr), mNc(a.mNc)
-, mNDBlock(new NDataBlock<T>(*(a.DataBlock())))
+// Constructeur par copie (partage des donnees).
+: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
 {
 }
@@ -48 +47 @@
 TMatrix<T>::TMatrix(const TMatrix<T>& a,bool share)
 // Constructeur par copie.
-: mNr(a.mNr), mNc(a.mNc)
-, mNDBlock(new NDataBlock<T>(*(a.DataBlock()),share))
+: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
 {
 }
@@ -57 +55 @@
 // Destructeur
 {
-Delete();
-}
-
-template <class T>
-void TMatrix<T>::Delete(void)
-// Pour desallouer
-{
-mNr = mNc = 0;
-if(mNDBlock) delete mNDBlock; mNDBlock = NULL;
 }
 
@@ -71 +60 @@
 template <class T>
 void TMatrix<T>::Clone(const TMatrix<T>& a)
-// Clone (copie de donnee) a partir de "a"
-{
-if(a.mNr==0 || a.mNc==0 || a.DataBlock()==NULL) {Delete(); return;}
-if(mNDBlock) Delete();
-else NDataBlock<T>* mNDBlock = new NDataBlock<T>(1); // cas "a" vide
-mNr = a.mNr; mNc = a.mNc;
-mNDBlock->Clone(*(a.DataBlock()));
+// Clone (copie de donnees) a partir de "a"
+// Partage des donnees si "a" est temporaire.
+{
+mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;
 }
 
@@ -84 +70 @@
 // Reset de la matrice a "v"
 {
-if(mNDBlock==NULL)
-  throw(NullPtrError("TMatrix::Reset mNDBlock==NULL\n"));
-if(mNr==0 || mNc==0)
-  throw(SzMismatchError("TMatrix::Reset mNr==0 || mNc==0\n"));
-mNDBlock->Reset(v);
+mNDBlock.Reset(v);
 }
 
@@ -96 +78 @@
 {
 if(r==0 || c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
-if(mNDBlock==NULL) NDataBlock<T>* mNDBlock = new NDataBlock<T>(1);
 if(!force_alloc && mNr==r && mNc==c)  return;
-mNr = r; mNc = c;
-mNDBlock->ReSize(r*c);
+mNr = r; mNc = c; mNDBlock.ReSize(r*c,force_alloc);
 }
 
@@ -109 +89 @@
 // Operateur d'affectation depuis scalaire : identite * scalaire.
 {
-if(mNr!=mNc || mNr==0 || mNDBlock==NULL)
-  throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0 ou mNDBlock==NULL\n"));
-for(int r=0;r<mNr;r++) for(int c=0;c<mNc;c++) (*this)(r,c) = (r==c)? x: 0;
+if(mNr!=mNc || mNr==0)
+  throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0\n"));
+for(uint_4 r=0;r<mNr;r++) for(uint_4 c=0;c<mNc;c++)
+                               (*this)(r,c) = (r==c)? x: 0;
 return *this;
 }
@@ -117 +98 @@
 template <class T>
 TMatrix<T>& TMatrix<T>::operator = (const TMatrix<T>& a)
-// Operateur d'affectation A=B (sans reaffectation)
+// Operateur d'affectation A=a
+// (Re)allocation de la place memoire
+// sauf si "a" est cree temporairement (partage)
 {
 if(this == &a) return *this;
-if(mNr!=a.mNc || mNc!=a.mNc || mNr==0 || mNc==0)
-  throw(SzMismatchError("TMatrix::operator= mNc!=a.mNc ou mNr!=a.mNr ou ==0\n"));
-if(mNDBlock==NULL)
-  throw(SzMismatchError("TMatrix::operator= mNDBlock==NULL\n"));
-*(mNDBlock) = *(a.DataBlock());
+Clone(a);
 return *this;
 }
@@ -138 +117 @@
   os<<"TMatrix::Print("<<mNr<<","<<mNc<<")"<<endl;
 if(lp>0)
-  {os<<"       this="<<this<<" mNDBlock="<<mNDBlock<<endl;
-   if(mNDBlock) mNDBlock->Print(0,0);}
+  {os<<"       this="<<this<<endl; mNDBlock.Print(0,0);}
 if(mNr==0 || mNc==0) return;
 if(i0>=mNr || j0>=mNc || ni==0 || nj==0) return;
 uint_4 i1 = i0+ni; if(i1>mNr) i1 = mNr;
 uint_4 j1 = j0+nj; if(j1>mNc) j1 = mNc;
-for(int i=i0;i<i1;i++) {
-  for(int j=j0;j<j1;j++) cout<<" "<<(*this)(i,j)<<endl;
+for(uint_4 i=i0;i<i1;i++) {
+  for(uint_4 j=j0;j<j1;j++) cout<<" "<<(*this)(i,j);
+  cout<<endl;
 }
 }
@@ -155 +134 @@
 TMatrix<T>& TMatrix<T>::operator += (T b)
 {
-if(mNr==0 || mNc==0 || !mNDBlock)
-  throw(SzMismatchError("TMatrix::operator+=v mNr==0 || mNc==0 || !mNDBlock\n"));
-*mNDBlock += b;
+mNDBlock += b;
 return *this;
 }
@@ -164 +141 @@
 TMatrix<T>& TMatrix<T>::operator -= (T b)
 {
-if(mNr==0 || mNc==0 || !mNDBlock)
-  throw(SzMismatchError("TMatrix::operator-=v mNr==0 || mNc==0 || !mNDBlock\n"));
-*mNDBlock -= b;
+mNDBlock -= b;
 return *this;
 }
@@ -173 +148 @@
 TMatrix<T>& TMatrix<T>::operator *= (T b)
 {
-if(mNr==0 || mNc==0 || !mNDBlock)
-  throw(SzMismatchError("TMatrix::operator*=v mNr==0 || mNc==0 || !mNDBlock\n"));
-*mNDBlock *= b;
+mNDBlock *= b;
 return *this;
 }
@@ -182 +155 @@
 TMatrix<T>& TMatrix<T>::operator /= (T b)
 {
-if(b==(T) 0) throw(ParmError("TMatrix::operator/=v divide by zero\n"));
-if(mNr==0 || mNc==0 || !mNDBlock)
-  throw(SzMismatchError("TMatrix::operator/=v mNr==0 || mNc==0 || !mNDBlock\n"));
-*mNDBlock /= b;
+mNDBlock /= b;
 return *this;
 }
@@ -195 +165 @@
 TMatrix<T>& TMatrix<T>::operator += (const TMatrix<T>& a)
 {
-if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc || !mNDBlock || !a.mNDBlock)
-  throw(SzMismatchError("TMatrix::operator+=A size mismatch/null"));
-*(mNDBlock) += *(a.mNDBlock);
+if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
+  throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
+mNDBlock += a.mNDBlock;
 return *this;
 }
@@ -204 +174 @@
 TMatrix<T>& TMatrix<T>::operator -= (const TMatrix<T>& a)
 {
-if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc || !mNDBlock || !a.mNDBlock)
-  throw(SzMismatchError("TMatrix::operator-=A size mismatch/null"));
-*(mNDBlock) -= *(a.mNDBlock);
+if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
+  throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
+mNDBlock -= a.mNDBlock;
 return *this;
 }
@@ -215 +185 @@
 {
 uint_4 ndata = mNr*mNc;
-if(ndata==0 || mNc!=a.mNr || a.mNr!=a.mNc || !mNDBlock || !a.mNDBlock)
-  throw(SzMismatchError("TMatrix::operator*=A size mismatch/null"));
+if(ndata==0 || mNc!=a.mNr || a.mNr!=a.mNc)
+  throw(SzMismatchError("TMatrix::operator*=A size mismatch"));
 // A(i,j) = Sum(k) A(i,k)*B(k,j) ... il faut sauver la ligne "i" de A
 // Vecteur oi   : vecteur ou est sauve la ligne i de la matrice *this
@@ -238 +208 @@
 }
 
+////////////////////////////////////////////////////////////////
+//**** Surcharge de +,-,*: TMatrix = TMatrix + TMatrix;
+
+template <class T>
+TMatrix<T> TMatrix<T>::Add(const TMatrix<T>& b) const
+// C = A(this)+B
+{
+if(mNr!=b.mNr || mNc!=b.mNc)
+  throw(SzMismatchError("NDataBlock operator C=A+B size mismatch\n"));
+TMatrix<T> result; result.SetTemp(true);
+if(b.IsTemp()) {
+  result.Clone(b);
+  result += *this;
+} else {
+  result.Clone(*this);
+  result += b;
+}
+return result;
+}
+
+template <class T>
+TMatrix<T> TMatrix<T>::Sub(const TMatrix<T>& b) const
+// C = A(this)-B
+{
+if(mNr!=b.mNr || mNc!=b.mNc)
+  throw(SzMismatchError("NDataBlock operator C=A-B size mismatch\n"));
+TMatrix<T> result; result.SetTemp(true);
+if(b.IsTemp()) {
+  result.Clone(b);
+  result.mNDBlock = (*this).mNDBlock - result.mNDBlock;
+} else {
+  result.Clone(*this);
+  result -= b;
+}
+return result;
+}
+
+template <class T>
+TMatrix<T> TMatrix<T>::Mul(const TMatrix<T>& b) const
+// C = A(this)*B : Cij = Aik Bkj (allocation forcee dans tous les cas)
+{
+if(mNr==0 || mNc==0 || b.mNc==0 || mNc!=b.mNr)
+  throw(SzMismatchError("NDataBlock operator C=A*B size mismatch\n"));
+TMatrix<T> r; r.SetTemp(true); r.Realloc(mNr,b.mNc,true);
+T *ai,*aik,*bj,*bkj,*ri,*rij;
+for(ri=const_cast<T *>(r.Data()),ai=const_cast<T *>(Data());
+    ri<r.Data()+r.mNr*r.mNc;ri+=r.mNc,ai+=mNc) {
+  for(rij=ri,bj=const_cast<T *>(b.Data());rij<ri+r.mNc;rij++,bj++) {
+    *rij = 0;
+    for(aik=ai,bkj=bj;aik<ai+mNc;aik++,bkj+=b.mNc) *rij += *aik * *bkj;
+  }
+}
+return r;
+}
+
+////////////////////////////////////////////////////////////////
+// -------------------------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// -------------------------------------------------------------------------
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix()
+{
+dobj=new TMatrix<T>;
+ownobj=true;
+}
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix(string const & filename) 
+{
+dobj=new TMatrix<T>;
+ownobj=true; 
+Read(filename);
+}
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix(const TMatrix<T> & obj) 
+{ 
+dobj = new TMatrix<T>(obj);
+ownobj=true; 
+}
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix(TMatrix<T> * obj) 
+{ 
+dobj = obj;
+ownobj=false; 
+}
+
+template <class T>
+FIO_TMatrix<T>::~FIO_TMatrix()
+{
+if (ownobj && dobj) delete dobj;
+}
+
+template <class T>
+AnyDataObj* FIO_TMatrix<T>::DataObj()
+{
+return(dobj);
+}
+
+
+template <class T>
+void FIO_TMatrix<T>::ReadSelf(PInPersist& is)
+{
+// On lit les 3 premiers uint_8
+uint_4 itab[2];
+is.Get(itab,2);
+if (dobj == NULL) dobj = new TMatrix<T>(itab[0],itab[1]);
+else dobj->Realloc(itab[0],itab[1],false);
+// On lit le NDataBlock
+//cmv   encore des problemes avec les consteries
+//FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+//fio_nd.ReadSelf(is);
+}
+
+template <class T>
+void FIO_TMatrix<T>::WriteSelf(POutPersist& os) const
+{
+if (dobj == NULL)   return;
+//  On ecrit 2 uint_4 .... 0 : NRows,  1 : NCol
+uint_4 itab[2];
+itab[0] = dobj->NRows();
+itab[1] = dobj->NCol();
+os.Put(itab,2);
+// On ecrit le NDataBlock
+//cmv   encore des problemes avec les consteries
+//FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+//fio_nd.WriteSelf(os);
+}
+
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__
-/*
 #pragma define_template TMatrix<uint_1>
 #pragma define_template TMatrix<uint_2>
@@ -249 +349 @@
 #pragma define_template TMatrix<uint_8>
 #pragma define_template TMatrix<r_4>
-*/
 #pragma define_template TMatrix<r_8>
-/*
-#pragma define_template TMatrix< complex<float> >
-#pragma define_template TMatrix< complex<double> >
-*/
+#pragma define_template TMatrix< complex<float> > 
+#pragma define_template TMatrix< complex<double> > 
 // Instances des delegues FileIO (PPersist)
+#pragma define_template FIO_TMatrix<uint_1>
+#pragma define_template FIO_TMatrix<uint_2>
+#pragma define_template FIO_TMatrix<int_2>
+#pragma define_template FIO_TMatrix<int_4>
+#pragma define_template FIO_TMatrix<int_8>
+#pragma define_template FIO_TMatrix<uint_4>
+#pragma define_template FIO_TMatrix<uint_8>
+#pragma define_template FIO_TMatrix<r_8>
+#pragma define_template FIO_TMatrix<r_4>
+#pragma define_template FIO_TMatrix< complex<float> >
+#pragma define_template FIO_TMatrix< complex<double> >
 #endif
-
 
 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
@@ -272 +379 @@
 template class TMatrix< complex<double> >;
 // Instances des delegues FileIO (PPersist)
+template class FIO_TMatrix<uint_1>;
+template class FIO_TMatrix<uint_2>;
+template class FIO_TMatrix<int_2>;
+template class FIO_TMatrix<int_4>;
+template class FIO_TMatrix<int_8>;
+template class FIO_TMatrix<uint_4>;
+template class FIO_TMatrix<uint_8>;
+template class FIO_TMatrix<r_8>;
+template class FIO_TMatrix<r_4>;
+template class FIO_TMatrix< complex<float> >;
+template class FIO_TMatrix< complex<double> >;
 #endif