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.
|
Do you think we could add a test that checks that this continues to work? |
|
@peterrum Ping? |
1 similar comment
|
@peterrum Ping? |
|
I have updated the comment.
What do you mean? I have added a new test. |
... and early return if norm is zero. This might happen if all degrees of freedoms (on a coarse grid/single cell) are constrained.