source: Sophya/trunk/SophyaLib/NTools/tmatrix.h@ 303

// 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"
class GeneralFit;

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();

  // Temporaire?
  inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
  inline void SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}

  // Gestion taille/Remplissage
  inline void Clone(const TMatrix<T>& a)  // Clone: copie des donnees de "a"
              {mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;}
  inline void Reset(T v=0) {mNDBlock.Reset(v);}
  inline void ReSize(uint_4 r,uint_4 c)  // Reallocation de place
       {if(r==0||c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
        mNr = r; mNc = c; mNDBlock.ReSize(r*c);}

  // Informations pointeur/data
  inline uint_4 NRows() const {return mNr;}
  inline uint_4 NCols() const {return mNc;}
  inline T const& operator()(uint_4 r,uint_4 c) const
                            {return *(mNDBlock.Begin()+r*mNc+c);}
  inline T&       operator()(uint_4 r,uint_4 c)
                            {return *(mNDBlock.Begin()+r*mNc+c);}
  inline T const& operator[](uint_4 ip) const
                            {return *(mNDBlock.Begin()+ip);}
  inline T&       operator[](uint_4 ip)
                            {return *(mNDBlock.Begin()+ip);}
  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;}

  // Operations matricielles
  TMatrix<T> Transpose(void) const;

  // Operateur d'affectation
  // A = x (matrice diagonale x*Identite)
  inline TMatrix<T>& operator = (T x)
    {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;}
  // A = B : partage les donnees si "a" est temporaire, clone sinon.
  inline TMatrix<T>& operator = (const TMatrix<T>& a)
                {if(this == &a) return *this; CloneOrShare(a); return *this;}

  // 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 INPLACE: A (+=,-=,*=,/=) (T) x
  inline TMatrix<T>& operator += (T b) {mNDBlock += b; return *this;}
  inline TMatrix<T>& operator -= (T b) {mNDBlock -= b; return *this;}
  inline TMatrix<T>& operator *= (T b) {mNDBlock *= b; return *this;}
  inline TMatrix<T>& operator /= (T b) {mNDBlock /= b; return *this;}

  // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) B
  inline TMatrix<T>& operator += (const TMatrix<T>& a)
              {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
                  throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
               mNDBlock += a.mNDBlock; return *this;}
  inline TMatrix<T>& operator -= (const TMatrix<T>& a)
              {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
                  throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
               mNDBlock -= a.mNDBlock; return *this;}
  TMatrix<T>& operator *= (const TMatrix<T>& a);

  // Pour surcharge d'operateurs C = A (+,-,*) B
  TMatrix<T> Add(const TMatrix<T>& b) const;
  TMatrix<T> Sub(const TMatrix<T>& b) const;
  TMatrix<T> Mul(const TMatrix<T>& b) const;

  // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
  // operations sur la matrice B
  TMatrix<T> Inverse() const;
  static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);

  // Residus et fonction fittees.
  TMatrix<T> FitResidus(GeneralFit& gfit
            ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);
  TMatrix<T> FitFunction(GeneralFit& gfit
            ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);

protected:
  // partage les donnees si "a" temporaire, clone sinon.
  inline void CloneOrShare(const TMatrix<T>& a)
              {mNDBlock.CloneOrShare(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
  // Share: partage les donnees de "a"
  inline void Share(const TMatrix<T>& a)
              {mNDBlock.Share(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}

  uint_4 mNr,mNc;
  NDataBlock<T> mNDBlock;
};

////////////////////////////////////////////////////////////////
// Impression

template <class T>
inline ostream& operator << (ostream& os, const TMatrix<T>& a)
                            {a.Print(os); return(os);}

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs A (+=,-=,*=,/=) (T) x

template <class T> inline TMatrix<T> operator + (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}

template <class T> inline TMatrix<T> operator + (T b,const TMatrix<T>& a)
    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}

template <class T> inline TMatrix<T> operator - (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result -= b; return result;}

template <class T> inline TMatrix<T> operator - (T b,const TMatrix<T>& a)
    {TMatrix<T> result(a); result.SetTemp(true);
     result.DataBlock() = b-result.DataBlock(); return result;}

template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}

template <class T> inline TMatrix<T> operator * (T b,const TMatrix<T>& a)
    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}

template <class T> inline TMatrix<T> operator / (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result /= b; return result;}

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs C = A (+,-,*,/) B

template <class T>
inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
                  {return a.Add(b);}

template <class T>
inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
                  {return a.Sub(b);}

template <class T>
inline TMatrix<T> operator * (const TMatrix<T>& a,const TMatrix<T>& b)
                  {return a.Mul(b);}


/////////////////////////////////////////////////////////////////////////
// Classe pour la gestion de persistance
template <class T>
class FIO_TMatrix : public  PPersist  {
public:
  FIO_TMatrix();
  FIO_TMatrix(string const & filename); 
  FIO_TMatrix(const TMatrix<T> & obj);
  FIO_TMatrix(TMatrix<T> * obj);
  virtual ~FIO_TMatrix();
  virtual AnyDataObj* DataObj();
  inline operator TMatrix<T>() { return(*dobj); }
protected :
  virtual void ReadSelf(PInPersist&);           
  virtual void WriteSelf(POutPersist&) const;  
  TMatrix<T> * dobj;
  bool ownobj;
};

} // Fin du namespace

#endif