Compare the 4x4 cofactor determinant against its own scale - #7279
Conversation
GitPaean
added
the
manual:irrelevant
5d1481c to
9786636
Compare
| template <template<class K> class Matrix, typename K> | ||
| static inline K invertMatrix4(const Matrix<K>& matrix, Matrix<K>& inverse) | ||
| { | ||
| constexpr K singularLimit = K(1e-12); |
There was a problem hiding this comment.
the previous PR dropped a note that the gpu code needs to be synced. as this is now quite different, so @kjetilly @multitalentloes please comment.
There was a problem hiding this comment.
suggestion for both GPU and float, maybe it can be
constexpr K singularLimit = K(4.5e3) * std::numeric_limits<K>::epsilon();There was a problem hiding this comment.
With @multitalentloes 's comment, I think after the above suggestion goes in for the purpose of flow_blackoil_float, I think the PR can go in.
|
I think this is a genuine improvement and awaiting the people more familiar with the matter to review/discuss more. |
|
jenkins build this failure_report please |
9786636 to
eec27b8
Compare
|
I do not think it is necessary. We are waiting for some comments from the GPU development side. @kjetilly @multitalentloes |
Sorry for missing the earlier mention! |
invertMatrix4() tested |det| < 1e-40. The determinant is homogeneous of degree one in every row and in every column, so that measures the units a block is written in rather than how close to singular it is, and which blocks take the pivoted fallback moves when the model is rescaled. Compare it against max_i sum_j |a_ij * adj_ji| instead. Those products have the same homogeneity, so the ratio is left unchanged by any diagonal scaling A -> R*A*C, and each of the four sums is a cofactor expansion of the determinant, so the ratio never exceeds one. The block reported in OPM#7266 scores 0.20 and the well conditioned case in test_invert.cpp scores 0.66, so both take the cofactor path rather than the fallback. Formed from the inverse rather than the adjugate the same quantity is the reciprocal of that ratio, so the fallback result is judged by it as well. Dune reports a matrix as singular only for an exactly zero pivot, which a merely nearly dependent block does not produce; without the check a rank three block inverts to max |A*inv(A) - I| = 20 at ordinary scale and that inverse is returned silently. The four row sums are accumulated separately so the test does not sit on the critical path. It costs 4% on a single inversion and 22% streaming over blocks (clang -O3 -march=native, Apple M4 Pro), against the 10 to 20 4x4 block mat-muls ILU0 does per row for each inversion. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A rank three block, whose cofactor determinant comes out at roundoff rather than at zero and whose elimination reaches no exactly zero pivot, so that it is caught neither by an absolute threshold nor by the pivoting. That case fails against the previous code, which returned an inverse with max |A*inv(A) - I| = 20. And a well conditioned block together with a rank two block under row and column scalings spanning nearly two hundred orders of magnitude, which have to be inverted and rejected respectively at every one of them. The scalings are powers of two so that they round nothing of their own and so that the rank two block keeps its exactly zero pivot. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
eec27b8 to
54acaac
Compare
4 participants
Follow-up to #7266.
|det| < 1e-40tests the units a block is written in rather than its conditioning, so which blocks take the pivoted fallback moves when the model is rescaled. This compares the determinant againstmax_i sum_j |a_ij * adj_ji|instead, which has the same homogeneity in every row and column and so is invariant underA -> R*A*C. The block from #7266 scores 0.20 and the well conditioned case intest_invert.cppscores 0.66, so both take the cofactor path directly.Formed from the inverse rather than the adjugate the same quantity is the norm-relative guard offered in the Known trade-off of #7266, so the fallback result is judged by it too. That catches a rank three block, which on master inverts silently to
max |A*inv(A) - I| = 20at ordinary scale, no rescaling needed — the added test fails without the change.Cost is 4% on a single inversion and 22% streaming (clang -O3, M4 Pro), against the 10 to 20 4x4 block mat-muls ILU0 does per row, so a few percent of preconditioner setup. Not run on a deck — I have nothing that reproduces #7266.
Open: whether 1e-12 is the right limit, and whether the fallback still earns its place now that the trigger is scale free.