Index: trunk/SophyaLib/TArray/tmatrix.cc
===================================================================
--- trunk/SophyaLib/TArray/tmatrix.cc	(revision 772)
+++ trunk/SophyaLib/TArray/tmatrix.cc	(revision 804)
@@ -1,37 +1,30 @@
-// $Id: tmatrix.cc,v 1.2 2000-03-10 15:13:22 ansari Exp $
+// $Id: tmatrix.cc,v 1.3 2000-04-03 17:35:58 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
 #include <stdio.h>
 #include <stdlib.h>
-#include <iostream.h>
-#include <values.h>
-#include <complex>
 #include "pexceptions.h"
 #include "tmatrix.h"
-#include "objfio.h"
+
+
+
 
 ////////////////////////////////////////////////////////////////
 //**** Createur, Destructeur
-
 template <class T>
 TMatrix<T>::TMatrix()
 // Constructeur par defaut.
-: mNr(0), mNc(0), mNDBlock()
-{
-}
-
-template <class T>
-TMatrix<T>::TMatrix(uint_4 r,uint_4 c)
+  : TArray<T>()
+{
+}
+
+template <class T>
+TMatrix<T>::TMatrix(uint_4 r,uint_4 c, short mm)
 // Construit une matrice de r lignes et c colonnes.
-: mNr(r), mNc(c), mNDBlock(r*c)
-{
-}
-
-template <class T>
-TMatrix<T>::TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br)
-// Construit une matrice de r lignes et c colonnes. On fournit
-// le tableau des valeurs et eventuellement un Bridge.
-: mNr(r), mNc(c), mNDBlock(r*c,values,br)
-{
+  :  TArray<T>() 
+{
+  if ( (r == 0) || (c == 0) )
+    throw ParmError("TMatrix<T>::TMatrix(uint_4 r,uint_4 c) NRows or NCols = 0");
+  ReSize(r, c, mm);
 }
 
@@ -39,13 +32,31 @@
 TMatrix<T>::TMatrix(const TMatrix<T>& a)
 // Constructeur par copie (partage si "a" temporaire).
-: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
-{
-}
-
-template <class T>
-TMatrix<T>::TMatrix(const TMatrix<T>& a,bool share)
+  : TArray<T>(a)
+{
+}
+
+template <class T>
+TMatrix<T>::TMatrix(const TMatrix<T>& a, bool share)
 // Constructeur par copie avec possibilite de forcer le partage ou non.
-: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
-{
+: TArray<T>(a, share)
+{
+}
+
+
+template <class T>
+TMatrix<T>::TMatrix(const TArray<T>& a)
+: TArray<T>(a)
+{
+  if (a.NbDimensions() != 2) 
+    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions() != 2 ");
+}
+
+template <class T>
+TMatrix<T>::TMatrix(const TArray<T>& a, bool share, short mm )
+: TArray<T>(a, share)
+{
+  if (a.NbDimensions() != 2) 
+    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions() != 2 ");
+  UpdateMemoryMapping(a, mm);
 }
 
@@ -54,256 +65,193 @@
 // Destructeur
 {
-}
-
-////////////////////////////////////////////////////////////////
-// Operations matricielles
-template <class T>
-TMatrix<T> TMatrix<T>::Transpose(void) const
+
+}
+
+template <class T>
+TArray<T>& TMatrix<T>::Set(const TArray<T>& a)
+{
+  if (a.NbDimensions() != 2) 
+    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) a.NbDimensions() != 2 ");
+  return(TArray<T>::Set(a));
+}
+
+template <class T>
+void TMatrix<T>::ReSize(uint_4 r, uint_4 c, short mm)
+{
+  if(r==0||c==0) 
+    throw(SzMismatchError("TMatrix::ReSize r or c==0 "));
+  uint_4 size[BASEARRAY_MAXNDIMS];
+  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
+  size[0] = r;  size[1] = c;
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  UpdateMemoryMapping(mm, size[0], size[1]);
+  TArray<T>::ReSize(2, size, 1);
+  UpdateMemoryMapping(mm, size[0], size[1]);
+}
+
+template <class T>
+void TMatrix<T>::Realloc(uint_4 r,uint_4 c, short mm, bool force)
+{
+  if(r==0||c==0) 
+    throw(SzMismatchError("TMatrix::Realloc r or c==0 "));
+  uint_4 size[BASEARRAY_MAXNDIMS];
+  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
+  UpdateMemoryMapping(mm, size[0], size[1]);
+  TArray<T>::Realloc(2, size, 1, force);
+  UpdateMemoryMapping(mm, size[0], size[1]);
+}
+
+// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
+template <class T>
+TMatrix<T> TMatrix<T>::operator () (Range rline, Range rcol) const
+{
+  uint_4 nx, ny;
+  nx = ny = 0;
+  if (GetMemoryMapping() == CMemoryMapping ) {
+    TMatrix sm(SubArray(rcol, rline, 0, 0, 0),true,CMemoryMapping);
+    sm.UpdateMemoryMapping(CMemoryMapping, nx, ny);
+    sm.SetTemp(true);
+    return(sm);
+  }
+  else {
+    TMatrix sm(SubArray(rline, rcol, 0, 0, 0),true,CMemoryMapping);
+    sm.UpdateMemoryMapping(FortranMemoryMapping, nx, ny);
+    sm.SetTemp(true);
+    return(sm);
+  }
+}
+
+////////////////////////////////////////////////////////////////
 // Transposition
