[BACK]Return to zunbdb6.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/zunbdb6.f between versions 1.7 and 1.8

version 1.7, 2018/05/29 07:18:41 version 1.8, 2023/08/07 08:39:43
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 complex16OTHERcomputational  *> \ingroup complex16OTHERcomputational
 *  *
 *  =====================================================================  *  =====================================================================
       SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,        SUBROUTINE ZUNBDB6( 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 )
       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
Line 215 Line 221
          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
Line 242 Line 248
       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
Line 284 Line 294
       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
 *  *
Line 310 Line 318
 *     End of ZUNBDB6  *     End of ZUNBDB6
 *  *
       END        END
   
Removed from v.1.7  
changed lines
  Added in v.1.8

CVSweb interface <joel.bertrand@systella.fr>