[Getfem-commits] (no subject)

[Markus Bürg]

branch: mb-transInversion
commit d538ee7eef289736026d0d2fca81cd07b5fc02f7
Author: mb <address@hidden>
Date:   Tue Aug 22 15:02:21 2017 +0200

    Use analytic formula, if matrix is well-conditioned.
---
 src/gmm/gmm_opt.h | 56 +++++++++++++++++++++++++++----------------------------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/src/gmm/gmm_opt.h b/src/gmm/gmm_opt.h
index e73af41..4dbe05f 100644
--- a/src/gmm/gmm_opt.h
+++ b/src/gmm/gmm_opt.h
@@ -75,7 +75,6 @@ namespace gmm {
     T det(1);
     if (N) {
       T *p = &(A(0,0));
-      if (N <= 2) {
        switch (N) {
          case 1 : {
            det = *p;
@@ -90,34 +89,35 @@ namespace gmm {
            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;
-//         GMM_ASSERT1(det!=T(0), "non invertible matrix");
-//         *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;
-       }
-      }
-      else {
-       dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
-       std::vector<int> 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);
+         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;
+          break;
+        }
+         }
+      default : {
+           dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
+           std::vector<int> 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;
   }