version 1.7, 2018/05/29 07:18:02
|
version 1.8, 2023/08/07 08:39:01
|
Line 41
|
Line 41
|
*> 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 |
*> If the projection is zero according to Kahan's "twice is enough" |
*> lies entirely in the range of Q. |
*> 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 |
*>\endverbatim |
* |
* |
Line 146
|
Line 152
|
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date July 2012 |
|
* |
|
*> \ingroup doubleOTHERcomputational |
*> \ingroup doubleOTHERcomputational |
* |
* |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, |
SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, |
$ LDQ2, WORK, LWORK, INFO ) |
$ LDQ2, WORK, LWORK, INFO ) |
* |
* |
* -- LAPACK computational routine (version 3.7.1) -- |
* -- LAPACK computational routine -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* July 2012 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2, |
INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2, |
Line 170
|
Line 173
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. 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 |
Line 214
|
Line 220
|
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 |
Line 241
|
Line 247
|
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 |
Line 283
|
Line 293
|
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 |
* |
* |
Line 309
|
Line 317
|
* End of DORBDB6 |
* End of DORBDB6 |
* |
* |
END |
END |
|
|