source: Sophya/trunk/SophyaLib/TArray/tmatrix.cc@ 2574

// $Id: tmatrix.cc,v 1.25 2004-07-29 08:40:49 cmv Exp $
//                         C.Magneville          04/99
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include "pexceptions.h"
#include "tmatrix.h"

//#define DO_NOT_OPTIMIZE_PRODUCT

/*!
  \class SOPHYA::TMatrix
  \ingroup TArray
  
  The TMatrix class specializes the TArray class for representing
  two dimensional arrays as matrices. Matrix and vector operations,
  such as matrix multiplication or transposition is implemented.
  \b Matrix is a typedef for double precision floating point matrix ( TMatrix<r_8> ).

  \sa SOPHYA::TArray  SOPHYA::TVector
  \sa SOPHYA::Range  \sa SOPHYA::Sequence
  \sa SOPHYA::MathArray  \sa SOPHYA::SimpleMatrixOperation

  The following sample code illustrates vector-matrix multiplication 
  and matrix inversion, using simple gauss inversion.
  \code
  #include "array.h"
  // ....
  int n = 5;      // Size of matrix and vectors here
  Matrix a(n,n);  
  a = RandomSequence(RandomSequence::Gaussian, 0., 2.5);
  Vector x(n);
  x = RegularSequence(1.,3.);
  Vector b = a*x;
  cout << " ----- Vector x = \n " << x << endl;
  cout << " ----- Vector b = a*x = \n " << b << endl;
  SimpleMatrixOperation<r_8> smo;
  Matrix inva = smo.Inverse(a);
  cout << " ----- Matrix Inverse(a) = \n " << inva << endl;
  cout << " ----- Matrix a*Inverse(a) = \n " << inva*a << endl;
  cout << " ----- Matrix Inverse(a)*b (=Inv(a)*a*x) = \n " << inva*b << endl;
  cout << " ----- Matrix x-Inverse(a)*b = (=0 ?)\n " << x-inva*b << endl;  
  \endcode
  
 */

////////////////////////////////////////////////////////////////
//**** Createur, Destructeur
//! Default constructor
template <class T>
TMatrix<T>::TMatrix()
// Constructeur par defaut.
  : TArray<T>()
{
  arrtype_ = 1;   // Type = Matrix
}

//! constructor of a matrix with  r lines et c columns.
/*!
  \param r : number of rows
  \param c : number of columns
  \param mm : define the memory mapping type
  \sa ReSize
 */
template <class T>
TMatrix<T>::TMatrix(sa_size_t r,sa_size_t c, short mm)
// Construit une matrice de r lignes et c colonnes.
  :  TArray<T>() 
{
  if ( (r == 0) || (c == 0) )
    throw ParmError("TMatrix<T>::TMatrix(sa_size_t r,sa_size_t c) NRows or NCols = 0");
  arrtype_ = 1;   // Type = Matrix
  ReSize(r, c, mm);
}

//! Constructor by copy
/*!
  \warning datas are \b SHARED with \b a.
  \sa NDataBlock::NDataBlock(const NDataBlock<T>&)
*/
template <class T>
TMatrix<T>::TMatrix(const TMatrix<T>& a)
// Constructeur par copie
  : TArray<T>(a)
{
  arrtype_ = 1;   // Type = Matrix
  UpdateMemoryMapping(a, SameMemoryMapping);
}

//! Constructor by copy
/*!
  \param share : if true, share data. If false copy data
 */
template <class T>
TMatrix<T>::TMatrix(const TMatrix<T>& a, bool share)
// Constructeur par copie avec possibilite de forcer le partage ou non.
: TArray<T>(a, share)
{
  arrtype_ = 1;   // Type = Matrix
  UpdateMemoryMapping(a, SameMemoryMapping);
}

//! Constructor of a matrix from a TArray \b a
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 ");
  if (a.NbDimensions() == 1) {
    size_[1] = 1;
    step_[1] = size_[0]*step_[0];
    ndim_ = 2;
  }
  arrtype_ = 1;   // Type = Matrix
  UpdateMemoryMapping(a, SameMemoryMapping);
}

