Power iteration: check norm of initial vector #18178
Conversation
include/deal.II/lac/precondition.h
Outdated
| eigenvector /= eigenvector.l2_norm(); | ||
|
|
||
| const auto norm = eigenvector.l2_norm(); | ||
|
|
||
| if (norm == 0.0) | ||
| return 1.0; | ||
|
|
||
| eigenvector /= norm; |
There was a problem hiding this comment.
I think that's the wrong solution. If you pass in a zero vector by accident, the function pretends that the largest eigenvalue of the operator is one. That's just not true.
The function should probably document that one must provide a non-zero vector and that if the vector is zero throw an exception that one can catch.
One should note that this is the internal eigenvalue estimation of In the new commit, I have moved the if-statement from This is the best solution I can come up. I am open for suggestions how to alternatively fix this issue. |
| if (data.eigenvalue_algorithm == internal::EigenvalueAlgorithm::lanczos) | ||
| // note: temp_vector1 only contains zeroes, which happens if all dofs | ||
| // are constrained, there is nothing to do | ||
| if (temp_vector1.l2_norm() != 0.0) |
There was a problem hiding this comment.
Is it possible to use all_zero() here? We do require that for the vector space vector concept (which this template doesn't check for).
There was a problem hiding this comment.
Indeed, that should be cheaper.
| if (data.eigenvalue_algorithm == internal::EigenvalueAlgorithm::lanczos) | ||
| // note: temp_vector1 only contains zeroes, which happens if all dofs | ||
| // are constrained, there is nothing to do | ||
| if (temp_vector1.l2_norm() != 0.0) |
There was a problem hiding this comment.
Indeed, that should be cheaper.
| // note: temp_vector1 only contains zeroes, which happens if all dofs | ||
| // are constrained, there is nothing to do | ||
| if (temp_vector1.l2_norm() != 0.0) |
There was a problem hiding this comment.
I think the comment has the wrong grammar. "temp_vector1 only contains zeroes" is a statement, but the statement is not generally true. I think what you want to say is something like this:
| // note: temp_vector1 only contains zeroes, which happens if all dofs | |
| // are constrained, there is nothing to do | |
| if (temp_vector1.l2_norm() != 0.0) | |
| // We want to run the power iteration starting with temp_vector1, but for that it | |
| // needs to be a nonzero vector. However, there are situations where it only contains | |
| // zeroes (e.g., if all dofs are constrained). If that is the case, there is nothing to do. | |
| if (temp_vector1.l2_norm() != 0.0) |
|
Do you think we could add a test that checks that this continues to work? |
|
@peterrum Ping? |
... and early return if norm is zero. This might happen if all degrees of freedoms (on a coarse grid/single cell) are constrained.