Skip to content
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

Fix xORBDB6 performs numerically unadvisable operations #647

Conversation

christoph-conrads
Copy link
Contributor

Description

fixes #634

Checklist

  • The documentation has been updated.
  • If the PR solves a specific issue, it is set to be closed on merge.

@christoph-conrads
Copy link
Contributor Author

The fixes are only implemented for SORBDB6 so far.

langou
langou previously approved these changes Apr 12, 2022
@christoph-conrads
Copy link
Contributor Author

@langou The changes from SORBDB6 were applied to DORBDB6, CUNBDB6, ZUNBDB6. The pull request should be complete.

@christoph-conrads
Copy link
Contributor Author

@langou There are 28 compiler errors and at least one test executable is not present. Is this fixed on master (the tests build and pass for me)? If not, do you know where can I see the error messages?

SRC/zunbdb6.f Outdated Show resolved Hide resolved
The zeros are overwritten by the following SORGQR call.
This patch was authored by Brian D. Sutton and posted to the discussion
of LAPACK pull request Reference-LAPACK#406.

* fix indexing for vector increments different from one
* always set vectors that are numerically zero to zero

Previously SORBDB6 would only set vectors to zero if a second iteration
of Gram-Schmidt was necessary. This would cause problems on the caller
site if the test for a zero vector differed from the SORBDB6 test for
zero.
* Require unit-norm vector X for otherwise the following computations
  might underflow
* Avoid over- and underflows in the computation of the Euclidean norm of
  X
* Fix the Euclidean norm computation after the second Gram-Schmidt
  iteration
* Consider round-off errors when checking for zero vectors
* Update identifiers
This change makes it easier to port the recent changes to SORBDB6 to
(complex) double precision.
* Require unit-norm vector X for otherwise the following computations
  might underflow
* Avoid over- and underflows in the computation of the Euclidean norm of
  X
* Fix the Euclidean norm computation after the second Gram-Schmidt
  iteration
* Consider round-off errors when checking for zero vectors
* Update identifiers

Note that the caller DORBDB5 always passed unit-norm vectors X.
* Require unit-norm vector X for otherwise the following computations
  might underflow
* Avoid over- and underflows in the computation of the Euclidean norm of
  X
* Fix the Euclidean norm computation after the second Gram-Schmidt
  iteration
* Consider round-off errors when checking for zero vectors
* Update identifiers

Note that the caller CUNBDB5 always passed unit-norm vectors X.
* Require unit-norm vector X for otherwise the following computations
  might underflow
* Avoid over- and underflows in the computation of the Euclidean norm of
  X
* Fix the Euclidean norm computation after the second Gram-Schmidt
  iteration
* Consider round-off errors when checking for zero vectors
* Update identifiers

Note that the caller ZUNBDB5 always passed unit-norm vectors X.

Thank you @angsch for discovering typos in the function names (`CLASSQ`
was called instead of `ZLASSQ`).
@christoph-conrads christoph-conrads force-pushed the 634-xORBDB6-performs-numerically-unadvisable-operations branch from 9a68612 to 54b3964 Compare July 10, 2022 20:17
@christoph-conrads
Copy link
Contributor Author

christoph-conrads commented Jul 14, 2022

The code changes in all four variants are horribly broken: the variable NORMSQ2 is nowhere introduced, nowhere being assigned to but it is being read. gfortran 11 does not complain about this without -Wextra flag. I would love to know the rationale for not warning about such obvious nonsense automatically.

      IF( NORMSQ2 .LE. N * EPS * NORM ) THEN

There is no warning with GCC 11 about an undeclared and uninitialized
varible being read without the `-Wextra` flag. Why?
@codecov
Copy link

codecov bot commented Jul 14, 2022

Codecov Report

Merging #647 (4842829) into master (4f97df9) will not change coverage.
The diff coverage is 0.00%.

❗ Current head 4842829 differs from pull request most recent head 145b1c3. Consider uploading reports for the commit 145b1c3 to get more accurate results

@@           Coverage Diff           @@
##           master     #647   +/-   ##
=======================================
  Coverage    0.00%    0.00%           
=======================================
  Files        1894     1894           
  Lines      184140   184118   -22     
=======================================
+ Misses     184140   184118   -22     
Impacted Files Coverage Δ
SRC/cunbdb6.f 0.00% <0.00%> (ø)
SRC/dorbdb6.f 0.00% <0.00%> (ø)
SRC/sorbdb2.f 0.00% <ø> (ø)
SRC/sorbdb4.f 0.00% <ø> (ø)
SRC/sorbdb6.f 0.00% <0.00%> (ø)
SRC/sorcsd2by1.f 0.00% <0.00%> (ø)
SRC/zunbdb6.f 0.00% <0.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4f97df9...145b1c3. Read the comment docs.

*> criterion, then the zero vector is returned.
*> The projection is computed with at most two iterations of the
*> classical Gram-Schmidt algorithm, see
*> * L. Giraud, J. Langou, M. Rozložník. "On the round-off error
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for quoting this paper.

I think the following paper is a better quote for this fact:

L. Giraud, J. Langou, M. Rozložník, Jasper van den Eshof. "Rounding error analysis of the classical Gram-Schmidt orthogonalization process." Numerische Mathematik, volume volume 101, pages 87–100, 2005. DOI = 10.1007/s00211-005-0615-4.

@langou langou merged commit 8b3bbfa into Reference-LAPACK:master Jul 15, 2022
@julielangou julielangou added this to the LAPACK 3.11.0 milestone Nov 12, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

xORBDB6 performs numerically unadvisable operations
4 participants