source: Sophya/trunk/SophyaLib/TArray/sopemtx.cc@ 785

//                         C.Magneville, R.Ansari        03/2000

#include "machdefs.h"
#include <stdio.h>
#include <iostream.h>
#include <complex>
#include <math.h>
#include "sopemtx.h"
#include "smathconst.h"

////////////////////////////////////////////////////////////////
// -------------------------------------------------------------
//   La classe de gestion des lignes et colonnes d'une matrice
// -------------------------------------------------------------
////////////////////////////////////////////////////////////////

template <class T> TMatrixRC<T>::TMatrixRC()
: matrix(NULL), data(NULL), index(0), step(0)
{}

template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
: matrix(&m), data(Org(m,rckind,ind)),
  index(ind), step(Step(m,rckind)), kind(rckind)
{
if (kind == TmatrixDiag && m.NCols() != m.NRows())
  throw(SzMismatchError("TMatrixRC::TMatrixRC(...,TmatrixDiag,...) size mismatch\n"));
}

////////////////////////////////////////////////////////////////
// Acces aux rangees et colonnes de matrices

template <class T> 
TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, uint_4 r) 
{
TMatrixRC<T> rc(m, TmatrixRow, r);
return rc;
}

template <class T> 
TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, uint_4 c) 
{
TMatrixRC<T> rc(m, TmatrixCol, c);
return rc;
}

template <class T> 
TMatrixRC<T> TMatrixRC<T>::Diag(TMatrix<T> & m) 
{
TMatrixRC<T> rc(m, TmatrixDiag);
return rc;
}


template <class T> int_4 TMatrixRC<T>::Next()
{
if (!matrix || kind==TmatrixDiag) return -1;
index++;
if(kind == TmatrixRow) {
  if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
  data += matrix->NCols();
} else {
  if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
  data++;
}
return index;
}

template <class T> int_4 TMatrixRC<T>::Prev()
{
if (!matrix || kind == TmatrixDiag) return -1;
index--;
if(index < 0) {index = 0; return -1;}
if(kind == TmatrixRow) data -= matrix->NCols();
else data--;
return index;
}

template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
{
if(!matrix) return -1;
if(c<0 || c>(int_4)matrix->NCols()) return -1;
kind = TmatrixCol;
index = c;
step = Step(*matrix, TmatrixCol);
data = Org(*matrix, TmatrixCol, c);
return c;
}

template <class T> int_4 TMatrixRC<T>::SetRow(int_4 r)
{
if(!matrix) return -1;
if(r<0 && r>(int_4)matrix->NRows()) return -1;
kind = TmatrixRow;
index = r;
step = Step(*matrix, TmatrixRow);
data = Org(*matrix, TmatrixRow, r);
return r;
}

template <class T> int_4 TMatrixRC<T>::SetDiag()
{
if (!matrix) return -1;
if (matrix->NCols() != matrix->NRows())
  throw(SzMismatchError("TMatrixRC::SetDiag size mismatch\n"));
kind = TmatrixDiag;
index = 0;
step = Step(*matrix, TmatrixDiag);
data = Org(*matrix, TmatrixDiag);
return 0;
}


template <class T> TMatrixRC<T>& TMatrixRC<T>::operator = (const TMatrixRC<T>& rc)
{
matrix = rc.matrix;
data   = rc.data;
index  = rc.index;
step   = rc.step;
kind   = rc.kind;
return *this;
}

template <class T> TVector<T> TMatrixRC<T>::GetVect() const
{
TVector<T> v(NElts());
for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
return v;
}

template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
{
if ( NElts() != rc.NElts() )
  throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
if ( kind != rc.kind )
  throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
return *this;
}

template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
{
if( NElts() != rc.NElts() )
  throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
if( kind != rc.kind )
  throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
return *this;
}


template <class T> TMatrixRC<T>& TMatrixRC<T>::operator *= (T x)
{
for(uint_4 i=0; i<NElts(); i++) (*this)(i) *= x;
return *this;
}

template <class T> TMatrixRC<T>& TMatrixRC<T>::operator /= (T x)
{
for(uint_4 i=0; i<NElts(); i++) (*this)(i) /= x;
return *this;
}


template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
{
for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
return *this;
}

template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
{
for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
return *this;
}

