Index: /trunk/SophyaLib/NTools/tmatrix.cc
===================================================================
--- /trunk/SophyaLib/NTools/tmatrix.cc	(revision 279)
+++ /trunk/SophyaLib/NTools/tmatrix.cc	(revision 279)
@@ -0,0 +1,274 @@
+// $Id: tmatrix.cc,v 1.1 1999-04-29 08:57:15 ansari Exp $
+//                         C.Magneville          04/99
+#include "machdefs.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream.h>
+#include <complex>
+#include "pexceptions.h"
+#include "tmatrix.h"
+#include "objfio.h"
+
+using namespace PlanckDPC;
+
+////////////////////////////////////////////////////////////////
+//************ Createur, Destructeur, gestion des donnees
+
+template <class T>
+TMatrix<T>::TMatrix()
+// Constructeur par defaut.
+: mNr(0), mNc(0), mNDBlock(NULL)
+{
+}
+
+template <class T>
+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))
+{
+}
+
+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(new NDataBlock<T>(r*c,values,br))
+{
+}
+
+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())))
+{
+}
+
+template <class T>
+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))
+{
+}
+
+template <class T>
+TMatrix<T>::~TMatrix()
+// Destructeur
+{
+Delete();
+}
+
+template <class T>
+void TMatrix<T>::Delete(void)
+// Pour desallouer
+{
+mNr = mNc = 0;
+if(mNDBlock) delete mNDBlock; mNDBlock = NULL;
+}
+
+////////////////////////////////////////////////////////////////
+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()));
+}
+
+template <class T>
+void TMatrix<T>::Reset(T v)
+// 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);
+}
+
+template <class T>
+void TMatrix<T>::Realloc(uint_4 r,uint_4 c,bool force_alloc)
+// Reallocation de place
+{
+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);
+}
+
+////////////////////////////////////////////////////////////////
+//**** Surcharge de TMatrix=TMatrix; TMatrix=<T> b; 
+
+template <class T>
+TMatrix<T>& TMatrix<T>::operator = (T x)
+// 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;
+return *this;
+}
+
+template <class T>
+TMatrix<T>& TMatrix<T>::operator = (const TMatrix<T>& a)
+// Operateur d'affectation A=B (sans reaffectation)
+{
+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());
+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<<" mNDBlock="<<mNDBlock<<endl;
+   if(mNDBlock) 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;
+}
+}
+
+////////////////////////////////////////////////////////////////
+//**** Surcharge de +=,-=,*=,/= (INPLACE): TMatrix += <T> b;
+
+template <class T>
+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;
+return *this;
+}
+
+template <class T>
+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;
+return *this;
+}
+
+template <class T>
+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;
+return *this;
+}
+
+template <class T>
+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;
+return *this;
+}
+
+////////////////////////////////////////////////////////////////
+//**** Surcharge de +=,-=,*= (INPLACE): TMatrix += TMatrix;
+
+template <class T>
+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);
+return *this;
+}
+
+template <class T>
+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);
+return *this;
+}
+
+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 || !mNDBlock || !a.mNDBlock)
+  throw(SzMismatchError("TMatrix::operator*=A size mismatch/null"));
+// 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;
+}
+
+///////////////////////////////////////////////////////////////
+#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< complex<float> >
+#pragma define_template TMatrix< complex<double> >
+*/
+// Instances des delegues FileIO (PPersist)
+#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<float> >;
+template class TMatrix< complex<double> >;
+// Instances des delegues FileIO (PPersist)
+#endif
Index: /trunk/SophyaLib/NTools/tmatrix.h
===================================================================
--- /trunk/SophyaLib/NTools/tmatrix.h	(revision 279)
+++ /trunk/SophyaLib/NTools/tmatrix.h	(revision 279)
@@ -0,0 +1,75 @@
+// This may look like C code, but it is really -*- C++ -*-
+//                         C.Magneville          04/99
+#ifndef TMATRIX_SEEN
+#define TMATRIX_SEEN
+
+#include "machdefs.h"
+#include <stdio.h>
+#include <iostream.h>
+#include "ppersist.h"
+#include "anydataobj.h"
+#include "ndatablock.h"
+
+namespace PlanckDPC {
+
+template <class T>
+class TMatrix : public AnyDataObj {
+public:
+  // Creation / destruction 
+  TMatrix();
+  TMatrix(uint_4 r,uint_4 c);
+  TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br=NULL);
+  TMatrix(const TMatrix<T>& a);
+  TMatrix(const TMatrix<T>& a,bool share);
+  virtual ~TMatrix();
+
+  // Gestion taille/Remplissage
+  void Clone(const TMatrix<T>& a);
+  void Reset(T v=0);
+  void Realloc(uint_4 r,uint_4 c,bool force_alloc=true);
+
+  // Informations pointeur/data
+  inline int NRows() const {return mNr;}
+  inline int NCol()  const {return mNc;}
+  T const& operator()(uint_4 r,uint_4 c) const
+                     {return *(mNDBlock->Begin()+r*mNc+c);}
+  T&       operator()(uint_4 r,uint_4 c)
+                     {return *(mNDBlock->Begin()+r*mNc+c);}
+  inline       T* Data()       {return mNDBlock->Begin();}
+  inline const T* Data() const {return mNDBlock->Begin();}
+  inline       NDataBlock<T>* DataBlock()       {return mNDBlock;}
+  inline const NDataBlock<T>* DataBlock() const {return mNDBlock;}
+
+  // Operateur d'affectation
+  TMatrix<T>& operator = (T x);
+  TMatrix<T>& operator = (const TMatrix<T>& a);
+
+  // Impression
+  void Print(ostream& os,int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const;
+  inline void Print(int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const
+              {Print(cout,lp,i0,ni,j0,nj);}
+
+  // Surcharge d'operateurs
+  TMatrix<T>& operator += (T b);
+  TMatrix<T>& operator -= (T b);
+  TMatrix<T>& operator *= (T b);
+  TMatrix<T>& operator /= (T b);
+
+  TMatrix<T>& operator += (const TMatrix<T>& a);
+  TMatrix<T>& operator -= (const TMatrix<T>& a);
+  TMatrix<T>& operator *= (const TMatrix<T>& a);
+
+protected:
+  void Delete(void);
+
+  uint_4 mNr,mNc;
+  NDataBlock<T>* mNDBlock;
+};
+
+template<class T>
+inline ostream& operator << (ostream& os, const TMatrix<T>& a)
+                      {a.Print(os); return(os);}
+
+} // Fin du namespace
+
+#endif