-{
-TMatrix<T> a; a.Clone(*this); a.SetTemp(true);
-a.mNr = mNc; a.mNc = mNr;
-{for(uint_4 i=0; i<mNr; i++)
-  for(uint_4 j=0; j<mNc; j++) {
-    a(j,i) = (*this)(i,j);
-}}
-return a;
-}
+template <class T>
+TMatrix<T>& TMatrix<T>::Transpose() 
+{
+  uint_4 rci = macoli_;
+  macoli_ = marowi_;
+  marowi_ = rci;
+  return(*this);
+}
+
+
+template <class T>
+TMatrix<T> TMatrix<T>::Transpose(short mm)
+{
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  TMatrix<T> tm(NCols(), NRows(), mm);
+  for(uint_4 i=0; i<NRows(); i++)
+    for(uint_4 j=0; j<NCols(); j++) 
+      tm(j,i) = (*this)(i,j);
+  tm.SetTemp(true);
+  return tm;
+}
+
+// Rearrangement memoire
+template <class T>
+TMatrix<T> TMatrix<T>::Rearrange(short mm)
+{
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  TMatrix<T> tm(NRows(), NCols(), mm);
+  for(uint_4 i=0; i<NRows(); i++)
+    for(uint_4 j=0; j<NCols(); j++) 
+      tm(i,j) = (*this)(i,j);
+  tm.SetTemp(true);
+  return tm;
+}
+
+template <class T>
+TMatrix<T>& TMatrix<T>::SetIdentity(IdentityMatrix imx)
+{
+  if (ndim_ == 0) {
+    uint_4 sz = imx.Size();
+    if (sz < 1) sz = 1;
+    ReSize(sz, sz);
+  }
+  T diag = (T)imx.Diag();
+  if (NRows() != NCols()) 
+    throw SzMismatchError("TMatrix::operator= (IdentityMatrix) NRows() != NCols()") ;
+  for(uint_4 i=0; i<NRows(); i++) (*this)(i,i) = diag;
+
+  return (*this);
+}
+
+
 
 ////////////////////////////////////////////////////////////////
 //**** Impression
