Skip to content

Make vector orthogonalization more reliable #930

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.

Sign up for GitHub

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

Merged
merged 6 commits into from
Nov 11, 2023
Merged

Make vector orthogonalization more reliable #930

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
bc2b374
xORBDB6/xUNBDB6: do not require unit-norm vector
christoph-conrads Nov 7, 2023
ad3d9b8
xORBDB6/xUNBDB6: use safe termination criterion
christoph-conrads Nov 9, 2023
80d7c69
xORBDB5/xUNBDB5: introduce real zero parameter
christoph-conrads Nov 9, 2023
7eb4057
xORBDB5/xUNBDB5: detect numerically zero vectors
christoph-conrads Nov 9, 2023
5b0bc5e
xORBDB5/xUNBDB5: ensure xORBDB6 input has unit norm
christoph-conrads Nov 10, 2023
64897c4
xLARFGP: avoid overflows
christoph-conrads Nov 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 27 additions & 14 deletions SRC/cunbdb5.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,18 +169,21 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
* =====================================================================
*
* .. Parameters ..
REAL REALZERO
PARAMETER ( REALZERO = 0.0E0 )
COMPLEX ONE, ZERO
PARAMETER ( ONE = (1.0E0,0.0E0), ZERO = (0.0E0,0.0E0) )
* ..
* .. Local Scalars ..
INTEGER CHILDINFO, I, J
REAL EPS, NORM, SCL, SSQ
* ..
* .. External Subroutines ..
EXTERNAL CUNBDB6, XERBLA
EXTERNAL CLASSQ, CUNBDB6, XERBLA
* ..
* .. External Functions ..
REAL SCNRM2
EXTERNAL SCNRM2
REAL SLAMCH, SCNRM2
EXTERNAL SLAMCH, SCNRM2
* ..
* .. Intrinsic Function ..
INTRINSIC MAX
Expand Down Expand Up @@ -213,16 +216,26 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
RETURN
END IF
*
* Project X onto the orthogonal complement of Q
EPS = SLAMCH( 'Precision' )
*
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2,
$ WORK, LWORK, CHILDINFO )
* Project X onto the orthogonal complement of Q if X is nonzero
*
* If the projection is nonzero, then return
SCL = REALZERO
SSQ = REALZERO
CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
NORM = SCL * SQRT( SSQ )
*
IF( SCNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. SCNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
RETURN
IF( NORM .GT. N * EPS ) THEN
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
*
* If the projection is nonzero, then return
*
IF( SCNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. SCNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END IF
*
* Project each standard basis vector e_1,...,e_M1 in turn, stopping
Expand All @@ -238,8 +251,8 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
END DO
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
IF( SCNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. SCNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
IF( SCNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. SCNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END DO
Expand All @@ -257,8 +270,8 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
X2(I) = ONE
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
IF( SCNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. SCNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
IF( SCNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. SCNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END DO
Expand Down
21 changes: 11 additions & 10 deletions SRC/cunbdb6.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,9 +41,8 @@
*> with respect to the columns of
*> Q = [ Q1 ] .
*> [ Q2 ]
*> The Euclidean norm of X must be one and the columns of Q must be
*> orthonormal. The orthogonalized vector will be zero if and only if it
*> lies entirely in the range of Q.
*> The columns of Q must be orthonormal. The orthogonalized vector will
*> be zero if and only if it lies entirely in the range of Q.
*>
*> The projection is computed with at most two iterations of the
*> classical Gram-Schmidt algorithm, see
Expand Down Expand Up @@ -174,7 +173,7 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
*
* .. Parameters ..
REAL ALPHA, REALONE, REALZERO
PARAMETER ( ALPHA = 0.1E0, REALONE = 1.0E0,
PARAMETER ( ALPHA = 0.83E0, REALONE = 1.0E0,
$ REALZERO = 0.0E0 )
COMPLEX NEGONE, ONE, ZERO
PARAMETER ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0),
Expand Down Expand Up @@ -223,14 +222,16 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
*
EPS = SLAMCH( 'Precision' )
*
* Compute the Euclidean norm of X
*
SCL = REALZERO
SSQ = REALZERO
CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
NORM = SCL * SQRT( SSQ )
*
* First, project X onto the orthogonal complement of Q's column
* space
*
* Christoph Conrads: In debugging mode the norm should be computed
* and an assertion added comparing the norm with one. Alas, Fortran
* never made it into 1989 when assert() was introduced into the C
* programming language.
NORM = REALONE
*
IF( M1 .EQ. 0 ) THEN
DO I = 1, N
Expand Down
41 changes: 27 additions & 14 deletions SRC/dorbdb5.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,18 +169,21 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION REALZERO
PARAMETER ( REALZERO = 0.0D0 )
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
* ..
* .. Local Scalars ..
INTEGER CHILDINFO, I, J
DOUBLE PRECISION EPS, NORM, SCL, SSQ
* ..
* .. External Subroutines ..
EXTERNAL DORBDB6, XERBLA
EXTERNAL DLASSQ, DORBDB6, XERBLA
* ..
* .. External Functions ..
DOUBLE PRECISION DNRM2
EXTERNAL DNRM2
DOUBLE PRECISION DLAMCH, DNRM2
EXTERNAL DLAMCH, DNRM2
* ..
* .. Intrinsic Function ..
INTRINSIC MAX
Expand Down Expand Up @@ -213,16 +216,26 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
RETURN
END IF
*
* Project X onto the orthogonal complement of Q
EPS = DLAMCH( 'Precision' )
*
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2,
$ WORK, LWORK, CHILDINFO )
* Project X onto the orthogonal complement of Q if X is nonzero
*
* If the projection is nonzero, then return
SCL = REALZERO
SSQ = REALZERO
CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
NORM = SCL * SQRT( SSQ )
*
IF( DNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
RETURN
IF( NORM .GT. N * EPS ) THEN
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
*
* If the projection is nonzero, then return
*
IF( DNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. DNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END IF
*
* Project each standard basis vector e_1,...,e_M1 in turn, stopping
Expand All @@ -238,8 +251,8 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
END DO
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
IF( DNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
IF( DNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. DNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END DO
Expand All @@ -257,8 +270,8 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
X2(I) = ONE
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
IF( DNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
IF( DNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. DNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END DO
Expand Down
21 changes: 11 additions & 10 deletions SRC/dorbdb6.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,9 +41,8 @@
*> with respect to the columns of
*> Q = [ Q1 ] .
*> [ Q2 ]
*> The Euclidean norm of X must be one and the columns of Q must be
*> orthonormal. The orthogonalized vector will be zero if and only if it
*> lies entirely in the range of Q.
*> The columns of Q must be orthonormal. The orthogonalized vector will
*> be zero if and only if it lies entirely in the range of Q.
*>
*> The projection is computed with at most two iterations of the
*> classical Gram-Schmidt algorithm, see
Expand Down Expand Up @@ -174,7 +173,7 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
*
* .. Parameters ..
DOUBLE PRECISION ALPHA, REALONE, REALZERO
PARAMETER ( ALPHA = 0.1D0, REALONE = 1.0D0,
PARAMETER ( ALPHA = 0.83D0, REALONE = 1.0D0,
$ REALZERO = 0.0D0 )
DOUBLE PRECISION NEGONE, ONE, ZERO
PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
Expand Down Expand Up @@ -222,14 +221,16 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
*
EPS = DLAMCH( 'Precision' )
*
* Compute the Euclidean norm of X
*
SCL = REALZERO
SSQ = REALZERO
CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
NORM = SCL * SQRT( SSQ )
*
* First, project X onto the orthogonal complement of Q's column
* space
*
* Christoph Conrads: In debugging mode the norm should be computed
* and an assertion added comparing the norm with one. Alas, Fortran
* never made it into 1989 when assert() was introduced into the C
* programming language.
NORM = REALONE
*
IF( M1 .EQ. 0 ) THEN
DO I = 1, N
Expand Down
41 changes: 27 additions & 14 deletions SRC/sorbdb5.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,18 +169,21 @@ SUBROUTINE SORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
* =====================================================================
*
* .. Parameters ..
REAL REALZERO
PARAMETER ( REALZERO = 0.0E0 )
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E0, ZERO = 0.0E0 )
* ..
* .. Local Scalars ..
INTEGER CHILDINFO, I, J
REAL EPS, NORM, SCL, SSQ
* ..
* .. External Subroutines ..
EXTERNAL SORBDB6, XERBLA
EXTERNAL SLASSQ, SORBDB6, XERBLA
* ..
* .. External Functions ..
REAL SNRM2
EXTERNAL SNRM2
REAL SLAMCH, SNRM2
EXTERNAL SLAMCH, SNRM2
* ..
* .. Intrinsic Function ..
INTRINSIC MAX
Expand Down Expand Up @@ -213,16 +216,26 @@ SUBROUTINE SORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
RETURN
END IF
*
* Project X onto the orthogonal complement of Q
EPS = SLAMCH( 'Precision' )
*
CALL SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2,
$ WORK, LWORK, CHILDINFO )
* Project X onto the orthogonal complement of Q if X is nonzero
*
* If the projection is nonzero, then return
SCL = REALZERO
SSQ = REALZERO
CALL SLASSQ( M1, X1, INCX1, SCL, SSQ )
CALL SLASSQ( M2, X2, INCX2, SCL, SSQ )
NORM = SCL * SQRT( SSQ )
*
IF( SNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. SNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
RETURN
IF( NORM .GT. N * EPS ) THEN
CALL SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
*
* If the projection is nonzero, then return
*
IF( SNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. SNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END IF
*
* Project each standard basis vector e_1,...,e_M1 in turn, stopping
Expand All @@ -238,8 +251,8 @@ SUBROUTINE SORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
END DO
CALL SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
IF( SNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. SNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
IF( SNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. SNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END DO
Expand All @@ -257,8 +270,8 @@ SUBROUTINE SORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
X2(I) = ONE
CALL SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
$ LDQ2, WORK, LWORK, CHILDINFO )
IF( SNRM2(M1,X1,INCX1) .NE. ZERO
$ .OR. SNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
IF( SNRM2(M1,X1,INCX1) .NE. REALZERO
$ .OR. SNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
RETURN
END IF
END DO
Expand Down
21 changes: 11 additions & 10 deletions SRC/sorbdb6.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,9 +41,8 @@
*> with respect to the columns of
*> Q = [ Q1 ] .
*> [ Q2 ]
*> The Euclidean norm of X must be one and the columns of Q must be
*> orthonormal. The orthogonalized vector will be zero if and only if it
*> lies entirely in the range of Q.
*> The columns of Q must be orthonormal. The orthogonalized vector will
*> be zero if and only if it lies entirely in the range of Q.
*>
*> The projection is computed with at most two iterations of the
*> classical Gram-Schmidt algorithm, see
Expand Down Expand Up @@ -174,7 +173,7 @@ SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
*
* .. Parameters ..
REAL ALPHA, REALONE, REALZERO
PARAMETER ( ALPHA = 0.1E0, REALONE = 1.0E0,
PARAMETER ( ALPHA = 0.83E0, REALONE = 1.0E0,
$ REALZERO = 0.0E0 )
REAL NEGONE, ONE, ZERO
PARAMETER ( NEGONE = -1.0E0, ONE = 1.0E0, ZERO = 0.0E0 )
Expand Down Expand Up @@ -222,14 +221,16 @@ SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
*
EPS = SLAMCH( 'Precision' )
*
* Compute the Euclidean norm of X
*
SCL = REALZERO
SSQ = REALZERO
CALL SLASSQ( M1, X1, INCX1, SCL, SSQ )
CALL SLASSQ( M2, X2, INCX2, SCL, SSQ )
NORM = SCL * SQRT( SSQ )
*
* First, project X onto the orthogonal complement of Q's column
* space
*
* Christoph Conrads: In debugging mode the norm should be computed
* and an assertion added comparing the norm with one. Alas, Fortran
* never made it into 1989 when assert() was introduced into the C
* programming language.
NORM = REALONE
*
IF( M1 .EQ. 0 ) THEN
DO I = 1, N
Expand Down
Loading