template <class T>
TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
{
if ( NElts() != rc.NElts() )
  throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
if ( kind != rc.kind )
  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
for(uint_4 i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;
return *this;
}

template <class T>
TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
{
if ( NElts() != rc.NElts() )
  throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
if ( kind != rc.kind )
  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
for(uint_4 i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;
return *this;
}

template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
{
if (first>NElts())
  throw(SzMismatchError("TMatrixRC::IMaxAbs size mismatch\n"));
uint_4 imax = first;
double vmax = Abs_Value((*this)(first));
for(uint_4 i=first+1; i<NElts(); i++) {
  double v = Abs_Value((*this)(i));
  if(v > vmax) {vmax = v; imax = i;}
}
return imax;
}

template <class T>
void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
{
if(rc1.NElts() != rc2.NElts())
  throw(SzMismatchError("TMatrixRC::Swap size mismatch\n"));
if(rc1.kind != rc2.kind)
  throw(SzMismatchError("TMatrixRC::Swap type mismatch\n"));
if(rc1.data == rc2.data) return;
for(uint_4 i=0; i<rc1.NElts(); i++)
  {T tmp = rc1(i); rc1(i) = rc2(i); rc2(i) = tmp;}
}




////////////////////////////////////////////////////////////////
//**** Pour inversion
#ifndef M_LN2
#define M_LN2 0.69314718055994530942
#endif
#ifndef LN_MINDOUBLE
#define LN_MINDOUBLE  (M_LN2 * (DMINEXP - 1))
#endif
#ifndef LN_MAXDOUBLE
#define LN_MAXDOUBLE  (M_LN2 * DMAXEXP)
#endif

template <class T>
T SimpleMatrixOperation<T>::GausPiv(TMatrix<T>& a, TMatrix<T>& b)
// Pivot de Gauss
// * Attention: egcs impose que cette fonction soit mise dans le .cc
//              avant ::Inverse() (car Inverse() l'utilise)
// {TMatrix A(a); TMatrix B(b); return (T) TMatrix::GausPiv(A,B);}
{
uint_4 n = a.NRows();
if(n!=b.NRows())
  throw(SzMismatchError("TMatrix::GausPiv size mismatch\n"));
// On fait une normalisation un peu brutale...
double vmin=MAXDOUBLE;
double vmax=0;
for(uint_4 iii=0; iii<a.NRows(); iii++)
  for(uint_4 jjj=0; jjj<a.NCols(); jjj++) {
    double v = TMatrixRC<T>::Abs_Value(a(iii,jjj));
    if(v>vmax) vmax = v;
    if(v<vmin && v>0) vmin = v;
}
double nrm = sqrt(vmin*vmax);
if(nrm > 1.e5 || nrm < 1.e-5) {
  a /= nrm;
  b /= nrm;
  //cout << "normalisation matrice " << nrm << endl;
} else nrm=1;

double det = 1.0;
if(nrm != 1) {
  double ld = a.NRows() * log(nrm);
  if (ld <= LN_MINDOUBLE || ld >= LN_MAXDOUBLE) {
   // cerr << "TMatrix warning, overflow for det" << endl;
  } else {
    det = exp(ld);
  }
}

TMatrixRC<T> pivRowa(a,TmatrixRow);
TMatrixRC<T> pivRowb(b,TmatrixRow);

for(uint_4 k=0; k<n-1; k++) {
  uint_4 iPiv = TMatrixRC<T>::Col(a, k).IMaxAbs(k);
  if(iPiv != k) {
    TMatrixRC<T> aIPiv(TMatrixRC<T>::Row(a,iPiv));
    TMatrixRC<T> aK(TMatrixRC<T>::Row(a, k));
    TMatrixRC<T>::Swap(aIPiv,aK);
    TMatrixRC<T> bIPiv(TMatrixRC<T>::Row(b, iPiv));
    TMatrixRC<T> bK(TMatrixRC<T>::Row(b, k));
    TMatrixRC<T>::Swap(bIPiv,bK);
  }
  double pivot = a(k,k);
  if (fabs(pivot) < 1.e-50) return 0.0;
  //det *= pivot;
  pivRowa.SetRow(k); // to avoid constructors
  pivRowb.SetRow(k);
  for (uint_4 i=k+1; i<n; i++) {
    double r = -a(i,k)/pivot;
    TMatrixRC<T>::Row(a, i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
    TMatrixRC<T>::Row(b, i).LinComb(r, pivRowb);
  }
}
det *= a(n-1, n-1);

// on remonte
for(uint_4 kk=n-1; kk>0; kk--) {
  double pivot = a(kk,kk);
  if (fabs(pivot) <= 1.e-50) return 0.0;
  pivRowa.SetRow(kk); // to avoid constructors
  pivRowb.SetRow(kk);
  for(uint_4 jj=0; jj<kk; jj++) {
    double r = -a(jj,kk)/pivot;
    TMatrixRC<T>::Row(a, jj).LinComb(r, pivRowa);
    TMatrixRC<T>::Row(b, jj).LinComb(r, pivRowb);
  }
}

for(uint_4 l=0; l<n; l++) {
  if (fabs(a(l,l)) <= 1.e-50) return 0.0;
  TMatrixRC<T>::Row(b, l) /= a(l,l);
}

return det;
}

template <class T>
TMatrix<T> SimpleMatrixOperation<T>::Inverse(TMatrix<T> const & A) 
// Inversion
{
TMatrix<T> a(A);
TMatrix<T> b(a.NCols(),a.NRows());  b = 1.;
if(fabs(GausPiv(a,b)) < 1.e-50)
  throw(MathExc("TMatrix Inverse() Singular OMatrix"));
return b;
}



///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
// Instances gestion lignes/colonnes
#pragma define_template TMatrixRC<uint_1>
#pragma define_template TMatrixRC<uint_2>
#pragma define_template TMatrixRC<int_2>
#pragma define_template TMatrixRC<int_4>
#pragma define_template TMatrixRC<int_8>
#pragma define_template TMatrixRC<uint_4>
#pragma define_template TMatrixRC<uint_8>
#pragma define_template TMatrixRC<r_4>
#pragma define_template TMatrixRC<r_8>
#pragma define_template TMatrixRC< complex<r_4> > 
#pragma define_template TMatrixRC< complex<r_8> >
// #pragma define_template SimpleMatrixOperation<r_4>
#pragma define_template SimpleMatrixOperation<r_8>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
// Instances gestion lignes/colonnes
template class TMatrixRC<uint_1>;
template class TMatrixRC<uint_2>;
template class TMatrixRC<int_2>;
template class TMatrixRC<int_4>;
template class TMatrixRC<int_8>;
template class TMatrixRC<uint_4>;
template class TMatrixRC<uint_8>;
template class TMatrixRC<r_4>;
template class TMatrixRC<r_8>;
template class TMatrixRC< complex<r_4> >;
template class TMatrixRC< complex<r_8> >;
template class SimpleMatrixOperation<r_4>;
template class SimpleMatrixOperation<r_8>;
#endif