From 3d38da2bc401d56162a528547430b0b69664ea38 Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Sun, 12 Nov 2023 16:50:52 +0100 Subject: [PATCH] Make vector orthogonalization more reliable (Reference-LAPACK PR 930) --- lapack-netlib/SRC/clarfgp.f | 7 +++--- lapack-netlib/SRC/cunbdb5.f | 50 ++++++++++++++++++++++++++----------- lapack-netlib/SRC/cunbdb6.f | 21 ++++++++-------- lapack-netlib/SRC/dlarfgp.f | 9 ++++--- lapack-netlib/SRC/dorbdb5.f | 50 ++++++++++++++++++++++++++----------- lapack-netlib/SRC/dorbdb6.f | 21 ++++++++-------- lapack-netlib/SRC/slarfgp.f | 7 +++--- lapack-netlib/SRC/sorbdb5.f | 50 ++++++++++++++++++++++++++----------- lapack-netlib/SRC/sorbdb6.f | 21 ++++++++-------- lapack-netlib/SRC/zlarfgp.f | 7 +++--- lapack-netlib/SRC/zunbdb5.f | 50 ++++++++++++++++++++++++++----------- lapack-netlib/SRC/zunbdb6.f | 21 ++++++++-------- 12 files changed, 201 insertions(+), 113 deletions(-) diff --git a/lapack-netlib/SRC/clarfgp.f b/lapack-netlib/SRC/clarfgp.f index b584484c7f..47b5e47b07 100644 --- a/lapack-netlib/SRC/clarfgp.f +++ b/lapack-netlib/SRC/clarfgp.f @@ -97,7 +97,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexOTHERauxiliary +*> \ingroup larfgp * * ===================================================================== SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) @@ -122,7 +122,7 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) * .. * .. Local Scalars .. INTEGER J, KNT - REAL ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM + REAL ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM COMPLEX SAVEALPHA * .. * .. External Functions .. @@ -143,11 +143,12 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) RETURN END IF * + EPS = SLAMCH( 'Precision' ) XNORM = SCNRM2( N-1, X, INCX ) ALPHR = REAL( ALPHA ) ALPHI = AIMAG( ALPHA ) * - IF( XNORM.EQ.ZERO ) THEN + IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN * * H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0. * diff --git a/lapack-netlib/SRC/cunbdb5.f b/lapack-netlib/SRC/cunbdb5.f index d2ff4e7000..22513cf8b1 100644 --- a/lapack-netlib/SRC/cunbdb5.f +++ b/lapack-netlib/SRC/cunbdb5.f @@ -148,7 +148,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexOTHERcomputational +*> \ingroup unbdb5 * * ===================================================================== SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, @@ -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, CSCAL, XERBLA * .. * .. External Functions .. - REAL SCNRM2 - EXTERNAL SCNRM2 + REAL SLAMCH, SCNRM2 + EXTERNAL SLAMCH, SCNRM2 * .. * .. Intrinsic Function .. INTRINSIC MAX @@ -213,16 +216,33 @@ 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 +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL CSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL CSCAL( M2, ONE / NORM, X2, INCX2 ) + 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 @@ -238,8 +258,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 @@ -257,8 +277,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 diff --git a/lapack-netlib/SRC/cunbdb6.f b/lapack-netlib/SRC/cunbdb6.f index cd14d9295c..566fd76b7c 100644 --- a/lapack-netlib/SRC/cunbdb6.f +++ b/lapack-netlib/SRC/cunbdb6.f @@ -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 @@ -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), @@ -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 diff --git a/lapack-netlib/SRC/dlarfgp.f b/lapack-netlib/SRC/dlarfgp.f index 69845056d4..a8cf1b31e3 100644 --- a/lapack-netlib/SRC/dlarfgp.f +++ b/lapack-netlib/SRC/dlarfgp.f @@ -97,7 +97,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleOTHERauxiliary +*> \ingroup larfgp * * ===================================================================== SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) @@ -122,7 +122,7 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * .. * .. Local Scalars .. INTEGER J, KNT - DOUBLE PRECISION BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM + DOUBLE PRECISION BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 @@ -141,11 +141,12 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) RETURN END IF * + EPS = DLAMCH( 'Precision' ) XNORM = DNRM2( N-1, X, INCX ) * - IF( XNORM.EQ.ZERO ) THEN + IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN * -* H = [+/-1, 0; I], sign chosen so ALPHA >= 0 +* H = [+/-1, 0; I], sign chosen so ALPHA >= 0. * IF( ALPHA.GE.ZERO ) THEN * When TAU.eq.ZERO, the vector is special-cased to be diff --git a/lapack-netlib/SRC/dorbdb5.f b/lapack-netlib/SRC/dorbdb5.f index 6e057a05f8..cbd58ae547 100644 --- a/lapack-netlib/SRC/dorbdb5.f +++ b/lapack-netlib/SRC/dorbdb5.f @@ -148,7 +148,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleOTHERcomputational +*> \ingroup unbdb5 * * ===================================================================== SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, @@ -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, DSCAL, XERBLA * .. * .. External Functions .. - DOUBLE PRECISION DNRM2 - EXTERNAL DNRM2 + DOUBLE PRECISION DLAMCH, DNRM2 + EXTERNAL DLAMCH, DNRM2 * .. * .. Intrinsic Function .. INTRINSIC MAX @@ -213,16 +216,33 @@ 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 +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL DSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL DSCAL( M2, ONE / NORM, X2, INCX2 ) + 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 @@ -238,8 +258,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 @@ -257,8 +277,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 diff --git a/lapack-netlib/SRC/dorbdb6.f b/lapack-netlib/SRC/dorbdb6.f index 1428876847..3e356d0010 100644 --- a/lapack-netlib/SRC/dorbdb6.f +++ b/lapack-netlib/SRC/dorbdb6.f @@ -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 @@ -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 ) @@ -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 diff --git a/lapack-netlib/SRC/slarfgp.f b/lapack-netlib/SRC/slarfgp.f index df42980c4e..c28274c2c4 100644 --- a/lapack-netlib/SRC/slarfgp.f +++ b/lapack-netlib/SRC/slarfgp.f @@ -97,7 +97,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup realOTHERauxiliary +*> \ingroup larfgp * * ===================================================================== SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) @@ -122,7 +122,7 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * .. * .. Local Scalars .. INTEGER J, KNT - REAL BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM + REAL BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM * .. * .. External Functions .. REAL SLAMCH, SLAPY2, SNRM2 @@ -141,9 +141,10 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) RETURN END IF * + EPS = SLAMCH( 'Precision' ) XNORM = SNRM2( N-1, X, INCX ) * - IF( XNORM.EQ.ZERO ) THEN + IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN * * H = [+/-1, 0; I], sign chosen so ALPHA >= 0. * diff --git a/lapack-netlib/SRC/sorbdb5.f b/lapack-netlib/SRC/sorbdb5.f index 8c67aedfb1..8fb88876fe 100644 --- a/lapack-netlib/SRC/sorbdb5.f +++ b/lapack-netlib/SRC/sorbdb5.f @@ -148,7 +148,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup realOTHERcomputational +*> \ingroup unbdb5 * * ===================================================================== SUBROUTINE SORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, @@ -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, SSCAL, XERBLA * .. * .. External Functions .. - REAL SNRM2 - EXTERNAL SNRM2 + REAL SLAMCH, SNRM2 + EXTERNAL SLAMCH, SNRM2 * .. * .. Intrinsic Function .. INTRINSIC MAX @@ -213,16 +216,33 @@ 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 +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL SSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL SSCAL( M2, ONE / NORM, X2, INCX2 ) + 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 @@ -238,8 +258,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 @@ -257,8 +277,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 diff --git a/lapack-netlib/SRC/sorbdb6.f b/lapack-netlib/SRC/sorbdb6.f index d320c9e464..eac1777225 100644 --- a/lapack-netlib/SRC/sorbdb6.f +++ b/lapack-netlib/SRC/sorbdb6.f @@ -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 @@ -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 ) @@ -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 diff --git a/lapack-netlib/SRC/zlarfgp.f b/lapack-netlib/SRC/zlarfgp.f index 77eba8e869..6c9efb04c6 100644 --- a/lapack-netlib/SRC/zlarfgp.f +++ b/lapack-netlib/SRC/zlarfgp.f @@ -97,7 +97,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHERauxiliary +*> \ingroup larfgp * * ===================================================================== SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) @@ -122,7 +122,7 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) * .. * .. Local Scalars .. INTEGER J, KNT - DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM + DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM COMPLEX*16 SAVEALPHA * .. * .. External Functions .. @@ -143,11 +143,12 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) RETURN END IF * + EPS = DLAMCH( 'Precision' ) XNORM = DZNRM2( N-1, X, INCX ) ALPHR = DBLE( ALPHA ) ALPHI = DIMAG( ALPHA ) * - IF( XNORM.EQ.ZERO ) THEN + IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN * * H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0. * diff --git a/lapack-netlib/SRC/zunbdb5.f b/lapack-netlib/SRC/zunbdb5.f index 23174fe502..c451ae921a 100644 --- a/lapack-netlib/SRC/zunbdb5.f +++ b/lapack-netlib/SRC/zunbdb5.f @@ -148,7 +148,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHERcomputational +*> \ingroup unbdb5 * * ===================================================================== SUBROUTINE ZUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, @@ -169,18 +169,21 @@ SUBROUTINE ZUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * ===================================================================== * * .. Parameters .. + DOUBLE PRECISION REALZERO + PARAMETER ( REALZERO = 0.0D0 ) COMPLEX*16 ONE, ZERO PARAMETER ( ONE = (1.0D0,0.0D0), ZERO = (0.0D0,0.0D0) ) * .. * .. Local Scalars .. INTEGER CHILDINFO, I, J + DOUBLE PRECISION EPS, NORM, SCL, SSQ * .. * .. External Subroutines .. - EXTERNAL ZUNBDB6, XERBLA + EXTERNAL ZLASSQ, ZUNBDB6, ZSCAL, XERBLA * .. * .. External Functions .. - DOUBLE PRECISION DZNRM2 - EXTERNAL DZNRM2 + DOUBLE PRECISION DLAMCH, DZNRM2 + EXTERNAL DLAMCH, DZNRM2 * .. * .. Intrinsic Function .. INTRINSIC MAX @@ -213,16 +216,33 @@ SUBROUTINE ZUNBDB5( 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 ZUNBDB6( 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 ZLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL ZLASSQ( M2, X2, INCX2, SCL, SSQ ) + NORM = SCL * SQRT( SSQ ) * - IF( DZNRM2(M1,X1,INCX1) .NE. ZERO - $ .OR. DZNRM2(M2,X2,INCX2) .NE. ZERO ) THEN - RETURN + IF( NORM .GT. N * EPS ) THEN +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL ZSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL ZSCAL( M2, ONE / NORM, X2, INCX2 ) + CALL ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, + $ LDQ2, WORK, LWORK, CHILDINFO ) +* +* If the projection is nonzero, then return +* + IF( DZNRM2(M1,X1,INCX1) .NE. REALZERO + $ .OR. DZNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN + RETURN + END IF END IF * * Project each standard basis vector e_1,...,e_M1 in turn, stopping @@ -238,8 +258,8 @@ SUBROUTINE ZUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, END DO CALL ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, $ LDQ2, WORK, LWORK, CHILDINFO ) - IF( DZNRM2(M1,X1,INCX1) .NE. ZERO - $ .OR. DZNRM2(M2,X2,INCX2) .NE. ZERO ) THEN + IF( DZNRM2(M1,X1,INCX1) .NE. REALZERO + $ .OR. DZNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN RETURN END IF END DO @@ -257,8 +277,8 @@ SUBROUTINE ZUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, X2(I) = ONE CALL ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, $ LDQ2, WORK, LWORK, CHILDINFO ) - IF( DZNRM2(M1,X1,INCX1) .NE. ZERO - $ .OR. DZNRM2(M2,X2,INCX2) .NE. ZERO ) THEN + IF( DZNRM2(M1,X1,INCX1) .NE. REALZERO + $ .OR. DZNRM2(M2,X2,INCX2) .NE. REALZERO ) THEN RETURN END IF END DO diff --git a/lapack-netlib/SRC/zunbdb6.f b/lapack-netlib/SRC/zunbdb6.f index ac7fa4be3b..ddc9dfc61f 100644 --- a/lapack-netlib/SRC/zunbdb6.f +++ b/lapack-netlib/SRC/zunbdb6.f @@ -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 @@ -174,7 +173,7 @@ SUBROUTINE ZUNBDB6( 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 ) COMPLEX*16 NEGONE, ONE, ZERO PARAMETER ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0), @@ -223,14 +222,16 @@ SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, * EPS = DLAMCH( 'Precision' ) * +* Compute the Euclidean norm of X +* + SCL = REALZERO + SSQ = REALZERO + CALL ZLASSQ( M1, X1, INCX1, SCL, SSQ ) + CALL ZLASSQ( 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