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

Merge pull request #3831 from martin-frbg/lapack647+697+702 

Fix code and documentation for ?SORBDB?/?CUNBDB? (Reference-LAPACK PRs 647+697+702)
This commit is contained in:
Martin Kroeker 2022-11-20 19:34:41 +01:00 committed by GitHub
parent b1102fe250 aaea0804bc
commit f63c93274c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
12 changed files with 204 additions and 164 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
lapack-netlib/SRC/cunbdb2.f
View File

@ -122,14 +122,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is COMPLEX array, dimension (P) *> TAUP1 is COMPLEX array, dimension (P-1)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is COMPLEX array, dimension (M-P) *> TAUP2 is COMPLEX array, dimension (Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

4
lapack-netlib/SRC/cunbdb4.f
View File

@ -124,14 +124,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is COMPLEX array, dimension (P) *> TAUP1 is COMPLEX array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is COMPLEX array, dimension (M-P) *> TAUP2 is COMPLEX array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

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

@ -41,10 +41,16 @@
*> with respect to the columns of *> with respect to the columns of
*> Q = [ Q1 ] . *> Q = [ Q1 ] .
*> [ Q2 ] *> [ 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" *> The projection is computed with at most two iterations of the
*> criterion, then the zero vector is returned. *> 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 *>\endverbatim
* *
@ -167,16 +173,19 @@
* ===================================================================== * =====================================================================
* *
* .. Parameters .. * .. Parameters ..
REAL ALPHASQ, REALONE, REALZERO REAL ALPHA, REALONE, REALZERO
PARAMETER ( ALPHASQ = 0.01E0, REALONE = 1.0E0, PARAMETER ( ALPHA = 0.01E0, REALONE = 1.0E0,
$ REALZERO = 0.0E0 ) $ REALZERO = 0.0E0 )
COMPLEX NEGONE, ONE, ZERO COMPLEX NEGONE, ONE, ZERO
PARAMETER ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0), PARAMETER ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0),
$ ZERO = (0.0E0,0.0E0) ) $ ZERO = (0.0E0,0.0E0) )
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I INTEGER I, IX
REAL NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 REAL EPS, NORM, NORM_NEW, SCL, SSQ
* ..
* .. External Functions ..
REAL SLAMCH
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CGEMV, CLASSQ, XERBLA EXTERNAL CGEMV, CLASSQ, XERBLA
@ -211,17 +220,17 @@
CALL XERBLA( 'CUNBDB6', -INFO ) CALL XERBLA( 'CUNBDB6', -INFO )
RETURN RETURN
END IF END IF
*
EPS = SLAMCH( 'Precision' )
* *
* First, project X onto the orthogonal complement of Q's column * First, project X onto the orthogonal complement of Q's column
* space * space
* *
SCL1 = REALZERO * Christoph Conrads: In debugging mode the norm should be computed
SSQ1 = REALONE * and an assertion added comparing the norm with one. Alas, Fortran
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) * never made it into 1989 when assert() was introduced into the C
SCL2 = REALZERO * programming language.
SSQ2 = REALONE NORM = REALONE
CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
IF( M1 .EQ. 0 ) THEN IF( M1 .EQ. 0 ) THEN
DO I = 1, N DO I = 1, N
@ -239,27 +248,31 @@
CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL CLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If projection is sufficiently large in norm, then stop. * If projection is sufficiently large in norm, then stop.
* If projection is zero, then stop. * If projection is zero, then stop.
* Otherwise, project again. * Otherwise, project again.
* *
IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .GE. ALPHA * NORM ) THEN
RETURN RETURN
END IF 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 RETURN
END IF END IF
* *
NORMSQ1 = NORMSQ2 NORM = NORM_NEW
* *
DO I = 1, N DO I = 1, N
WORK(I) = ZERO WORK(I) = ZERO
@ -281,24 +294,22 @@
CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL CGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL CLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL CLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL CLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If second projection is sufficiently large in norm, then do * If second projection is sufficiently large in norm, then do
* nothing more. Alternatively, if it shrunk significantly, then * nothing more. Alternatively, if it shrunk significantly, then
* truncate it to zero. * truncate it to zero.
* *
IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .LT. ALPHA * NORM ) THEN
DO I = 1, M1 DO IX = 1, 1 + (M1-1)*INCX1, INCX1
X1(I) = ZERO X1(IX) = ZERO
END DO END DO
DO I = 1, M2 DO IX = 1, 1 + (M2-1)*INCX2, INCX2
X2(I) = ZERO X2(IX) = ZERO
END DO END DO
END IF END IF
* *
@ -307,4 +318,3 @@
* End of CUNBDB6 * End of CUNBDB6
* *
END END