-
-template <class T>
-void TMatrix<T>::Print(ostream& os,int lp
-                 ,uint_4 i0,uint_4 ni,uint_4 j0,uint_4 nj) const
-// Impression de la sous-matrice (i:i+ni-1,i:j+nj-1)
-{
-os<<"TMatrix::Print("<<mNr<<","<<mNc<<")"<<endl;
-if(lp>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(uint_4 i=i0;i<i1;i++) {
-  for(uint_4 j=j0;j<j1;j++) cout<<" "<<(*this)(i,j);
-  cout<<endl;
-}
-}
-
-////////////////////////////////////////////////////////////////
-//**** Surcharge de *= (INPLACE): TMatrix *= TMatrix;
-
-template <class T>
-TMatrix<T>& TMatrix<T>::operator *= (const TMatrix<T>& a)
-// A = A*B  -> A(n,m) = A(n,m)*B(m,m)
-{
-uint_4 ndata = mNr*mNc;
-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
-//                oi,oe = pointeur de debut et de fin du vecteur temporaire
-//                oij   = pointeur parcourant le vecteur oi
-// Matrice *this: i     = pointeur du debut de la ligne i
-//                ij    = pointeur parcourant la ligne i
-// Matrice a    : aj    = pointeur de debut de la colonne j
-//                aji   = pointeur parcourant la colonne j
-T* oi = new T[mNc]; T* oe = oi+mNc;
-for(T *i=Data(); i<Data()+ndata; i+=mNc) {
-  {for(T *oij=oi,*ij=i; oij<oe;) *oij++ = *ij++;}
-  {for(T *ij=i,*aj=const_cast<T *>(a.Data()); aj<a.Data()+a.mNc; ij++,aj++) {
-    T sum = 0;
-    for(T *oij=oi,*aji=aj; oij<oe; oij++,aji+=a.mNc) sum += *oij * *aji;
-    *ij = sum;
-  }}
-}
-delete [] oi;
-return *this;
-}
-
-////////////////////////////////////////////////////////////////
-//**** Pour surcharge d'operateurs C = A (+,-,*,/) B
-
-template <class T> TMatrix<T> TMatrix<T>::Add(const TMatrix<T>& b) const
-{
-if(mNr!=b.mNr || mNc!=b.mNc)
-    throw(SzMismatchError("TMatrix operator C=A+B size mismatch\n"));
-TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
-result.mNDBlock = mNDBlock+b.mNDBlock;
-return result;
-}
-
-template <class T> TMatrix<T> TMatrix<T>::Sub(const TMatrix<T>& b) const
-{
-if(mNr!=b.mNr || mNc!=b.mNc)
-    throw(SzMismatchError("TMatrix operator C=A-B size mismatch\n"));
-TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
-result.mNDBlock = mNDBlock-b.mNDBlock;
-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("TMatrix operator C=A*B size mismatch\n"));
-TMatrix<T> r; r.SetTemp(true); r.ReSize(mNr,b.mNc);
-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>::SetDataObj(AnyDataObj & o)
-{
-TMatrix<T> * po = dynamic_cast< TMatrix<T> * >(&o);
-if (po == NULL) return;
-if (ownobj && dobj) delete dobj;
-dobj = po; ownobj = false;
-} 
-
-template <class T>
-void FIO_TMatrix<T>::ReadSelf(PInPersist& is)
-{
-// On lit les 3 premiers uint_4
-//  0: Numero de version,  1 : NRows,  2 : NCol
-uint_4 itab[3];
-is.Get(itab,3);
-if (dobj == NULL) dobj = new TMatrix<T>(itab[1],itab[2]);
-else dobj->ReSize(itab[1],itab[2]);
-// On lit le NDataBlock
-FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
-fio_nd.Read(is);
-}
-
-template <class T>
-void FIO_TMatrix<T>::WriteSelf(POutPersist& os) const
-{
-if (dobj == NULL)   return;
-//  On ecrit 3 uint_4 .... 
-//  0: Numero de version,  1 : NRows,  2 : NCol
-uint_4 itab[3];
- itab[0] = 1;  // Numero de version a 1
-itab[1] = dobj->NRows();
-itab[2] = dobj->NCols();
-os.Put(itab,3);
-// On ecrit le NDataBlock
-FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
-fio_nd.Write(os);
-}
-
-////////////////////////////////////////////////////////////////
-// -------------------------------------------------------------
-//   La classe de gestion des lignes et colonnes d'une matrice
-// -------------------------------------------------------------
-////////////////////////////////////////////////////////////////
-#include "tvector.h"
+template <class T>
+void TMatrix<T>::Print(ostream& os, int_4 maxprt, bool si) const
+{
+  os << "... TMatrix<T>(NRows=" << NRows() << " x NCols=" << NCols() 
+     << ")::Print() (MemOrg=" << GetMemoryMapping() << ")" << endl;
+  os << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() << " VectKA=" << VectKA() << endl;
+  if (maxprt < 0)  maxprt = max_nprt_;
+  int_4 npr = 0;
+  Show(os, si);
+  uint_4 kc,kr;  
+  for(kr=0; kr<size_[marowi_]; kr++) {
+    if ( (size_[marowi_] > 1) && (size_[macoli_] > 10) ) cout << "----- Ligne Line= " << kr << endl;
+    for(kc=0; kc<size_[macoli_]; kc++) {
+      if(kc > 0) os << ", ";  
+      os << (*this)(kr, kc);   npr++; 
+      if (npr >= maxprt) {
+	if (npr < totsize_)  os << "\n     .... " << endl; return;
+      }
+    }
+    os << endl;
+  }
+  os << "\n" << endl;
+}
+
+////////////////////////////////////////////////////////////////
+//****  Multiplication matricielle  *****
+////////////////////////////////////////////////////////////////
+
+template <class T>
+TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
+{
+  if (NCols() != b.NRows()) 
+    throw(SzMismatchError("TMatrix<T>::Multiply(b) NCols() != b.NRows() ") );
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  TMatrix<T> rm(NRows(), b.NCols(), mm);
+
+  const T * pea;
+  const T * peb;
+  T sum;
+  uint_4 r,c,k;
+  uint_4 stepa = Step(ColsKA());
+  uint_4 stepb = b.Step(RowsKA());
+  // Calcul de C=rm = A*B   (A=*this)
+  for(r=0; r<rm.NRows(); r++)      // Boucle sur les lignes de A
+    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
+      sum = 0;
+      pea = &((*this)(r,0));       // 1er element de la ligne r de A
+      peb = &(b(0,c));                // 1er element de la colonne c de B
+      for(k=0; k<NCols(); k++)  sum += pea[k*stepa]*peb[k*stepb];
+      rm(r,c) = sum;
+    }
+
+  rm.SetTemp(true);
+  return rm;
+}
 
 
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__
-#pragma define_template TMatrix<uint_1>
 #pragma define_template TMatrix<uint_2>
-#pragma define_template TMatrix<int_2>
 #pragma define_template TMatrix<int_4>
 #pragma define_template TMatrix<int_8>
-#pragma define_template TMatrix<uint_4>
-#pragma define_template TMatrix<uint_8>
 #pragma define_template TMatrix<r_4>
-#pragma define_template TMatrix<r_8>
+#pragma define_template TMatrix<r_8> 
 #pragma define_template TMatrix< complex<r_4> > 
 #pragma define_template TMatrix< complex<r_8> > 
-// 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<r_4> >
-#pragma define_template FIO_TMatrix< complex<r_8> >
 #endif
 
 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
-template class TMatrix<uint_1>;
 template class TMatrix<uint_2>;
-template class TMatrix<int_2>;
 template class TMatrix<int_4>;
 template class TMatrix<int_8>;
-template class TMatrix<uint_4>;
-template class TMatrix<uint_8>;
 template class TMatrix<r_4>;
 template class TMatrix<r_8>;
 template class TMatrix< complex<r_4> >;
 template class TMatrix< complex<r_8> >;
-// 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<r_4> >;
-template class FIO_TMatrix< complex<r_8> >;
 #endif
