floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Make vector orthogonalization more reliable (Reference-LAPACK PR 930)

This commit is contained in:
Martin Kroeker 2023-11-12 16:50:52 +01:00 committed by GitHub
parent d58c88cf42
commit 3d38da2bc4
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
12 changed files with 201 additions and 113 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
lapack-netlib/SRC/clarfgp.f
View File

@ -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 @@
* ..
* .. 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 @@
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.
*

50
lapack-netlib/SRC/cunbdb5.f
View File

@ -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 @@
* =====================================================================
*
* .. 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 @@
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 @@
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 @@
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

21
lapack-netlib/SRC/cunbdb6.f
View File

@ -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 @@
*
* .. 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 @@
*
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

9
lapack-netlib/SRC/dlarfgp.f
View File

@ -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 @@
* ..
* .. 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 @@
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

50
lapack-netlib/SRC/dorbdb5.f
View File

@ -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 @@
* =====================================================================
*
* .. 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 @@
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 @@
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 @@
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

21
lapack-netlib/SRC/dorbdb6.f
View File

@ -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 @@
*
* .. 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 @@
*
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

7
lapack-netlib/SRC/slarfgp.f
View File

@ -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 @@
* ..
* .. 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 @@
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.
*

50
lapack-netlib/SRC/sorbdb5.f
View File

@ -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 @@
* =====================================================================
*
* .. 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 @@
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 @@
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 @@
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

21
lapack-netlib/SRC/sorbdb6.f
View File

@ -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 @@
*
* .. 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 @@
*
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

7
lapack-netlib/SRC/zlarfgp.f
View File

@ -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 @@
* ..
* .. 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 @@
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.
*

50
lapack-netlib/SRC/zunbdb5.f
View File

@ -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 @@
* =====================================================================
*
* .. 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 @@
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 @@
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 @@
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

21
lapack-netlib/SRC/zunbdb6.f
View File

@ -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 @@
*
* .. 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 @@
*
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