-
Notifications
You must be signed in to change notification settings - Fork 46
Sparse matrix error: detect ill conditions on-the-fly #1205
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Conversation
Signed-off-by: Martijn Govers <Martijn.Govers@Alliander.com>
| .binaryExpr( | ||
| (matrix.col(pivot).tail(size - pivot - 1) * matrix.row(pivot).tail(size - pivot - 1)), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
TODO: also add the contribution from the accumulated_error.col(pivot).tail(size - pivot - 1) * matrix.row(pivot).tail(size - pivot - 1)) and (matrix.col(pivot).tail(size - pivot - 1) * accumulated_error.row(pivot).tail(size - pivot - 1)
| // check stability | ||
| Scalar const piv_size = cabs(matrix(pivot, pivot)); | ||
| double const piv_error = accumulated_error(pivot, pivot); | ||
| if (cabs(matrix(pivot, pivot)) <= std::max(accumulated_error.row(pivot).tail(size - pivot).maxCoeff(), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
this is actually a very weak criterion. by the time your accumulated error even becomes close to your pivot element size, the rest of your calculation may already become garbage. E.g., if you want a precision of 1E-5, odds are that your calculation will not be as precise anymore when your matrix element precision is something like 1E-1. This criterion probably should be stricter:
| if (cabs(matrix(pivot, pivot)) <= std::max(accumulated_error.row(pivot).tail(size - pivot).maxCoeff(), | |
| if (cabs(matrix(pivot, pivot)) <= threshold * std::max(accumulated_error.row(pivot).tail(size - pivot).maxCoeff(), |
Signed-off-by: Martijn Govers <martygovers@hotmail.com> Signed-off-by: Martijn Govers <Martijn.Govers@Alliander.com> Apply suggestion from @mgovers Signed-off-by: Martijn Govers <martygovers@hotmail.com>
888321e to
912af80
Compare
This is an experiment. Do not merge yet!
Rather than using some form of condition numbers to detect unsolvable matrix decompositions, we now detect them on-the-fly
NOTE: some ill-conditioned cases now actually fail even though they did not do so before!!! It is dependent on the calculation method and whether the calculation is symmetric or asymmetric whether the system raises errors or not.