// $Id: tmatrix.cc,v 1.3 2000-04-03 17:35:58 ansari Exp $
//                         C.Magneville          04/99
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include "pexceptions.h"
#include "tmatrix.h"




////////////////////////////////////////////////////////////////
//**** Createur, Destructeur
template <class T>
TMatrix<T>::TMatrix()
// Constructeur par defaut.
  : TArray<T>()
{
}

template <class T>
TMatrix<T>::TMatrix(uint_4 r,uint_4 c, short mm)
// Construit une matrice de r lignes et c colonnes.
  :  TArray<T>() 
{
  if ( (r == 0) || (c == 0) )
    throw ParmError("TMatrix<T>::TMatrix(uint_4 r,uint_4 c) NRows or NCols = 0");
  ReSize(r, c, mm);
}

template <class T>
TMatrix<T>::TMatrix(const TMatrix<T>& a)
// Constructeur par copie (partage si "a" temporaire).
  : TArray<T>(a)
{
}

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)
{
}


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 ");
}

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 ");
  UpdateMemoryMapping(a, mm);
}

template <class T>
TMatrix<T>::~TMatrix()
// Destructeur
{

}

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 ");
  return(TArray<T>::Set(a));
}

template <class T>
void TMatrix<T>::ReSize(uint_4 r, uint_4 c, short mm)
{
  if(r==0||c==0) 
    throw(SzMismatchError("TMatrix::ReSize r or c==0 "));
  uint_4 size[BASEARRAY_MAXNDIMS];
  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
  size[0] = r;  size[1] = c;
  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
  UpdateMemoryMapping(mm, size[0], size[1]);
  TArray<T>::ReSize(2, size, 1);
  UpdateMemoryMapping(mm, size[0], size[1]);
}

template <class T>
void TMatrix<T>::Realloc(uint_4 r,uint_4 c, short mm, bool force)
{
  if(r==0||c==0) 
    throw(SzMismatchError("TMatrix::Realloc r or c==0 "));
  uint_4 size[BASEARRAY_MAXNDIMS];
  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
  UpdateMemoryMapping(mm, size[0], size[1]);
  TArray<T>::Realloc(2, size, 1, force);
  UpdateMemoryMapping(mm, size[0], size[1]);
}

// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
template <class T>
TMatrix<T> TMatrix<T>::operator () (Range rline, Range rcol) const
{
  uint_4 nx, ny;
  nx = ny = 0;
  if (GetMemoryMapping() == CMemoryMapping ) {
    TMatrix sm(SubArray(rcol, rline, 0, 0, 0),true,CMemoryMapping);
    sm.UpdateMemoryMapping(CMemoryMapping, nx, ny);
    sm.SetTemp(true);
    return(sm);
  }
  else {
    TMatrix sm(SubArray(rline, rcol, 0, 0, 0),true,CMemoryMapping);
    sm.UpdateMemoryMapping(FortranMemoryMapping, nx, ny);
    sm.SetTemp(true);
    return(sm);
  }
}

////////////////////////////////////////////////////////////////
// Transposition
template <class T>
TMatrix<T>& TMatrix<T>::Transpose() 
{
  uint_4 rci = macoli_;
  macoli_ = marowi_;
  marowi_ = rci;
  return(*this);
}


template <class T>
TMatrix<T> TMatrix<T>::Transpose(short mm)
{
  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
  TMatrix<T> tm(NCols(), NRows(), mm);
  for(uint_4 i=0; i<NRows(); i++)
    for(uint_4 j=0; j<NCols(); j++) 
      tm(j,i) = (*this)(i,j);
  tm.SetTemp(true);
  return tm;
}

// Rearrangement memoire
template <class T>
TMatrix<T> TMatrix<T>::Rearrange(short mm)
{
  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
  TMatrix<T> tm(NRows(), NCols(), mm);
  for(uint_4 i=0; i<NRows(); i++)
    for(uint_4 j=0; j<NCols(); j++) 
      tm(i,j) = (*this)(i,j);
  tm.SetTemp(true);
  return tm;
}

template <class T>
TMatrix<T>& TMatrix<T>::SetIdentity(IdentityMatrix imx)
{
  if (ndim_ == 0) {
    uint_4 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()") ;
  for(uint_4 i=0; i<NRows(); i++) (*this)(i,i) = diag;

  return (*this);
}



////////////////////////////////////////////////////////////////
//**** Impression
template <class T>
void TMatrix<T>::Print(ostream& os, int_4 maxprt, bool si) const
{
  os << "... TMatrix<T>(NRows=" << NRows() << " x NCols=" << NCols() 
     << ")::Print() (MemOrg=" << GetMemoryMapping() << ")" << endl;
  os << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() << " VectKA=" << VectKA() << endl;
  if (maxprt < 0)  maxprt = max_nprt_;
  int_4 npr = 0;
  Show(os, si);
  uint_4 kc,kr;  
  for(kr=0; kr<size_[marowi_]; kr++) {
    if ( (size_[marowi_] > 1) && (size_[macoli_] > 10) ) cout << "----- Ligne Line= " << kr << endl;
    for(kc=0; kc<size_[macoli_]; kc++) {
      if(kc > 0) os << ", ";  
      os << (*this)(kr, kc);   npr++; 
      if (npr >= maxprt) {
	if (npr < totsize_)  os << "\n     .... " << endl; return;
      }
    }
    os << endl;
  }
  os << "\n" << endl;
}

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

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;
  uint_4 r,c,k;
  uint_4 stepa = Step(ColsKA());
  uint_4 stepb = b.Step(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;
    }

  rm.SetTemp(true);
  return rm;
}


///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template TMatrix<uint_2>
#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<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