Changeset 304 in Sophya for trunk/SophyaLib/NTools/tmatrix.cc

Timestamp:

May 18, 1999, 7:19:37 PM (26 years ago)

Author:

ansari

Message:

TMatrixRC<T> + GausPiv interne a TMatrix (sans Matrix)
Resolution pb des friend template (namesapce!)
createur Vector a partir de TVector<r_8> cmv 18/5/99

File:

• trunk/SophyaLib/NTools/tmatrix.cc (modified) (8 diffs)

trunk/SophyaLib/NTools/tmatrix.cc

@@ -1 @@
-// $Id: tmatrix.cc,v 1.7 1999-05-18 12:58:24 ansari Exp $
+// $Id: tmatrix.cc,v 1.8 1999-05-18 17:19:36 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -5 +5 @@
 #include <stdlib.h>
 #include <iostream.h>
+#include <values.h>
 #include <complex>
 #include "pexceptions.h"
@@ -162 +163 @@
 }
 
+////////////////////////////////////////////////////////////////
+//**** Pour gestion des lignes et des colonnes
+
+template <class T> TMatrixRC<T> TMatrix<T>::Row(uint_4 r) const
+{
+TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixRow, r);
+return rc;
+}
+
+template <class T> TMatrixRC<T> TMatrix<T>::Col(uint_4 c) const
+{
+TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixCol, c);
+return rc;
+}
+
+template <class T> TMatrixRC<T> TMatrix<T>::Diag() const
+{
+TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixDiag);
+return rc;
+}
+
 #include "matrix.h"
 ////////////////////////////////////////////////////////////////
 //**** 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
+
 r_8 TMatrix<r_8>::GausPiv(TMatrix<r_8>& a, TMatrix<r_8>& b)
 // Pivot de Gauss
 // * Attention: egcs impose que cette fonction soit mise dans le .cc
 //              avant ::Inverse() (car Inverse() l'utilise)
-{
-Matrix A(a);
-Matrix B(b);
-return (r_8) Matrix::GausPiv(A,B);
+// {Matrix A(a); Matrix B(b); return (r_8) Matrix::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<r_8>::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 << "Matrix warning, overflow for det" << endl;
+  } else {
+    det = exp(ld);
+  }
+}
+
+TMatrixRC<r_8> pivRowa(a,TmatrixRow);
+TMatrixRC<r_8> pivRowb(b,TmatrixRow);
+
+for(uint_4 k=0; k<n-1; k++) {
+  uint_4 iPiv = a.Col(k).IMaxAbs(k);
+  if(iPiv != k) {
+    TMatrixRC<r_8> aIPiv(a.Row(iPiv));
+    TMatrixRC<r_8> aK(a.Row(k));
+    TMatrixRC<r_8>::Swap(aIPiv,aK);
+    TMatrixRC<r_8> bIPiv(b.Row(iPiv));
+    TMatrixRC<r_8> bK(b.Row(k));
+    TMatrixRC<r_8>::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;
+    a.Row(i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
+    b.Row(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;
+    a.Row(jj).LinComb(r, pivRowa);
+    b.Row(jj).LinComb(r, pivRowb);
+  }
+}
+
+for(uint_4 l=0; l<n; l++) {
+  if (fabs(a(l,l)) <= 1.e-50) return 0.0;
+  b.Row(l) /= a(l,l);
+}
+
+return det;
 }
 
@@ -178 +282 @@
 // Inversion
 {
-TMatrix<r_8> b(mNc,mNr);
-TMatrix<r_8> a(*this);
-b = 1.;
-if (fabs(TMatrix<r_8>::GausPiv(a,b)) < 1.e-50)
+TMatrix<r_8> b(mNc,mNr); TMatrix<r_8> a(*this); b = 1.;
+if(fabs(TMatrix<r_8>::GausPiv(a,b)) < 1.e-50)
   throw(MathExc("TMatrix Inverse() Singular Matrix"));
 return b;
@@ -241 +343 @@
 }
 
-////////////////////////////////////////////////////////////////
-// -------------------------------------------------------------------------
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
 //   Les objets delegues pour la gestion de persistance 
-// -------------------------------------------------------------------------
+// --------------------------------------------------------
+///////////////////////////////////////////////////////////
 
 template <class T>
@@ -316 +419 @@
 FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
 fio_nd.Write(os);
+}
+
+////////////////////////////////////////////////////////////////
+// -------------------------------------------------------------
+//   La classe de gestion des lignes et colonnes d'une matrice
+// -------------------------------------------------------------
+////////////////////////////////////////////////////////////////
+#include "tvector.h"
+
+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.mNc != m.mNr)
+  throw(SzMismatchError("TMatrixRC::TMatrixRC(...,TmatrixDiag,...) size mismatch\n"));
+}
+
+template <class T> int_4 TMatrixRC<T>::Next()
+{
+if (!matrix || kind==TmatrixDiag) return -1;
+index++;
+if(kind == TmatrixRow) {
+  if(index > (int_4)matrix->mNr) {index = (int_4)matrix->mNr; return -1;}
+  data += matrix->mNc;
+} else {
+  if (index > (int_4)matrix->mNc) {index = (int_4)matrix->mNc; 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->mNc;
+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->mNc) 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->mNr) 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->mNc != matrix->mNr)
+  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;}
 }
 
@@ -343 +634 @@
 #pragma define_template FIO_TMatrix< complex<float> >
 #pragma define_template FIO_TMatrix< complex<double> >
+// 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<float> > 
+#pragma define_template TMatrixRC< complex<double> >
 #endif
 
@@ -369 +672 @@
 template class FIO_TMatrix< complex<float> >;
 template class FIO_TMatrix< complex<double> >;
+// 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<float> >;
+template class TMatrixRC< complex<double> >;
 #endif