diff --git a/SRC/cunbdb6.f b/SRC/cunbdb6.f index 7acc99cb8b..b93a389d6b 100644 --- a/SRC/cunbdb6.f +++ b/SRC/cunbdb6.f @@ -41,10 +41,16 @@ *> with respect to the columns of *> Q = [ Q1 ] . *> [ Q2 ] -*> The columns of Q must be orthonormal. +*> 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. *> -*> If the projection is zero according to Kahan's "twice is enough" -*> 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 +*> analysis of the Gram-Schmidt algorithm with reorthogonalization." +*> 2002. CERFACS Technical Report No. TR/PA/02/33. URL: +*> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf *> *>\endverbatim * @@ -167,16 +173,19 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * ===================================================================== * * .. Parameters .. - REAL ALPHASQ, REALONE, REALZERO - PARAMETER ( ALPHASQ = 0.01E0, REALONE = 1.0E0, + REAL ALPHA, REALONE, REALZERO + PARAMETER ( ALPHA = 0.01E0, REALONE = 1.0E0, $ REALZERO = 0.0E0 ) COMPLEX NEGONE, ONE, ZERO PARAMETER ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0), $ ZERO = (0.0E0,0.0E0) ) * .. * .. Local Scalars .. - INTEGER I - REAL NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 + INTEGER I, IX + REAL EPS, NORM, NORM_NEW, SCL, SSQ +* .. +* .. External Functions .. + REAL SLAMCH * .. * .. External Subroutines .. EXTERNAL CGEMV, CLASSQ, XERBLA @@ -211,17 +220,17 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL XERBLA( 'CUNBDB6', -INFO ) RETURN END IF +* + EPS = SLAMCH( 'Precision' ) * * First, project X onto the orthogonal complement of Q's column * space * - SCL1 = REALZERO - SSQ1 = REALONE - CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2 +* 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 @@ -239,27 +248,31 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL CLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL CLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If projection is sufficiently large in norm, then stop. * If projection is zero, then stop. * Otherwise, project again. * - IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN + IF( NORM_NEW .GE. ALPHA * NORM ) THEN RETURN END IF * - IF( NORMSQ2 .EQ. ZERO ) THEN + IF( NORM_NEW .LE. N * EPS * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1( IX ) = ZERO + END DO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2( IX ) = ZERO + END DO RETURN END IF * - NORMSQ1 = NORMSQ2 + NORM = NORM_NEW * DO I = 1, N WORK(I) = ZERO @@ -281,24 +294,22 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL CLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL CLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If second projection is sufficiently large in norm, then do * nothing more. Alternatively, if it shrunk significantly, then * truncate it to zero. * - IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN - DO I = 1, M1 - X1(I) = ZERO + IF( NORM_NEW .LT. ALPHA * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1(IX) = ZERO END DO - DO I = 1, M2 - X2(I) = ZERO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2(IX) = ZERO END DO END IF * @@ -307,4 +318,3 @@ SUBROUTINE CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * End of CUNBDB6 * END - diff --git a/SRC/dorbdb6.f b/SRC/dorbdb6.f index fac52f760d..80bcc06e27 100644 --- a/SRC/dorbdb6.f +++ b/SRC/dorbdb6.f @@ -41,10 +41,16 @@ *> with respect to the columns of *> Q = [ Q1 ] . *> [ Q2 ] -*> The columns of Q must be orthonormal. +*> 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. *> -*> If the projection is zero according to Kahan's "twice is enough" -*> 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 +*> analysis of the Gram-Schmidt algorithm with reorthogonalization." +*> 2002. CERFACS Technical Report No. TR/PA/02/33. URL: +*> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf *> *>\endverbatim * @@ -167,15 +173,18 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ALPHASQ, REALONE, REALZERO - PARAMETER ( ALPHASQ = 0.01D0, REALONE = 1.0D0, + DOUBLE PRECISION ALPHA, REALONE, REALZERO + PARAMETER ( ALPHA = 0.01D0, REALONE = 1.0D0, $ REALZERO = 0.0D0 ) DOUBLE PRECISION NEGONE, ONE, ZERO PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 ) * .. * .. Local Scalars .. - INTEGER I - DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 + INTEGER I, IX + REAL EPS, NORM, NORM_NEW, SCL, SSQ +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH * .. * .. External Subroutines .. EXTERNAL DGEMV, DLASSQ, XERBLA @@ -210,17 +219,17 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL XERBLA( 'DORBDB6', -INFO ) RETURN END IF +* + EPS = DLAMCH( 'Precision' ) * * First, project X onto the orthogonal complement of Q's column * space * - SCL1 = REALZERO - SSQ1 = REALONE - CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2 +* 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 @@ -238,27 +247,31 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL DLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL DLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If projection is sufficiently large in norm, then stop. * If projection is zero, then stop. * Otherwise, project again. * - IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN + IF( NORM_NEW .GE. ALPHA * NORM ) THEN RETURN END IF * - IF( NORMSQ2 .EQ. ZERO ) THEN + IF( NORM_NEW .LE. N * EPS * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1( IX ) = ZERO + END DO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2( IX ) = ZERO + END DO RETURN END IF * - NORMSQ1 = NORMSQ2 + NORM = NORM_NEW * DO I = 1, N WORK(I) = ZERO @@ -280,24 +293,22 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL DLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL DLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If second projection is sufficiently large in norm, then do * nothing more. Alternatively, if it shrunk significantly, then * truncate it to zero. * - IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN - DO I = 1, M1 - X1(I) = ZERO + IF( NORM_NEW .LT. ALPHA * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1(IX) = ZERO END DO - DO I = 1, M2 - X2(I) = ZERO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2(IX) = ZERO END DO END IF * @@ -306,4 +317,3 @@ SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * End of DORBDB6 * END - diff --git a/SRC/sorbdb2.f b/SRC/sorbdb2.f index ad3eb269dd..484d352f8c 100644 --- a/SRC/sorbdb2.f +++ b/SRC/sorbdb2.f @@ -122,14 +122,14 @@ *> *> \param[out] TAUP1 *> \verbatim -*> TAUP1 is REAL array, dimension (P) +*> TAUP1 is REAL array, dimension (P-1) *> The scalar factors of the elementary reflectors that define *> P1. *> \endverbatim *> *> \param[out] TAUP2 *> \verbatim -*> TAUP2 is REAL array, dimension (M-P) +*> TAUP2 is REAL array, dimension (Q) *> The scalar factors of the elementary reflectors that define *> P2. *> \endverbatim diff --git a/SRC/sorbdb4.f b/SRC/sorbdb4.f index b18ed3b270..4d664435c1 100644 --- a/SRC/sorbdb4.f +++ b/SRC/sorbdb4.f @@ -131,7 +131,7 @@ *> *> \param[out] TAUP2 *> \verbatim -*> TAUP2 is REAL array, dimension (M-P) +*> TAUP2 is REAL array, dimension (M-Q) *> The scalar factors of the elementary reflectors that define *> P2. *> \endverbatim diff --git a/SRC/sorbdb6.f b/SRC/sorbdb6.f index a23b42bebc..b2449e3bed 100644 --- a/SRC/sorbdb6.f +++ b/SRC/sorbdb6.f @@ -41,10 +41,16 @@ *> with respect to the columns of *> Q = [ Q1 ] . *> [ Q2 ] -*> The columns of Q must be orthonormal. +*> 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. *> -*> If the projection is zero according to Kahan's "twice is enough" -*> 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 +*> analysis of the Gram-Schmidt algorithm with reorthogonalization." +*> 2002. CERFACS Technical Report No. TR/PA/02/33. URL: +*> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf *> *>\endverbatim * @@ -167,15 +173,18 @@ SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * ===================================================================== * * .. Parameters .. - REAL ALPHASQ, REALONE, REALZERO - PARAMETER ( ALPHASQ = 0.01E0, REALONE = 1.0E0, + REAL ALPHA, REALONE, REALZERO + PARAMETER ( ALPHA = 0.01E0, REALONE = 1.0E0, $ REALZERO = 0.0E0 ) REAL NEGONE, ONE, ZERO PARAMETER ( NEGONE = -1.0E0, ONE = 1.0E0, ZERO = 0.0E0 ) * .. * .. Local Scalars .. - INTEGER I - REAL NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 + INTEGER I, IX + REAL EPS, NORM, NORM_NEW, SCL, SSQ +* .. +* .. External Functions .. + REAL SLAMCH * .. * .. External Subroutines .. EXTERNAL SGEMV, SLASSQ, XERBLA @@ -210,17 +219,17 @@ SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL XERBLA( 'SORBDB6', -INFO ) RETURN END IF +* + EPS = SLAMCH( 'Precision' ) * * First, project X onto the orthogonal complement of Q's column * space * - SCL1 = REALZERO - SSQ1 = REALONE - CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL SLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2 +* 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 @@ -238,27 +247,31 @@ SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL SLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL SLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL SLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If projection is sufficiently large in norm, then stop. * If projection is zero, then stop. * Otherwise, project again. * - IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN + IF( NORM_NEW .GE. ALPHA * NORM ) THEN RETURN END IF * - IF( NORMSQ2 .EQ. ZERO ) THEN + IF( NORM_NEW .LE. N * EPS * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1( IX ) = ZERO + END DO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2( IX ) = ZERO + END DO RETURN END IF * - NORMSQ1 = NORMSQ2 + NORM = NORM_NEW * DO I = 1, N WORK(I) = ZERO @@ -280,24 +293,22 @@ SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL SLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL SLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If second projection is sufficiently large in norm, then do * nothing more. Alternatively, if it shrunk significantly, then * truncate it to zero. * - IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN - DO I = 1, M1 - X1(I) = ZERO + IF( NORM_NEW .LT. ALPHA * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1(IX) = ZERO END DO - DO I = 1, M2 - X2(I) = ZERO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2(IX) = ZERO END DO END IF * @@ -306,4 +317,3 @@ SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * End of SORBDB6 * END - diff --git a/SRC/sorcsd2by1.f b/SRC/sorcsd2by1.f index 25c317f6f6..a2d0163da1 100644 --- a/SRC/sorcsd2by1.f +++ b/SRC/sorcsd2by1.f @@ -672,13 +672,10 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * Accumulate Householder reflectors * IF( WANTU2 .AND. M-P .GT. 0 ) THEN - CALL SCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 ) + CALL SCOPY( M-Q, WORK(IORBDB+P), 1, U2, 1 ) END IF IF( WANTU1 .AND. P .GT. 0 ) THEN CALL SCOPY( P, WORK(IORBDB), 1, U1, 1 ) - DO J = 2, P - U1(1,J) = ZERO - END DO CALL SLACPY( 'L', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2), $ LDU1 ) CALL SORGQR( P, P, M-Q, U1, LDU1, WORK(ITAUP1), diff --git a/SRC/zunbdb6.f b/SRC/zunbdb6.f index ec681b597a..1fa7622ee3 100644 --- a/SRC/zunbdb6.f +++ b/SRC/zunbdb6.f @@ -41,10 +41,16 @@ *> with respect to the columns of *> Q = [ Q1 ] . *> [ Q2 ] -*> The columns of Q must be orthonormal. +*> 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. *> -*> If the projection is zero according to Kahan's "twice is enough" -*> 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 +*> analysis of the Gram-Schmidt algorithm with reorthogonalization." +*> 2002. CERFACS Technical Report No. TR/PA/02/33. URL: +*> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf *> *>\endverbatim * @@ -167,16 +173,19 @@ SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ALPHASQ, REALONE, REALZERO - PARAMETER ( ALPHASQ = 0.01D0, REALONE = 1.0D0, + DOUBLE PRECISION ALPHA, REALONE, REALZERO + PARAMETER ( ALPHA = 0.01D0, REALONE = 1.0D0, $ REALZERO = 0.0D0 ) COMPLEX*16 NEGONE, ONE, ZERO PARAMETER ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0), $ ZERO = (0.0D0,0.0D0) ) * .. * .. Local Scalars .. - INTEGER I - DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 + INTEGER I, IX + DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ +* .. +* .. External Functions .. + REAL DLAMCH * .. * .. External Subroutines .. EXTERNAL ZGEMV, ZLASSQ, XERBLA @@ -211,17 +220,17 @@ SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL XERBLA( 'ZUNBDB6', -INFO ) RETURN END IF +* + EPS = DLAMCH( 'Precision' ) * * First, project X onto the orthogonal complement of Q's column * space * - SCL1 = REALZERO - SSQ1 = REALONE - CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2 +* 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 @@ -239,27 +248,31 @@ SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL ZLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL ZLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If projection is sufficiently large in norm, then stop. * If projection is zero, then stop. * Otherwise, project again. * - IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN + IF( NORM_NEW .GE. ALPHA * NORM ) THEN RETURN END IF * - IF( NORMSQ2 .EQ. ZERO ) THEN + IF( NORM_NEW .LE. N * EPS * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1( IX ) = ZERO + END DO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2( IX ) = ZERO + END DO RETURN END IF * - NORMSQ1 = NORMSQ2 + NORM = NORM_NEW * DO I = 1, N WORK(I) = ZERO @@ -281,24 +294,22 @@ SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, $ INCX2 ) * - SCL1 = REALZERO - SSQ1 = REALONE - CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - SCL2 = REALZERO - SSQ2 = REALONE - CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) - NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2 + SCL = REALZERO + SSQ = REALZERO + CALL ZLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL ZLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM_NEW = SCL * SQRT(SSQ) * * If second projection is sufficiently large in norm, then do * nothing more. Alternatively, if it shrunk significantly, then * truncate it to zero. * - IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN - DO I = 1, M1 - X1(I) = ZERO + IF( NORM_NEW .LT. ALPHA * NORM ) THEN + DO IX = 1, 1 + (M1-1)*INCX1, INCX1 + X1(IX) = ZERO END DO - DO I = 1, M2 - X2(I) = ZERO + DO IX = 1, 1 + (M2-1)*INCX2, INCX2 + X2(IX) = ZERO END DO END IF * @@ -307,4 +318,3 @@ SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * End of ZUNBDB6 * END -