4
lapack-netlib/SRC/dorbdb2.f
View File

@ -122,14 +122,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is DOUBLE PRECISION array, dimension (P) *> TAUP1 is DOUBLE PRECISION array, dimension (P-1)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is DOUBLE PRECISION array, dimension (M-P) *> TAUP2 is DOUBLE PRECISION array, dimension (Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

4
lapack-netlib/SRC/dorbdb4.f
View File

@ -124,14 +124,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is DOUBLE PRECISION array, dimension (P) *> TAUP1 is DOUBLE PRECISION array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is DOUBLE PRECISION array, dimension (M-P) *> TAUP2 is DOUBLE PRECISION array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

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

@ -41,10 +41,16 @@
*> with respect to the columns of *> with respect to the columns of
*> Q = [ Q1 ] . *> Q = [ Q1 ] .
*> [ Q2 ] *> [ 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" *> The projection is computed with at most two iterations of the
*> criterion, then the zero vector is returned. *> 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 *>\endverbatim
* *
@ -167,15 +173,18 @@
* ===================================================================== * =====================================================================
* *
* .. Parameters .. * .. Parameters ..
DOUBLE PRECISION ALPHASQ, REALONE, REALZERO DOUBLE PRECISION ALPHA, REALONE, REALZERO
PARAMETER ( ALPHASQ = 0.01D0, REALONE = 1.0D0, PARAMETER ( ALPHA = 0.01D0, REALONE = 1.0D0,
$ REALZERO = 0.0D0 ) $ REALZERO = 0.0D0 )
DOUBLE PRECISION NEGONE, ONE, ZERO DOUBLE PRECISION NEGONE, ONE, ZERO
PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 ) PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I INTEGER I, IX
DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL DGEMV, DLASSQ, XERBLA EXTERNAL DGEMV, DLASSQ, XERBLA
@ -210,17 +219,17 @@
CALL XERBLA( 'DORBDB6', -INFO ) CALL XERBLA( 'DORBDB6', -INFO )
RETURN RETURN
END IF END IF
*
EPS = DLAMCH( 'Precision' )
* *
* First, project X onto the orthogonal complement of Q's column * First, project X onto the orthogonal complement of Q's column
* space * space
* *
SCL1 = REALZERO * Christoph Conrads: In debugging mode the norm should be computed
SSQ1 = REALONE * and an assertion added comparing the norm with one. Alas, Fortran
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) * never made it into 1989 when assert() was introduced into the C
SCL2 = REALZERO * programming language.
SSQ2 = REALONE NORM = REALONE
CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
IF( M1 .EQ. 0 ) THEN IF( M1 .EQ. 0 ) THEN
DO I = 1, N DO I = 1, N
@ -238,27 +247,31 @@
CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If projection is sufficiently large in norm, then stop. * If projection is sufficiently large in norm, then stop.
* If projection is zero, then stop. * If projection is zero, then stop.
* Otherwise, project again. * Otherwise, project again.
* *
IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .GE. ALPHA * NORM ) THEN
RETURN RETURN
END IF 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 RETURN
END IF END IF
* *
NORMSQ1 = NORMSQ2 NORM = NORM_NEW
* *
DO I = 1, N DO I = 1, N
WORK(I) = ZERO WORK(I) = ZERO
@ -280,24 +293,22 @@
CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL DLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If second projection is sufficiently large in norm, then do * If second projection is sufficiently large in norm, then do
* nothing more. Alternatively, if it shrunk significantly, then * nothing more. Alternatively, if it shrunk significantly, then
* truncate it to zero. * truncate it to zero.
* *
IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .LT. ALPHA * NORM ) THEN
DO I = 1, M1 DO IX = 1, 1 + (M1-1)*INCX1, INCX1
X1(I) = ZERO X1(IX) = ZERO
END DO END DO
DO I = 1, M2 DO IX = 1, 1 + (M2-1)*INCX2, INCX2
X2(I) = ZERO X2(IX) = ZERO
END DO END DO
END IF END IF
* *
@ -306,4 +317,3 @@
* End of DORBDB6 * End of DORBDB6
* *
END END

4
lapack-netlib/SRC/sorbdb2.f
View File

