Skip to content

Commit 6261d62

Browse files
langoulangou
authored
Merge pull request #930 from christoph-conrads/917-sorcsd2by1-computes-inaccurate-result
Make vector orthogonalization more reliable
2 parents d469b91 + 64897c4 commit 6261d62

File tree

12 files changed

+193
-105
lines changed

12 files changed

+193
-105
lines changed

‎SRC/clarfgp.f‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
122122
* ..
123123
* .. Local Scalars ..
124124
INTEGER J, KNT
125-
REAL ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
125+
REAL ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM
126126
COMPLEX SAVEALPHA
127127
* ..
128128
* .. External Functions ..
@@ -143,11 +143,12 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
143143
RETURN
144144
END IF
145145
*
146+
EPS = SLAMCH( 'Precision' )
146147
XNORM = SCNRM2( N-1, X, INCX )
147148
ALPHR = REAL( ALPHA )
148149
ALPHI = AIMAG( ALPHA )
149150
*
150-
IF( XNORM.EQ.ZERO ) THEN
151+
IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
151152
*
152153
* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
153154
*

‎SRC/cunbdb5.f‎

Lines changed: 34 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -169,18 +169,21 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
169169
* =====================================================================
170170
*
171171
* .. Parameters ..
172+
REAL REALZERO
173+
PARAMETER ( REALZERO = 0.0E0 )
172174
COMPLEX ONE, ZERO
173175
PARAMETER ( ONE = (1.0E0,0.0E0), ZERO = (0.0E0,0.0E0) )
174176
* ..
175177
* .. Local Scalars ..
176178
INTEGER CHILDINFO, I, J
179+
REAL EPS, NORM, SCL, SSQ
177180
* ..
178181
* .. External Subroutines ..
179-
EXTERNAL CUNBDB6, XERBLA
182+
EXTERNAL CLASSQ, CUNBDB6, CSCAL, XERBLA
180183
* ..
181184
* .. External Functions ..
182-
REAL SCNRM2
183-
EXTERNAL SCNRM2
185+
REAL SLAMCH, SCNRM2
186+
EXTERNAL SLAMCH, SCNRM2
184187
* ..
185188
* .. Intrinsic Function ..
186189
INTRINSIC MAX
@@ -213,16 +216,33 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
213216
RETURN
214217
END IF
215218
*
216-
* Project X onto the orthogonal complement of Q
219+
EPS = SLAMCH( 'Precision' )
217220
*
218-
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2,
219-
$ WORK, LWORK, CHILDINFO )
221+
* Project X onto the orthogonal complement of Q if X is nonzero
220222
*
221-
* If the projection is nonzero, then return
223+
SCL = REALZERO
224+
SSQ = REALZERO
225+
CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
226+
CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
227+
NORM = SCL * SQRT( SSQ )
222228
*
223-
IF( SCNRM2(M1,X1,INCX1) .NE. ZERO
224-
$ .OR. SCNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
225-
RETURN
229+
IF( NORM .GT. N * EPS ) THEN
230+
* Scale vector to unit norm to avoid problems in the caller code.
231+
* Computing the reciprocal is undesirable but
232+
* * xLASCL cannot be used because of the vector increments and
233+
* * the round-off error has a negligible impact on
234+
* orthogonalization.
235+
CALL CSCAL( M1, ONE / NORM, X1, INCX1 )
236+
CALL CSCAL( M2, ONE / NORM, X2, INCX2 )
237+
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
238+
$ LDQ2, WORK, LWORK, CHILDINFO )
239+
*
240+
* If the projection is nonzero, then return
241+
*
242+
IF( SCNRM2(M1,X1,INCX1) .NE. REALZERO
243+
$ .OR. SCNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
244+
RETURN
245+
END IF
226246
END IF
227247
*
228248
* Project each standard basis vector e_1,...,e_M1 in turn, stopping
@@ -238,8 +258,8 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
238258
END DO
239259
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
240260
$ LDQ2, WORK, LWORK, CHILDINFO )
241-
IF( SCNRM2(M1,X1,INCX1) .NE. ZERO
242-
$ .OR. SCNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
261+
IF( SCNRM2(M1,X1,INCX1) .NE. REALZERO
262+
$ .OR. SCNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
243263
RETURN
244264
END IF
245265
END DO
@@ -257,8 +277,8 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
257277
X2(I) = ONE
258278
CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
259279
$ LDQ2, WORK, LWORK, CHILDINFO )
260-
IF( SCNRM2(M1,X1,INCX1) .NE. ZERO
261-
$ .OR. SCNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
280+
IF( SCNRM2(M1,X1,INCX1) .NE. REALZERO
281+
$ .OR. SCNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
262282
RETURN
263283
END IF
264284
END DO

‎SRC/cunbdb6.f‎

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -41,9 +41,8 @@
4141
*> with respect to the columns of
4242
*> Q = [ Q1 ] .
4343
*> [ Q2 ]
44-
*> The Euclidean norm of X must be one and the columns of Q must be
45-
*> orthonormal. The orthogonalized vector will be zero if and only if it
46-
*> lies entirely in the range of Q.
44+
*> The columns of Q must be orthonormal. The orthogonalized vector will
45+
*> be zero if and only if it lies entirely in the range of Q.
4746
*>
4847
*> The projection is computed with at most two iterations of the
4948
*> classical Gram-Schmidt algorithm, see
@@ -174,7 +173,7 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
174173
*
175174
* .. Parameters ..
176175
REAL ALPHA, REALONE, REALZERO
177-
PARAMETER ( ALPHA = 0.1E0, REALONE = 1.0E0,
176+
PARAMETER ( ALPHA = 0.83E0, REALONE = 1.0E0,
178177
$ REALZERO = 0.0E0 )
179178
COMPLEX NEGONE, ONE, ZERO
180179
PARAMETER ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0),
@@ -223,14 +222,16 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
223222
*
224223
EPS = SLAMCH( 'Precision' )
225224
*
225+
* Compute the Euclidean norm of X
226+
*
227+
SCL = REALZERO
228+
SSQ = REALZERO
229+
CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
230+
CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
231+
NORM = SCL * SQRT( SSQ )
232+
*
226233
* First, project X onto the orthogonal complement of Q's column
227234
* space
228-
*
229-
* Christoph Conrads: In debugging mode the norm should be computed
230-
* and an assertion added comparing the norm with one. Alas, Fortran
231-
* never made it into 1989 when assert() was introduced into the C
232-
* programming language.
233-
NORM = REALONE
234235
*
235236
IF( M1 .EQ. 0 ) THEN
236237
DO I = 1, N

‎SRC/dlarfgp.f‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
122122
* ..
123123
* .. Local Scalars ..
124124
INTEGER J, KNT
125-
DOUBLE PRECISION BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
125+
DOUBLE PRECISION BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM
126126
* ..
127127
* .. External Functions ..
128128
DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
@@ -141,11 +141,12 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
141141
RETURN
142142
END IF
143143
*
144+
EPS = DLAMCH( 'Precision' )
144145
XNORM = DNRM2( N-1, X, INCX )
145146
*
146-
IF( XNORM.EQ.ZERO ) THEN
147+
IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
147148
*
148-
* H = [+/-1, 0; I], sign chosen so ALPHA >= 0
149+
* H = [+/-1, 0; I], sign chosen so ALPHA >= 0.
149150
*
150151
IF( ALPHA.GE.ZERO ) THEN
151152
* When TAU.eq.ZERO, the vector is special-cased to be

‎SRC/dorbdb5.f‎

Lines changed: 34 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -169,18 +169,21 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
169169
* =====================================================================
170170
*
171171
* .. Parameters ..
172+
DOUBLE PRECISION REALZERO
173+
PARAMETER ( REALZERO = 0.0D0 )
172174
DOUBLE PRECISION ONE, ZERO
173175
PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
174176
* ..
175177
* .. Local Scalars ..
176178
INTEGER CHILDINFO, I, J
179+
DOUBLE PRECISION EPS, NORM, SCL, SSQ
177180
* ..
178181
* .. External Subroutines ..
179-
EXTERNAL DORBDB6, XERBLA
182+
EXTERNAL DLASSQ, DORBDB6, DSCAL, XERBLA
180183
* ..
181184
* .. External Functions ..
182-
DOUBLE PRECISION DNRM2
183-
EXTERNAL DNRM2
185+
DOUBLE PRECISION DLAMCH, DNRM2
186+
EXTERNAL DLAMCH, DNRM2
184187
* ..
185188
* .. Intrinsic Function ..
186189
INTRINSIC MAX
@@ -213,16 +216,33 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
213216
RETURN
214217
END IF
215218
*
216-
* Project X onto the orthogonal complement of Q
219+
EPS = DLAMCH( 'Precision' )
217220
*
218-
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2,
219-
$ WORK, LWORK, CHILDINFO )
221+
* Project X onto the orthogonal complement of Q if X is nonzero
220222
*
221-
* If the projection is nonzero, then return
223+
SCL = REALZERO
224+
SSQ = REALZERO
225+
CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
226+
CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
227+
NORM = SCL * SQRT( SSQ )
222228
*
223-
IF( DNRM2(M1,X1,INCX1) .NE. ZERO
224-
$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
225-
RETURN
229+
IF( NORM .GT. N * EPS ) THEN
230+
* Scale vector to unit norm to avoid problems in the caller code.
231+
* Computing the reciprocal is undesirable but
232+
* * xLASCL cannot be used because of the vector increments and
233+
* * the round-off error has a negligible impact on
234+
* orthogonalization.
235+
CALL DSCAL( M1, ONE / NORM, X1, INCX1 )
236+
CALL DSCAL( M2, ONE / NORM, X2, INCX2 )
237+
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
238+
$ LDQ2, WORK, LWORK, CHILDINFO )
239+
*
240+
* If the projection is nonzero, then return
241+
*
242+
IF( DNRM2(M1,X1,INCX1) .NE. REALZERO
243+
$ .OR. DNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
244+
RETURN
245+
END IF
226246
END IF
227247
*
228248
* Project each standard basis vector e_1,...,e_M1 in turn, stopping
@@ -238,8 +258,8 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
238258
END DO
239259
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
240260
$ LDQ2, WORK, LWORK, CHILDINFO )
241-
IF( DNRM2(M1,X1,INCX1) .NE. ZERO
242-
$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
261+
IF( DNRM2(M1,X1,INCX1) .NE. REALZERO
262+
$ .OR. DNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
243263
RETURN
244264
END IF
245265
END DO
@@ -257,8 +277,8 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
257277
X2(I) = ONE
258278
CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
259279
$ LDQ2, WORK, LWORK, CHILDINFO )
260-
IF( DNRM2(M1,X1,INCX1) .NE. ZERO
261-
$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
280+
IF( DNRM2(M1,X1,INCX1) .NE. REALZERO
281+
$ .OR. DNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN
262282
RETURN
263283
END IF
264284
END DO

‎SRC/dorbdb6.f‎

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -41,9 +41,8 @@
4141
*> with respect to the columns of
4242
*> Q = [ Q1 ] .
4343
*> [ Q2 ]
44-
*> The Euclidean norm of X must be one and the columns of Q must be
45-
*> orthonormal. The orthogonalized vector will be zero if and only if it
46-
*> lies entirely in the range of Q.
44+
*> The columns of Q must be orthonormal. The orthogonalized vector will
45+
*> be zero if and only if it lies entirely in the range of Q.
4746
*>
4847
*> The projection is computed with at most two iterations of the
4948
*> classical Gram-Schmidt algorithm, see
@@ -174,7 +173,7 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
174173
*
175174
* .. Parameters ..
176175
DOUBLE PRECISION ALPHA, REALONE, REALZERO
177-
PARAMETER ( ALPHA = 0.1D0, REALONE = 1.0D0,
176+
PARAMETER ( ALPHA = 0.83D0, REALONE = 1.0D0,
178177
$ REALZERO = 0.0D0 )
179178
DOUBLE PRECISION NEGONE, ONE, ZERO
180179
PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
@@ -222,14 +221,16 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
222221
*
223222
EPS = DLAMCH( 'Precision' )
224223
*
224+
* Compute the Euclidean norm of X
225+
*
226+
SCL = REALZERO
227+
SSQ = REALZERO
228+
CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
229+
CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
230+
NORM = SCL * SQRT( SSQ )
231+
*
225232
* First, project X onto the orthogonal complement of Q's column
226233
* space
227-
*
228-
* Christoph Conrads: In debugging mode the norm should be computed
229-
* and an assertion added comparing the norm with one. Alas, Fortran
230-
* never made it into 1989 when assert() was introduced into the C
231-
* programming language.
232-
NORM = REALONE
233234
*
234235
IF( M1 .EQ. 0 ) THEN
235236
DO I = 1, N

‎SRC/slarfgp.f‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
122122
* ..
123123
* .. Local Scalars ..
124124
INTEGER J, KNT
125-
REAL BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
125+
REAL BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM
126126
* ..
127127
* .. External Functions ..
128128
REAL SLAMCH, SLAPY2, SNRM2
@@ -141,9 +141,10 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
141141
RETURN
142142
END IF
143143
*
144+
EPS = SLAMCH( 'Precision' )
144145
XNORM = SNRM2( N-1, X, INCX )
145146
*
146-
IF( XNORM.EQ.ZERO ) THEN
147+
IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
147148
*
148149
* H = [+/-1, 0; I], sign chosen so ALPHA >= 0.
149150
*

0 commit comments

Comments
 (0)