Changeset 1007 in Sophya

Timestamp:

May 16, 2000, 7:54:20 PM (25 years ago)

Author:

ansari

Message:

gestion det dans GausPiv et new norm cmv 16/5/00

Location:

trunk/SophyaLib/TArray

Files:

• sopemtx.cc (modified) (10 diffs)
• sopemtx.h (modified) (6 diffs)

trunk/SophyaLib/TArray/sopemtx.cc

@@ -438 +438 @@
 #define M_LN2 0.69314718055994530942
 #endif
+//// CMV BUG BUG : sur OSF 5.0 DMINEXP est deconnant (~1.e+19 !!!)
 #ifndef LN_MINDOUBLE
 #define LN_MINDOUBLE  (M_LN2 * (DMINEXP - 1))
@@ -445 +446 @@
 #endif
 
+template <class T>
+int SimpleMatrixOperation<T>::GausPivScaling = PIV_GLOB_SCALE;
+
 //! Gaussian pivoting
 /*!
   Do Gauss pivoting of \b a, doing the same operations on matrix \b b
-  \return determinant of \b a.
+  \param computedet = true : return determimant of \b a (beware of over/underfloat)
+  \param computedet = false : return 1 if OK, 0 if not.
   \verbatim
   Solve linear system A(n,n) * X(n,m) = B(n,m)
@@ -457 +462 @@
  */  
 template <class T>
-T SimpleMatrixOperation<T>::GausPiv(TMatrix<T>& a, TMatrix<T>& b)
+T SimpleMatrixOperation<T>::GausPiv(TMatrix<T>& a, TMatrix<T>& b,bool computedet)
 // Pivot de Gauss
 // * Attention: egcs impose que cette fonction soit mise dans le .cc
@@ -466 +471 @@
 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 /= (T) nrm;
-  b /= (T) nrm;
-  //cout << "normalisation matrice " << nrm << endl;
-} else nrm=1;
 
 T det = 1;
-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 = (T) exp(ld);
-  }
-}
+
+//////////////////
+// Data scaling //
+//////////////////
+
+// Pas de normalisation des donnees
+if(GausPivScaling==PIV_NO_SCALE) {
+  if(computedet) det = (T) 1;
+// normalisation des donnees ligne par ligne
+} else if(GausPivScaling==PIV_ROW_SCALE) {
+  double nrm = 0.;
+  for(uint_4 iii=0; iii<a.NRows(); iii++) {
+    uint_4 jjj; double vmax = -1.;
+    for(jjj=0; jjj<a.NCols(); jjj++) {
+      double v = TMatrixRC<T>::Abs_Value(a(iii,jjj));
+      if(v>vmax) vmax = v;
+    }
+    if(vmax>0.) {
+      for(jjj=0; jjj<a.NCols(); jjj++) a(iii,jjj) /= (T) vmax;
+      for(jjj=0; jjj<b.NCols(); jjj++) b(iii,jjj) /= (T) vmax;
+      nrm += log(vmax);
+    } else return (T) 0;
+  }
+  if(computedet) {
+    if(nrm <= LN_MINDOUBLE || nrm >= LN_MAXDOUBLE) {
+      cerr<<"GausPiv_Row: normalisation failure, "
+          <<"determinant has to be multiplied by exp("<<nrm<<")"<<endl;
+    } else det = (T) exp(nrm);
+  }
+// On fait une normalisation un peu brutale globale
+} else {
+  double vmin=MAXDOUBLE, 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>0.) { a /= (T) nrm; b /= (T) nrm;} else return (T) 0;
+  if(computedet) {
+    nrm = a.NRows() * log(nrm);
+    if (nrm <= LN_MINDOUBLE || nrm >= LN_MAXDOUBLE) {
+      cerr<<"GausPiv_Glob: normalisation failure, "
+          <<"determinant has to be multiplied by exp("<<nrm<<")"<<endl;
+    } else det = (T) exp(nrm);
+  }
+}
+
+////////////////////////////////////////
+// Gaussian elimination with pivoting //
+////////////////////////////////////////
 
 TMatrixRC<T> pivRowa(a,TMatrixRC<T>::TmatrixRow);
@@ -507 +541 @@
   T pivot = a(k,k);
   if( TMatrixRC<T>::Abs_Value(pivot) < 1.e-50 ) return (T) 0;
-  //det *= pivot;
+  if(computedet) det *= pivot;
   pivRowa.SetRow(k); // to avoid constructors
   pivRowb.SetRow(k);
@@ -516 +550 @@
   }
 }
-det *= a(n-1, n-1);
+if(computedet) det *= a(n-1, n-1);
 
 // on remonte
@@ -545 +579 @@
 TMatrix<T> a(A,false);
 TMatrix<T> b(a.NCols(),a.NRows());  b = IdentityMatrix(1.);
