Skip to content

Commit 8b3bbfa

Browse files
langoulangou
authored
Merge pull request #647 from christoph-conrads/634-xORBDB6-performs-numerically-unadvisable-operations
Fix xORBDB6 performs numerically unadvisable operations
2 parents 3995eac + 145b1c3 commit 8b3bbfa

File tree

7 files changed

+192
-155
lines changed

7 files changed

+192
-155
lines changed

‎SRC/cunbdb6.f‎

Lines changed: 47 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -41,10 +41,16 @@
4141
*> with respect to the columns of
4242
*> Q = [ Q1 ] .
4343
*> [ Q2 ]
44-
*> The columns of Q must be orthonormal.
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.
4547
*>
46-
*> If the projection is zero according to Kahan's "twice is enough"
47-
*> criterion, then the zero vector is returned.
48+
*> The projection is computed with at most two iterations of the
49+
*> classical Gram-Schmidt algorithm, see
50+
*> * L. Giraud, J. Langou, M. Rozložník. "On the round-off error
51+
*> analysis of the Gram-Schmidt algorithm with reorthogonalization."
52+
*> 2002. CERFACS Technical Report No. TR/PA/02/33. URL:
53+
*> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf
4854
*>
4955
*>\endverbatim
5056
*
@@ -167,16 +173,19 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
167173
* =====================================================================
168174
*
169175
* .. Parameters ..
170-
REAL ALPHASQ, REALONE, REALZERO
171-
PARAMETER ( ALPHASQ = 0.01E0, REALONE = 1.0E0,
176+
REAL ALPHA, REALONE, REALZERO
177+
PARAMETER ( ALPHA = 0.01E0, REALONE = 1.0E0,
172178
$ REALZERO = 0.0E0 )
173179
COMPLEX NEGONE, ONE, ZERO
174180
PARAMETER ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0),
175181
$ ZERO = (0.0E0,0.0E0) )
176182
* ..
177183
* .. Local Scalars ..
178-
INTEGER I
179-
REAL NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
184+
INTEGER I, IX
185+
REAL EPS, NORM, NORM_NEW, SCL, SSQ
186+
* ..
187+
* .. External Functions ..
188+
REAL SLAMCH
180189
* ..
181190
* .. External Subroutines ..
182191
EXTERNAL CGEMV, CLASSQ, XERBLA
@@ -211,17 +220,17 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
211220
CALL XERBLA( 'CUNBDB6', -INFO )
212221
RETURN
213222
END IF
223+
*
224+
EPS = SLAMCH( 'Precision' )
214225
*
215226
* First, project X onto the orthogonal complement of Q's column
216227
* space
217228
*
218-
SCL1 = REALZERO
219-
SSQ1 = REALONE
220-
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
221-
SCL2 = REALZERO
222-
SSQ2 = REALONE
223-
CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
224-
NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
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
225234
*
226235
IF( M1 .EQ. 0 ) THEN
227236
DO I = 1, N
@@ -239,27 +248,31 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
239248
CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
240249
$ INCX2 )
241250
*
242-
SCL1 = REALZERO
243-
SSQ1 = REALONE
244-
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
245-
SCL2 = REALZERO
246-
SSQ2 = REALONE
247-
CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
248-
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
251+
SCL = REALZERO
252+
SSQ = REALZERO
253+
CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
254+
CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
255+
NORM_NEW = SCL * SQRT(SSQ)
249256
*
250257
* If projection is sufficiently large in norm, then stop.
251258
* If projection is zero, then stop.
252259
* Otherwise, project again.
253260
*
254-
IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN
261+
IF( NORM_NEW .GE. ALPHA * NORM ) THEN
255262
RETURN
256263
END IF
257264
*
258-
IF( NORMSQ2 .EQ. ZERO ) THEN
265+
IF( NORM_NEW .LE. N * EPS * NORM ) THEN
266+
DO IX = 1, 1 + (M1-1)*INCX1, INCX1
267+
X1( IX ) = ZERO
268+
END DO
269+
DO IX = 1, 1 + (M2-1)*INCX2, INCX2
270+
X2( IX ) = ZERO
271+
END DO
259272
RETURN
260273
END IF
261274
*
262-
NORMSQ1 = NORMSQ2
275+
NORM = NORM_NEW
263276
*
264277
DO I = 1, N
265278
WORK(I) = ZERO
@@ -281,24 +294,22 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
281294
CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
282295
$ INCX2 )
283296
*
284-
SCL1 = REALZERO
285-
SSQ1 = REALONE
286-
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
287-
SCL2 = REALZERO
288-
SSQ2 = REALONE
289-
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
290-
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
297+
SCL = REALZERO
298+
SSQ = REALZERO
299+
CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
300+
CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
301+
NORM_NEW = SCL * SQRT(SSQ)
291302
*
292303
* If second projection is sufficiently large in norm, then do
293304
* nothing more. Alternatively, if it shrunk significantly, then
294305
* truncate it to zero.
295306
*
296-
IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN
297-
DO I = 1, M1
298-
X1(I) = ZERO
307+
IF( NORM_NEW .LT. ALPHA * NORM ) THEN
308+
DO IX = 1, 1 + (M1-1)*INCX1, INCX1
309+
X1(IX) = ZERO
299310
END DO
300-
DO I = 1, M2
301-
X2(I) = ZERO
311+
DO IX = 1, 1 + (M2-1)*INCX2, INCX2
312+
X2(IX) = ZERO
302313
END DO
303314
END IF
304315
*
@@ -307,4 +318,3 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
307318
* End of CUNBDB6
308319
*
309320
END
310-