//! Constructor of a matrix from a TArray \b a
/*!
  \param a : TArray to be copied or shared
  \param share : if true, share data. If false copy data
  \param mm : define the memory mapping type
 */
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");
  if (a.NbDimensions() == 1) {
    size_[1] = 1;
    step_[1] = size_[0]*step_[0];
    ndim_ = 2;
  }
  arrtype_ = 1;   // Type = Matrix
  UpdateMemoryMapping(a, mm);
}

//! Constructor of a matrix from a TArray \b a , with a different data type
template <class T>
TMatrix<T>::TMatrix(const BaseArray& a)
: TArray<T>()
{
  arrtype_ = 1;   // Type = Matrix
  UpdateMemoryMapping(a, SameMemoryMapping);
  SetBA(a);
}



//! Destructor
template <class T>
TMatrix<T>::~TMatrix()
{
}

//! Set matrix equal to \b a and return *this
/*!
  \warning Datas are copied (cloned) from \b a.
  \sa NDataBlock::operator=(const NDataBlock<T>&)
*/
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");
  if ((arrtype_ == 2) && (a.NbDimensions() > 1) && (a.Size(0) > 1) && (a.Size(1) > 1) )
    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) Size(0,1)>1 for Vector");
  TArray<T>::Set(a);
  if (NbDimensions() == 1) {
    size_[1] = 1;
    step_[1] = size_[0]*step_[0];
    ndim_ = 2;
  }
  UpdateMemoryMapping(*this, SameMemoryMapping);
  return(*this);
}

template <class T>
TArray<T>& TMatrix<T>::SetBA(const BaseArray& a)
{
  if (a.NbDimensions() > 2) 
    throw SzMismatchError("TMatrix<T>::SetBA(const BaseArray& a) a.NbDimensions() > 2");
  if ((arrtype_ == 2) && (a.NbDimensions() > 1) && (a.Size(0) > 1) && (a.Size(1) > 1) )
    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) Size(0,1)>1 for Vector");
  TArray<T>::SetBA(a);
  if (NbDimensions() == 1) {
    size_[1] = 1;
    step_[1] = size_[0]*step_[0];
    ndim_ = 2;
  }
  UpdateMemoryMapping(*this, SameMemoryMapping);
  return(*this);
}



//! Resize the matrix
/*!
  \param r : number of rows
  \param c : number of columns
  \param mm : define the memory mapping type
         (SameMemoryMapping,CMemoryMapping
         ,FortranMemoryMapping,DefaultMemoryMapping)
 */
template <class T>
void TMatrix<T>::ReSize(sa_size_t r, sa_size_t c, short mm)
{
  if(r==0||c==0) 
    throw(SzMismatchError("TMatrix::ReSize r or c==0 "));
  if ((arrtype_ == 2) && (r > 1) && (c > 1))
    throw(SzMismatchError("TMatrix::ReSize r>1&&c>1 for Vector "));
  sa_size_t size[BASEARRAY_MAXNDIMS];
  for(int_4 kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) 
    mm = GetDefaultMemoryMapping();
  if (mm == CMemoryMapping) {
    size[0] = c;  size[1] = r;
  }
  else {
    size[0] = r;  size[1] = c;
  }
  TArray<T>::ReSize(2, size, 1);
  UpdateMemoryMapping(mm);
}

//! Re-allocate space for the matrix
/*!
  \param r : number of rows
  \param c : number of columns
  \param mm : define the memory mapping type
  \param force : if true re-allocation is forced, if not it occurs
          only if the required space is greater than the old one.
  \sa ReSize
 */
