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

// $Id: tmatrix.cc,v 1.22 2001-08-06 22:48:43 cmv 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
  Class of matrixes
  \sa TArray
 */

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