@ -122,14 +122,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is REAL array, dimension (P) *> TAUP1 is REAL array, dimension (P-1)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is REAL array, dimension (M-P) *> TAUP2 is REAL array, dimension (Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

4
lapack-netlib/SRC/sorbdb4.f
View File

@ -124,14 +124,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is REAL array, dimension (P) *> TAUP1 is REAL array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is REAL array, dimension (M-P) *> TAUP2 is REAL array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

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

@ -41,10 +41,16 @@
*> with respect to the columns of *> with respect to the columns of
*> Q = [ Q1 ] . *> Q = [ Q1 ] .
*> [ Q2 ] *> [ 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" *> The projection is computed with at most two iterations of the
*> criterion, then the zero vector is returned. *> 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 *>\endverbatim
* *
@ -167,15 +173,18 @@
* ===================================================================== * =====================================================================
* *
* .. Parameters .. * .. Parameters ..
REAL ALPHASQ, REALONE, REALZERO REAL ALPHA, REALONE, REALZERO
PARAMETER ( ALPHASQ = 0.01E0, REALONE = 1.0E0, PARAMETER ( ALPHA = 0.01E0, REALONE = 1.0E0,
$ REALZERO = 0.0E0 ) $ REALZERO = 0.0E0 )
REAL NEGONE, ONE, ZERO REAL NEGONE, ONE, ZERO
PARAMETER ( NEGONE = -1.0E0, ONE = 1.0E0, ZERO = 0.0E0 ) PARAMETER ( NEGONE = -1.0E0, ONE = 1.0E0, ZERO = 0.0E0 )
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I INTEGER I, IX
REAL NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 REAL EPS, NORM, NORM_NEW, SCL, SSQ
* ..
* .. External Functions ..
REAL SLAMCH
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL SGEMV, SLASSQ, XERBLA EXTERNAL SGEMV, SLASSQ, XERBLA
@ -210,17 +219,17 @@
CALL XERBLA( 'SORBDB6', -INFO ) CALL XERBLA( 'SORBDB6', -INFO )
RETURN RETURN
END IF END IF
*
EPS = SLAMCH( 'Precision' )
* *
* First, project X onto the orthogonal complement of Q's column * First, project X onto the orthogonal complement of Q's column
* space * space
* *
SCL1 = REALZERO * Christoph Conrads: In debugging mode the norm should be computed
SSQ1 = REALONE * and an assertion added comparing the norm with one. Alas, Fortran
CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) * never made it into 1989 when assert() was introduced into the C
SCL2 = REALZERO * programming language.
SSQ2 = REALONE NORM = REALONE
CALL SLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
IF( M1 .EQ. 0 ) THEN IF( M1 .EQ. 0 ) THEN
DO I = 1, N DO I = 1, N
@ -238,27 +247,31 @@
CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL SLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL SLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL SLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If projection is sufficiently large in norm, then stop. * If projection is sufficiently large in norm, then stop.
* If projection is zero, then stop. * If projection is zero, then stop.
* Otherwise, project again. * Otherwise, project again.
* *
IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .GE. ALPHA * NORM ) THEN
RETURN RETURN
END IF 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 RETURN
END IF END IF
* *
NORMSQ1 = NORMSQ2 NORM = NORM_NEW
* *
DO I = 1, N DO I = 1, N
WORK(I) = ZERO WORK(I) = ZERO
@ -280,24 +293,22 @@
CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL SLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL SLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If second projection is sufficiently large in norm, then do * If second projection is sufficiently large in norm, then do
* nothing more. Alternatively, if it shrunk significantly, then * nothing more. Alternatively, if it shrunk significantly, then
* truncate it to zero. * truncate it to zero.
* *
IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .LT. ALPHA * NORM ) THEN
DO I = 1, M1 DO IX = 1, 1 + (M1-1)*INCX1, INCX1
X1(I) = ZERO X1(IX) = ZERO
END DO END DO
DO I = 1, M2 DO IX = 1, 1 + (M2-1)*INCX2, INCX2
X2(I) = ZERO X2(IX) = ZERO
END DO END DO
END IF END IF
* *
@ -306,4 +317,3 @@
* End of SORBDB6 * End of SORBDB6
* *
END END

4
lapack-netlib/SRC/zunbdb2.f
View File

@ -122,14 +122,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is COMPLEX*16 array, dimension (P) *> TAUP1 is COMPLEX*16 array, dimension (P-1)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is COMPLEX*16 array, dimension (M-P) *> TAUP2 is COMPLEX*16 array, dimension (Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

4
lapack-netlib/SRC/zunbdb4.f
View File

@ -124,14 +124,14 @@
*> *>
*> \param[out] TAUP1 *> \param[out] TAUP1
*> \verbatim *> \verbatim
*> TAUP1 is COMPLEX*16 array, dimension (P) *> TAUP1 is COMPLEX*16 array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P1. *> P1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] TAUP2 *> \param[out] TAUP2
*> \verbatim *> \verbatim
*> TAUP2 is COMPLEX*16 array, dimension (M-P) *> TAUP2 is COMPLEX*16 array, dimension (M-Q)
*> The scalar factors of the elementary reflectors that define *> The scalar factors of the elementary reflectors that define
*> P2. *> P2.
*> \endverbatim *> \endverbatim

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

