Skip to content

Feature/#778 failing pseudo inverse

Elias Pipping requested to merge pipping/dune-geometry:feature/#778-failing-pseudo-inverse into master

This is an attempt to fix flyspray/FS#778 and flyspray/FS#866.

3x3 matrices can now be inverted via rightInvA() that would previously yield an error (this is not tested anywhere).

2x3 matrices can now be pseudo-inverted via rightInvA() that would previously yield an error (this is tested in the testmatrix test). What is the quality of said 3x2 inverse (called invB)? Here's the output that the test produces for me:

2: Running tests for type: double
2: sqrtDetAAT = 2.22045e-16
2: detA       = 2.22045e-16
2: invA        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 0
2: A        = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity    A*invA*A = A    satisfied up to error: 0
2: sqrtDetBBT = 2.22045e-16
2: detB       = 2.22045e-16
2: invB        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 0
2: B        = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity    B*invB*B = B    satisfied up to error: 0
2: 
2: Running tests for type: long double
2: sqrtDetAAT = 2.22045e-16
2: detA       = 2.22045e-16
2: invA        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 0
2: A        = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity    A*invA*A = A    satisfied up to error: 0
2: sqrtDetBBT = 2.22045e-16
2: detB       = 2.22045e-16
2: invB        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 0
2: B        = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity    B*invB*B = B    satisfied up to error: 0
2: 
2: Running tests for type: GMPField<72>
2: sqrtDetAAT = 2.22045e-16
2: detA       = 2.22045e-16
2: invA        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 1.50463e-38
2: A        = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity    A*invA*A = A    satisfied up to error: 1.61633e-41
2: sqrtDetBBT = 2.22045e-16
2: detB       = 2.22045e-16
2: invB        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 1.17549e-38
2: B        = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity    B*invB*B = B    satisfied up to error: 1.61633e-41
2: 
2: Running tests for type: GMPField<160>
2: sqrtDetAAT = 2.22045e-16
2: detA       = 2.22045e-16
2: invA        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: invA*A*invA = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14 ]
2: Identity invA*A*invA = invA satisfied up to error: 8.44269e-58
2: A        = [ 0.1 -0.01, 0.1 -0.01 ]
2: A*invA*A = [ 0.1 -0.01, 0.1 -0.01 ]
2: Identity    A*invA*A = A    satisfied up to error: 1.20888e-60
2: sqrtDetBBT = 2.22045e-16
2: detB       = 2.22045e-16
2: invB        = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: invB*B*invB = [ -4.5036e+13 4.5036e+13, -4.5036e+14 4.5036e+14, 0 0 ]
2: Identity invB*B*invB = invB satisfied up to error: 6.37237e-58
2: B        = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: B*invB*B = [ 0.1 -0.01 0, 0.1 -0.01 0 ]
2: Identity    B*invB*B = B    satisfied up to error: 1.20888e-60

(Update: redid the 2x3 case. Relying on maple was not a good idea. The generated code was unstable and incomprehensible. Now it behaves just as intended.)

In summary: the 2x3 pseudo-inverse is perfectly usable for the test case, also with the double field type.

Merge request reports