Changeset 1007 in Sophya for trunk/SophyaLib/TArray/sopemtx.cc

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

File:

• trunk/SophyaLib/TArray/sopemtx.cc (modified) (10 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
 
-
-
-
-
-
-
-
-
-
-
-
-