@ -41,10 +41,16 @@
*> with respect to the columns of *> with respect to the columns of
*> Q = [ Q1 ] . *> Q = [ Q1 ] .
*> [ Q2 ] *> [ 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" *> The projection is computed with at most two iterations of the
*> criterion, then the zero vector is returned. *> 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 *>\endverbatim
* *
@ -167,16 +173,19 @@
* ===================================================================== * =====================================================================
* *
* .. Parameters .. * .. Parameters ..
DOUBLE PRECISION ALPHASQ, REALONE, REALZERO DOUBLE PRECISION ALPHA, REALONE, REALZERO
PARAMETER ( ALPHASQ = 0.01D0, REALONE = 1.0D0, PARAMETER ( ALPHA = 0.01D0, REALONE = 1.0D0,
$ REALZERO = 0.0D0 ) $ REALZERO = 0.0D0 )
COMPLEX*16 NEGONE, ONE, ZERO COMPLEX*16 NEGONE, ONE, ZERO
PARAMETER ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0), PARAMETER ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0),
$ ZERO = (0.0D0,0.0D0) ) $ ZERO = (0.0D0,0.0D0) )
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I INTEGER I, IX
DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2 DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZGEMV, ZLASSQ, XERBLA EXTERNAL ZGEMV, ZLASSQ, XERBLA
@ -211,17 +220,17 @@
CALL XERBLA( 'ZUNBDB6', -INFO ) CALL XERBLA( 'ZUNBDB6', -INFO )
RETURN RETURN
END IF END IF
*
EPS = DLAMCH( 'Precision' )
* *
* First, project X onto the orthogonal complement of Q's column * First, project X onto the orthogonal complement of Q's column
* space * space
* *
SCL1 = REALZERO * Christoph Conrads: In debugging mode the norm should be computed
SSQ1 = REALONE * and an assertion added comparing the norm with one. Alas, Fortran
CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) * never made it into 1989 when assert() was introduced into the C
SCL2 = REALZERO * programming language.
SSQ2 = REALONE NORM = REALONE
CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
IF( M1 .EQ. 0 ) THEN IF( M1 .EQ. 0 ) THEN
DO I = 1, N DO I = 1, N
@ -239,27 +248,31 @@
CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL ZLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL ZLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If projection is sufficiently large in norm, then stop. * If projection is sufficiently large in norm, then stop.
* If projection is zero, then stop. * If projection is zero, then stop.
* Otherwise, project again. * Otherwise, project again.
* *
IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .GE. ALPHA * NORM ) THEN
RETURN RETURN
END IF 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 RETURN
END IF END IF
* *
NORMSQ1 = NORMSQ2 NORM = NORM_NEW
* *
DO I = 1, N DO I = 1, N
WORK(I) = ZERO WORK(I) = ZERO
@ -281,24 +294,22 @@
CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2, CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
$ INCX2 ) $ INCX2 )
* *
SCL1 = REALZERO SCL = REALZERO
SSQ1 = REALONE SSQ = REALZERO
CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 ) CALL ZLASSQ( M1, X1, INCX1, SCL, SSQ )
SCL2 = REALZERO CALL ZLASSQ( M2, X2, INCX2, SCL, SSQ )
SSQ2 = REALONE NORM_NEW = SCL * SQRT(SSQ)
CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
* *
* If second projection is sufficiently large in norm, then do * If second projection is sufficiently large in norm, then do
* nothing more. Alternatively, if it shrunk significantly, then * nothing more. Alternatively, if it shrunk significantly, then
* truncate it to zero. * truncate it to zero.
* *
IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN IF( NORM_NEW .LT. ALPHA * NORM ) THEN
DO I = 1, M1 DO IX = 1, 1 + (M1-1)*INCX1, INCX1
X1(I) = ZERO X1(IX) = ZERO
END DO END DO
DO I = 1, M2 DO IX = 1, 1 + (M2-1)*INCX2, INCX2
X2(I) = ZERO X2(IX) = ZERO
END DO END DO
END IF END IF
* *
@ -307,4 +318,3 @@
* End of ZUNBDB6 * End of ZUNBDB6
* *
END END