From b2547e01173bb4010887c1c0b3f49eaf687139ca Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Wed, 14 Dec 2011 08:55:24 +0000 Subject: [PATCH] Math inversion used to return zero (math_inv3x3) or error (math_invert3x3) for negative determinant. Now checking whether the absolute(!) value of the determinant is close to zero to avoid singularities, negative determinants are very well allowed. --- code/math.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/code/math.f90 b/code/math.f90 index 9fd1734a9..eb1146f1b 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -757,7 +757,7 @@ pure function math_transpose3x3(A) - A(1,2) * (A(2,1) * A(3,3) - A(2,3) * A(3,1))& + A(1,3) * (A(2,1) * A(3,2) - A(2,2) * A(3,1)) - if (DetA > tiny(DetA)) then + if (abs(DetA) > tiny(abs(DetA))) then math_inv3x3(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2)) / DetA math_inv3x3(2,1) = (-A(2,1) * A(3,3) + A(2,3) * A(3,1)) / DetA math_inv3x3(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1)) / DetA @@ -796,7 +796,7 @@ pure function math_transpose3x3(A) - A(1,2) * (A(2,1) * A(3,3) - A(2,3) * A(3,1))& + A(1,3) * (A(2,1) * A(3,2) - A(2,2) * A(3,1)) - if (DetA <= tiny(DetA)) then + if (abs(DetA) <= tiny(abs(DetA))) then error = .true. else InvA(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2)) / DetA