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

// $Id: tmatrix.cc,v 1.23 2002-11-15 09:35:44 ansari Exp $
//                         C.Magneville          04/99
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include "pexceptions.h"
#include "tmatrix.h"

/*!
  \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()
{
  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)
{
  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)
{
  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
 */
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;
}

///////////////////////////////////////////////////////////////
#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