‎SRC/dorbdb6.f‎

Lines changed: 47 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -41,10 +41,16 @@
4141
*> with respect to the columns of
4242
*> Q = [ Q1 ] .
4343
*> [ Q2 ]
44-
*> The columns of Q must be orthonormal.
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.
4547
*>
46-
*> If the projection is zero according to Kahan's "twice is enough"
47-
*> criterion, then the zero vector is returned.
48+
*> The projection is computed with at most two iterations of the
49+
*> classical Gram-Schmidt algorithm, see
50+
*> * L. Giraud, J. Langou, M. Rozložník. "On the round-off error
51+
*> analysis of the Gram-Schmidt algorithm with reorthogonalization."
52+
*> 2002. CERFACS Technical Report No. TR/PA/02/33. URL:
53+
*> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf
4854
*>
4955
*>\endverbatim
5056
*
@@ -167,15 +173,18 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
167173
* =====================================================================
168174
*
169175
* .. Parameters ..
170-
DOUBLE PRECISION ALPHASQ, REALONE, REALZERO
171-
PARAMETER ( ALPHASQ = 0.01D0, REALONE = 1.0D0,
176+
DOUBLE PRECISION ALPHA, REALONE, REALZERO
177+
PARAMETER ( ALPHA = 0.01D0, REALONE = 1.0D0,
172178
$ REALZERO = 0.0D0 )
173179
DOUBLE PRECISION NEGONE, ONE, ZERO
174180
PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
175181
* ..
176182
* .. Local Scalars ..
177-
INTEGER I
178-
DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
183+
INTEGER I, IX
184+
REAL EPS, NORM, NORM_NEW, SCL, SSQ
185+
* ..
186+
* .. External Functions ..
187+
DOUBLE PRECISION DLAMCH
179188
* ..
180189
* .. External Subroutines ..
181190
EXTERNAL DGEMV, DLASSQ, XERBLA
@@ -210,17 +219,17 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
210219
CALL XERBLA( 'DORBDB6', -INFO )
211220
RETURN
212221
END IF
222+
*
223+
EPS = DLAMCH( 'Precision' )
213224
*
214225
* First, project X onto the orthogonal complement of Q's column
215226
* space
216227
*
217-
SCL1 = REALZERO
218-
SSQ1 = REALONE
219-
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
220-
SCL2 = REALZERO
221-
SSQ2 = REALONE
222-
CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
223-
NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
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
224233
*
225234
IF( M1 .EQ. 0 ) THEN
226235
DO I = 1, N
@@ -238,27 +247,31 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
238247
CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
239248
$ INCX2 )
240249
*
241-
SCL1 = REALZERO
242-
SSQ1 = REALONE
243-
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
244-
SCL2 = REALZERO
245-
SSQ2 = REALONE
246-
CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
247-
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
250+
SCL = REALZERO
251+
SSQ = REALZERO
252+
CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
253+
CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
254+
NORM_NEW = SCL * SQRT(SSQ)
248255
*
249256
* If projection is sufficiently large in norm, then stop.
250257
* If projection is zero, then stop.
251258
* Otherwise, project again.
252259
*
253-
IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN
260+
IF( NORM_NEW .GE. ALPHA * NORM ) THEN
254261
RETURN
255262
END IF
256263
*
257-
IF( NORMSQ2 .EQ. ZERO ) THEN
264+
IF( NORM_NEW .LE. N * EPS * NORM ) THEN
265+
DO IX = 1, 1 + (M1-1)*INCX1, INCX1
266+
X1( IX ) = ZERO
267+
END DO
268+
DO IX = 1, 1 + (M2-1)*INCX2, INCX2
269+
X2( IX ) = ZERO
270+
END DO
258271
RETURN
259272
END IF
260273
*
261-
NORMSQ1 = NORMSQ2
274+
NORM = NORM_NEW
262275
*
263276
DO I = 1, N
264277
WORK(I) = ZERO
@@ -280,24 +293,22 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
280293
CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
281294
$ INCX2 )
282295
*
283-
SCL1 = REALZERO
284-
SSQ1 = REALONE
285-
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
286-
SCL2 = REALZERO
287-
SSQ2 = REALONE
288-
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
289-
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
296+
SCL = REALZERO
297+
SSQ = REALZERO
298+
CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
299+
CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
300+
NORM_NEW = SCL * SQRT(SSQ)
290301
*
291302
* If second projection is sufficiently large in norm, then do
292303
* nothing more. Alternatively, if it shrunk significantly, then
293304
* truncate it to zero.
294305
*
295-
IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN
296-
DO I = 1, M1
297-
X1(I) = ZERO
306+
IF( NORM_NEW .LT. ALPHA * NORM ) THEN
307+
DO IX = 1, 1 + (M1-1)*INCX1, INCX1
308+
X1(IX) = ZERO
298309
END DO
299-
DO I = 1, M2
300-
X2(I) = ZERO
310+
DO IX = 1, 1 + (M2-1)*INCX2, INCX2
311+
X2(IX) = ZERO
301312
END DO
302313
END IF
303314
*
@@ -306,4 +317,3 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
306317
* End of DORBDB6
307318
*
308319
END
309-

‎SRC/sorbdb2.f‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -122,14 +122,14 @@
122122
*>
123123
*> \param[out] TAUP1
124124
*> \verbatim
125-
*> TAUP1 is REAL array, dimension (P)
125+
*> TAUP1 is REAL array, dimension (P-1)
126126
*> The scalar factors of the elementary reflectors that define
127127
*> P1.
128128
*> \endverbatim
129129
*>
130130
*> \param[out] TAUP2
131131
*> \verbatim
132-
*> TAUP2 is REAL array, dimension (M-P)
132+
*> TAUP2 is REAL array, dimension (Q)
133133
*> The scalar factors of the elementary reflectors that define
134134
*> P2.
135135
*> \endverbatim

‎SRC/sorbdb4.f‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -131,7 +131,7 @@
131131
*>
132132
*> \param[out] TAUP2
133133
*> \verbatim
134-
*> TAUP2 is REAL array, dimension (M-P)
134+
*> TAUP2 is REAL array, dimension (M-Q)
135135
*> The scalar factors of the elementary reflectors that define
136136
*> P2.
137137
*> \endverbatim

0 commit comments

Comments
 (0)