[Getfem-commits] (no subject)

[Tetsuo Koyama]

branch: devel-tetsuo-fix-error
commit d57b0188ec3f34038a4139f83bb3eaf70f21fe46
Author: Tetsuo Koyama <address@hidden>
Date:   Thu Jan 31 12:53:21 2019 +0900

    Fix test error
---
 src/gmm/gmm_opt.h | 101 +++++++++++++++++++++++++++---------------------------
 1 file changed, 51 insertions(+), 50 deletions(-)

diff --git a/src/gmm/gmm_opt.h b/src/gmm/gmm_opt.h
index 87b1929..0d2d0b2 100644
--- a/src/gmm/gmm_opt.h
+++ b/src/gmm/gmm_opt.h
@@ -55,16 +55,16 @@ namespace gmm {
       case 2 : return (*p) * (*(p+3)) - (*(p+1)) * (*(p+2));
 // Not stable for nearly singular matrices
 //       case 3 : return (*p) * ((*(p+4)) * (*(p+8)) - (*(p+5)) * (*(p+7)))
-//              - (*(p+1)) * ((*(p+3)) * (*(p+8)) - (*(p+5)) * (*(p+6)))
-//              + (*(p+2)) * ((*(p+3)) * (*(p+7)) - (*(p+4)) * (*(p+6)));
+//                  - (*(p+1)) * ((*(p+3)) * (*(p+8)) - (*(p+5)) * (*(p+6)))
+//                  + (*(p+2)) * ((*(p+3)) * (*(p+7)) - (*(p+4)) * (*(p+6)));
       default :
-       {
-         dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
-         lapack_ipvt ipvt(mat_nrows(A));
-         gmm::copy(A, B);
-         lu_factor(B, ipvt);
-         return lu_det(B, ipvt);       
-       }
+        {
+          dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
+          lapack_ipvt ipvt(mat_nrows(A));
+          gmm::copy(A, B);
+          lu_factor(B, ipvt);
+          return lu_det(B, ipvt);
+        }
       }
     }
     return T(1);
@@ -72,55 +72,56 @@ namespace gmm {
 
 
   template <typename T> T lu_inverse(const dense_matrix<T> &A_,
-                                    bool doassert = true) {
+                                     bool doassert = true) {
     dense_matrix<T>& A = const_cast<dense_matrix<T> &>(A_);
     size_type N = mat_nrows(A);
     T det(1);
     if (N) {
       T *p = &(A(0,0));
-       switch (N) {
-         case 1 : {
-           det = *p;
-           if (doassert) GMM_ASSERT1(det!=T(0), "non invertible matrix");
-            if (det == T(0)) break;
-           *p = T(1) / det; 
-         } break;
-         case 2 : {
-           det = (*p) * (*(p+3)) - (*(p+1)) * (*(p+2));
-           if (doassert) GMM_ASSERT1(det!=T(0), "non invertible matrix");
-            if (det == T(0)) break;
-           std::swap(*p, *(p+3));
-           *p++ /= det; *p++ /= -det; *p++ /= -det; *p++ /= det; 
-         } break;
-         case 3 : { // not stable for nearly singular matrices
-           T a, b, c, d, e, f, g, h, i;
-           a =   (*(p+4)) * (*(p+8)) - (*(p+5)) * (*(p+7));
-           b = - (*(p+1)) * (*(p+8)) + (*(p+2)) * (*(p+7));
-           c =   (*(p+1)) * (*(p+5)) - (*(p+2)) * (*(p+4));
-           d = - (*(p+3)) * (*(p+8)) + (*(p+5)) * (*(p+6));
-           e =   (*(p+0)) * (*(p+8)) - (*(p+2)) * (*(p+6));
-           f = - (*(p+0)) * (*(p+5)) + (*(p+2)) * (*(p+3));
-           g =   (*(p+3)) * (*(p+7)) - (*(p+4)) * (*(p+6));
-           h = - (*(p+0)) * (*(p+7)) + (*(p+1)) * (*(p+6));
-           i =   (*(p+0)) * (*(p+4)) - (*(p+1)) * (*(p+3));
-           det = (*p) * a + (*(p+1)) * d + (*(p+2)) * g;
-        if (std::abs(det) > 1e-5) {
-             *p++ = a / det; *p++ = b / det; *p++ = c / det;
-             *p++ = d / det; *p++ = e / det; *p++ = f / det;
-             *p++ = g / det; *p++ = h / det; *p++ = i / det;
+      switch (N) {
+        case 1 : {
+          det = *p;
+          if (doassert) GMM_ASSERT1(det!=T(0), "non invertible matrix");
+          if (det == T(0)) break;
+          *p = T(1) / det; 
+        } break;
+        case 2 : {
+          det = (*p) * (*(p+3)) - (*(p+1)) * (*(p+2));
+          if (doassert) GMM_ASSERT1(det!=T(0), "non invertible matrix");
+          if (det == T(0)) break;
+          std::swap(*p, *(p+3));
+          *p++ /= det; *p++ /= -det; *p++ /= -det; *p++ /= det; 
+        } break;
+        default : {
+          if (N == 3) { // not stable for nearly singular matrices
+            T a, b, c, d, e, f, g, h, i;
+            a =   (*(p+4)) * (*(p+8)) - (*(p+5)) * (*(p+7));
+            b = - (*(p+1)) * (*(p+8)) + (*(p+2)) * (*(p+7));
+            c =   (*(p+1)) * (*(p+5)) - (*(p+2)) * (*(p+4));
+            d = - (*(p+3)) * (*(p+8)) + (*(p+5)) * (*(p+6));
+            e =   (*(p+0)) * (*(p+8)) - (*(p+2)) * (*(p+6));
+            f = - (*(p+0)) * (*(p+5)) + (*(p+2)) * (*(p+3));
+            g =   (*(p+3)) * (*(p+7)) - (*(p+4)) * (*(p+6));
+            h = - (*(p+0)) * (*(p+7)) + (*(p+1)) * (*(p+6));
+            i =   (*(p+0)) * (*(p+4)) - (*(p+1)) * (*(p+3));
+            det = (*p) * a + (*(p+1)) * d + (*(p+2)) * g;
+            if (std::abs(det) > 1e-5) {
+              *p++ = a / det; *p++ = b / det; *p++ = c / det;
+              *p++ = d / det; *p++ = e / det; *p++ = f / det;
+              *p++ = g / det; *p++ = h / det; *p++ = i / det;
+              break;
+            }
+          }
+          dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
+          lapack_ipvt ipvt(mat_nrows(A));
+          gmm::copy(A, B);
+          size_type info = lu_factor(B, ipvt);
+          GMM_ASSERT1(!info, "non invertible matrix");
+          lu_inverse(B, ipvt, A);
+          det = lu_det(B, ipvt);
           break;
         }
-         }
-      default : {
-           dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
-           lapack_ipvt ipvt(mat_nrows(A));
-           gmm::copy(A, B);
-       size_type info = lu_factor(B, ipvt);
-           GMM_ASSERT1(!info, "non invertible matrix");
-           lu_inverse(B, ipvt, A);
-           return lu_det(B, ipvt);
       }
-       }
     }
     return det;
   }