template <class T>
void TMatrix<T>::Realloc(sa_size_t r,sa_size_t c, short mm, bool force)
{
  if(r==0||c==0) 
    throw(SzMismatchError("TMatrix::Realloc r or c==0 "));
  if ((arrtype_ == 2) && (r > 1) && (c > 1))
    throw(SzMismatchError("TMatrix::Realloc r>1&&c>1 for Vector "));
  sa_size_t size[BASEARRAY_MAXNDIMS];
  for(int_4 kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) 
    mm = GetDefaultMemoryMapping();
  if (mm == CMemoryMapping) {
    size[0] = c;  size[1] = r;
  }
  else {
    size[0] = r;  size[1] = c;
  }
  TArray<T>::Realloc(2, size, 1, force);
  UpdateMemoryMapping(mm);
}

// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
//! Return a submatrix define by \b Range \b rline and \b rcol
template <class T>
TMatrix<T> TMatrix<T>::SubMatrix(Range rline, Range rcol) const
{
  short mm = GetMemoryMapping();
  Range rx, ry;
  if (mm == CMemoryMapping)  { rx = rcol;  ry = rline; }
  else { ry = rcol;  rx = rline; }
  TMatrix sm(SubArray(rx, ry, Range(0), Range(0), Range(0)),true, mm);
  sm.UpdateMemoryMapping(mm);
  return(sm);
}

////////////////////////////////////////////////////////////////
// Transposition
//! Transpose matrix in place, by changing the memory mapping
template <class T>
TMatrix<T>& TMatrix<T>::TransposeSelf() 
{
  short vt = (marowi_ == veceli_) ? ColumnVector : RowVector;
  int_4 rci = macoli_;
  macoli_ = marowi_;
  marowi_ = rci;
  veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
  return(*this);
}


//! Returns the transpose of the original matrix. 
/*!
  The data is shared between the two matrices
  \return return a new matrix
 */
template <class T>
TMatrix<T> TMatrix<T>::Transpose() const
{
  TMatrix<T> tm(*this);
  tm.TransposeSelf();
  return tm;
}

//! Returns a new matrix, corresponding to the transpose of the original matrix
/*!
  \param mm : define the memory mapping type
    (SameMemoryMapping,CMemoryMapping,FortranMemoryMapping)
  \return return a new matrix
 */
template <class T>
TMatrix<T> TMatrix<T>::Transpose(short mm) const
{
  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
  TMatrix<T> tm(NCols(), NRows(), mm);
  for(sa_size_t i=0; i<NRows(); i++)
    for(sa_size_t j=0; j<NCols(); j++) 
      tm(j,i) = (*this)(i,j);
  return tm;
}

//! Rearrange data in memory according to \b mm
/*!
  \param mm : define the memory mapping type
    (SameMemoryMapping,CMemoryMapping,FortranMemoryMapping)
  \warning If identical, return a matrix that share the datas
 */
template <class T>
TMatrix<T> TMatrix<T>::Rearrange(short mm) const 
{
  if ( mm == SameMemoryMapping)  mm = GetMemoryMapping();
  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) 
    mm = GetDefaultMemoryMapping();
  
  if  (mm == GetMemoryMapping())
    return (TMatrix<T>(*this, true));
  
  TMatrix<T> tm(NRows(), NCols(), mm);
  for(sa_size_t i=0; i<NRows(); i++)
    for(sa_size_t j=0; j<NCols(); j++) 
      tm(i,j) = (*this)(i,j);
  return tm;
}