-if( TMatrixRC<T>::Abs_Value(GausPiv(a,b)) < 1.e-50)
+if(GausPiv(a,b)==(T) 0)
   throw(MathExc("TMatrix Inverse() Singular Matrix"));
 b.SetTemp(true);
@@ -620 +654 @@
      a(j,k) = a(k,j) = TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), j) 
                      * TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), k);
-                                /* $CHECK$ Reza 10/3/2000 */
 
 c = fx * y;
 
 if(SimpleMatrixOperation<T>::GausPiv(a,c)==(T) 0) THROW(singMatxErr);
-                   /* $CHECK$ Reza 10/3/2000 */
 
 r_8 xi2 = 0., ax;
@@ -718 +750 @@
 
 if(SimpleMatrixOperation<T>::GausPiv(a,d)==(T) 0) THROW(singMatxErr);
-                                            /* $CHECK$ Reza 10/3/2000 */
 
 for(uint_4 l=0; l<nf; l++) {
@@ -806 +837 @@
 #endif
 
-
-
-
-
-
-
-
-
-
-
-
-

trunk/SophyaLib/TArray/sopemtx.h

@@ -27 +27 @@
 class SimpleMatrixOperation {
 public:
+  //! To define the type of data scaling for Gauss pivoting method (type T)
+  enum GausPivScal {
+    PIV_NO_SCALE = 0,   //!< do not scale datas for gauss pivoting
+    PIV_GLOB_SCALE = 1, //!< do a global scaling of datas for gauss pivoting (default)
+    PIV_ROW_SCALE = 2   //!< do a scaling by row of datas for gauss pivoting
+  };
+  static inline int SetGausPivScal(int gps = PIV_GLOB_SCALE)
+              { GausPivScaling = gps; return GausPivScaling; }
+
   static TMatrix<T> Inverse(TMatrix<T> const & A);
-  static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
+  static T GausPiv(TMatrix<T>& A, TMatrix<T>& B, bool computedet=false);
+
+protected:
+  static int GausPivScaling;
 };
 
@@ -47 +59 @@
 //------------------------------------------------------------
 /*! \ingroup TArray \fn LinSolveInPlace(TMatrix<T>&,TVector<T>&)
-    \brief  Solve A*C = B for C in place and return determinant
+    \brief  Solve A*C = B for C in place and return 1 if OK, 0 if not.
 */
 template <class T>
@@ -62 +74 @@
 /*! \ingroup TArray 
     \fn LinSolve(const TMatrix<T>&, const TVector<T>&, TVector<T>&)
-    \brief Solve A*C = B and return C and determinant
+    \brief Solve A*C = B and return C. Return 1 if OK, 0 if not.
 */
 template <class T>
@@ -86 +98 @@
 /*! \ingroup TArray 
     \fn Inverse(TMatrix<T> const &)
-    \brief To inverse a TMatrix
+    \brief To inverse a TMatrix.
 */
 template <class T>
@@ -94 +106 @@
 /*! \ingroup TArray 
     \fn Determinant(TMatrix<T> const &)
-    \brief Give the TMatrix determinant
+    \brief Return the TMatrix determinant
 */
 template <class T>
@@ -100 +112 @@
   TMatrix<T> a(A,false);
   TMatrix<T> b(a.NCols(),a.NRows());  b = IdentityMatrix(1.);
-  return SimpleMatrixOperation<T>::GausPiv(a,b);
+  return SimpleMatrixOperation<T>::GausPiv(a,b,true);
 }
 
 /*! \ingroup TArray 
-    \fn GausPiv(TMatrix<T> const &,TMatrix<T> &)
+    \fn GausPiv(TMatrix<T> const &,TMatrix<T> &, bool)
     \brief Gauss pivoting on matrix \b A, doing the same operations
-           on matrix \b B and return determinant of \b A.
-    \sa SOPHYA::SimpleMatrixOperation::GausPiv(TMatrix<T>&,TMatrix<T>&)
+           on matrix \b B.
+    \param computedet = true : return the determinant of \b A
+           (beware of over/underfloat), default is false.
+    \return determinant if requested, or 1 if OK.
+    \sa SOPHYA::SimpleMatrixOperation::GausPiv(TMatrix<T>&,TMatrix<T>&,bool)
 */
 template <class T>
-inline T GausPiv(TMatrix<T> const & A,TMatrix<T> & B) {
+inline T GausPiv(TMatrix<T> const & A,TMatrix<T> & B, bool computedet=false) {
   TMatrix<T> a(A,false);
-  return SimpleMatrixOperation<T>::GausPiv(a,B);
+  return SimpleMatrixOperation<T>::GausPiv(a,B,computedet);
 }