[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Markus Bürg |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Tue, 22 Aug 2017 09:02:28 -0400 (EDT) |
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;
}