//! Set the matrix to the identity matrix \b imx
template <class T>
TMatrix<T>& TMatrix<T>::SetIdentity(IdentityMatrix imx)
{
  if (ndim_ == 0) {
    sa_size_t 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()") ;
  *this = (T) 0;
  for(sa_size_t i=0; i<NRows(); i++) (*this)(i,i) = diag;

  return (*this);
}



////////////////////////////////////////////////////////////////
//**** Impression
//! Return info on number of rows, column and type \b T
template <class T>
string TMatrix<T>::InfoString() const
{
  string rs = "TMatrix<";
  rs += typeid(T).name();
  char buff[64];
  sprintf(buff, ">(NRows=%ld, NCols=%ld)", (long)NRows(), (long)NCols());
  rs += buff;
  return(rs);  
}

//! Print matrix
/*!
  \param os : output stream
  \param maxprt : maximum numer of print
  \param si : if true,  display attached DvList
  \param ascd : if true, suppresses the display of line numbers,
  suitable for ascii dump format.
  \sa SetMaxPrint
 */
template <class T>
void TMatrix<T>::Print(ostream& os, sa_size_t maxprt, bool si, bool ascd) const
{
  if (maxprt < 0)  maxprt = max_nprt_;
  sa_size_t npr = 0;
  Show(os, si);
  if (ndim_ < 1)  return;
  sa_size_t kc,kr;  
  for(kr=0; kr<size_[marowi_]; kr++) {
    if ( (size_[marowi_] > 1) && (size_[macoli_] > 10) && ascd) cout << "----- Line= " << kr << endl;
    for(kc=0; kc<size_[macoli_]; kc++) {
      if(kc > 0) os << " ";  
      os << (*this)(kr, kc);   npr++; 
      if (npr >= (sa_size_t) maxprt) {
        if (npr < totsize_)  os << "\n     .... " << endl; return;
      }
    }
    os << endl;
  }
  os << endl;
}

////////////////////////////////////////////////////////////////
//****  Multiplication matricielle  *****
////////////////////////////////////////////////////////////////

//! Return the matrix product C = (*this)*B
/*!
  \param mm : define the memory mapping type for the return matrix
 */
////////////// Routine de base sans optimisation //////////////
/*
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;
  sa_size_t r,c,k;
  sa_size_t stepa = Step(ColsKA());
  sa_size_t stepb = b.Step(b.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;
    }

  return rm;
}
*/

////////////// Routine optimisee //////////////
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() ") );

  // Commentaire: pas de difference de vitesse notable selon les mappings choisis pour la matrice produit "rm"
  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
  TMatrix<T> rm(NRows(), b.NCols(), mm);
  rm = (T) 0;  // Rien ne garantit que rm soit mise a zero a la creation !

  // Les "steps" pour l'adressage des colonnes de A et des lignes de B
  sa_size_t stepa = Step(ColsKA());
  sa_size_t stepb = b.Step(b.RowsKA());

  // On decide si on optimise ou non selon les dimensions de A et B
  // (il semble que optimiser ou non ne degrade pas notablement la vitesse pour les petites matrices)
#ifdef DO_NOT_OPTIMIZE_PRODUCT
  bool no_optim = true;
#else
  bool no_optim = false;
  //NON juste pour memoire: if((stepa+stepb)*NCols()*sizeof(T)<100000) no_optim=true;
#endif

  // Calcul de C=rm = A*B   (A=*this)
  // Remember: C-like      matrices are column packed
  //           Fortan-like matrices are line   packed
  sa_size_t r,c,k;
  T sum;
  const T * pe;

  // Pas d'optimisation demandee
  if( no_optim ) {
    //cout<<"no_optim("<<no_optim<<") "<<stepa<<" "<<stepb<<endl;
    const T * pea;
    const T * peb;
    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));
        peb = &(b(0,c));
        // On gagne un peu en remplacant "pea[k*stepa]" par "pea+=stepa" pour les grosses matrices
        //for(k=0; k<NCols(); k++)  sum += pea[k*stepa]*peb[k*stepb];
        for(k=0; k<NCols(); k++) {sum += (*pea)*(*peb); pea+=stepa; peb+=stepb;}
        rm(r,c) = sum;
      }
    }
  }
  // A.col est packed  et  B.row est packed (on a interet a optimiser quand meme)
  else if(stepa==1 && stepb==1) {
    //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
     T * pea = new T[rm.NCols()];
    const T * peb;
    for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
      pe = &((*this)(r,0));
      for(k=0; k<NCols(); k++) {pea[k] = *(pe++);}
      for(c=0; c<rm.NCols(); c++) {   // Boucle sur les colonnes de B
        sum = 0;
        peb = &(b(0,c));
        for(k=0; k<NCols(); k++) sum += *(peb++)*pea[k];
        rm(r,c) = sum;
      }
    }
    delete [] pea;
 }
  // A.col est packed  et  B.row n'est pas packed
  else if(stepa==1 && stepb!=1) {
    //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
    const T * pea;
    T * peb = new T[rm.NCols()];
    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
      pe = &(b(0,c));
      for(k=0; k<NCols(); k++) {peb[k] = *pe; pe+=stepb;}
      for(r=0; r<rm.NRows(); r++) {   // Boucle sur les lignes de A
        sum = 0;
        pea = &((*this)(r,0));
        for(k=0; k<NCols(); k++) sum += pea[k]*peb[k];
        rm(r,c) = sum;
      }
    }
    delete [] peb;
  }
  // A.col n'est pas packed  et  B.row est packed
  else if(stepa!=1 && stepb==1) {
    //cout<<"A.col not packed && B.row packed "<<stepa<<" "<<stepb<<endl;
    T * pea = new T[rm.NCols()];
    const T * peb;
    for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
      pe = &((*this)(r,0));
      for(k=0; k<NCols(); k++) {pea[k] = *pe; pe+=stepa;}
      for(c=0; c<rm.NCols(); c++) {   // Boucle sur les colonnes de B
        sum = 0;
        peb = &(b(0,c));
        for(k=0; k<NCols(); k++) sum += pea[k]*peb[k];
        rm(r,c) = sum;
      }
    }
    delete [] pea;
  }
  // A.col n'est pas packed  et  B.row n'est pas packed, stepb>stepa
  else if(stepa!=1 && stepb!=1 && stepb>stepa) {
    //cout<<"A.col not packed && B.row not packed ==> optimize on B "<<stepa<<" "<<stepb<<endl;
    const T * pea;
    T * peb = new T[rm.NCols()];
    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
      pe = &(b(0,c));
      for(k=0; k<NCols(); k++) {peb[k] = *pe; pe+=stepb;}
      for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
        sum = 0;
        pea = &((*this)(r,0));
        for(k=0; k<NCols(); k++) {sum += (*pea)*peb[k]; pea+=stepa;}
        rm(r,c) = sum;
      }
    }
    delete [] peb;
  }
  // A.col n'est pas packed  et  B.row n'est pas packed, stepa>=stepb
  else if(stepa!=1 && stepb!=1) {
    //cout<<"A.col not packed && B.row not packed ==> optimize on A "<<stepa<<" "<<stepb<<endl;
    T * pea = new T[rm.NCols()];
    const T * peb;
    for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
      pe = &((*this)(r,0));
      for(k=0; k<NCols(); k++) {pea[k] = *pe; pe+=stepa;}
      for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
        sum = 0;
        peb = &(b(0,c));
        for(k=0; k<NCols(); k++) {sum += pea[k]*(*peb); peb+=stepb;}
        rm(r,c) = sum;
      }
    }
    delete [] pea;
  }
  else {
    cout<<"TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! "<<endl;
    throw(SzMismatchError("TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! ") );
  }

  return rm;
}

///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template TMatrix<uint_2>
#pragma define_template TMatrix<uint_8>
#pragma define_template TMatrix<int_4>
#pragma define_template TMatrix<int_8>
#pragma define_template TMatrix<r_4>
#pragma define_template TMatrix<r_8> 
#pragma define_template TMatrix< complex<r_4> > 
#pragma define_template TMatrix< complex<r_8> > 
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class TMatrix<uint_2>;
template class TMatrix<uint_8>;
template class TMatrix<int_4>;
template class TMatrix<int_8>;
template class TMatrix<r_4>;
template class TMatrix<r_8>;
template class TMatrix< complex<r_4> >;
template class TMatrix< complex<r_8> >;
#endif