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

Diff for /rpl/lapack/lapack/zheequb.f between versions 1.11 and 1.12

version 1.11, 2016/08/27 15:34:49 version 1.12, 2017/06/17 10:54:14
Line 2 Line 2
 *  *
 *  =========== DOCUMENTATION ===========  *  =========== DOCUMENTATION ===========
 *  *
 * Online html documentation available at   * Online html documentation available at
 *            http://www.netlib.org/lapack/explore-html/   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *> \htmlonly  *> \htmlonly
 *> Download ZHEEQUB + dependencies   *> Download ZHEEQUB + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheequb.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheequb.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheequb.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )  *       SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       INTEGER            INFO, LDA, N  *       INTEGER            INFO, LDA, N
 *       DOUBLE PRECISION   AMAX, SCOND  *       DOUBLE PRECISION   AMAX, SCOND
Line 29 Line 29
 *       COMPLEX*16         A( LDA, * ), WORK( * )  *       COMPLEX*16         A( LDA, * ), WORK( * )
 *       DOUBLE PRECISION   S( * )  *       DOUBLE PRECISION   S( * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 37 Line 37
 *> \verbatim  *> \verbatim
 *>  *>
 *> ZHEEQUB computes row and column scalings intended to equilibrate a  *> ZHEEQUB computes row and column scalings intended to equilibrate a
 *> Hermitian matrix A and reduce its condition number  *> Hermitian matrix A (with respect to the Euclidean norm) and reduce
 *> (with respect to the two-norm).  S contains the scale factors,  *> its condition number. The scale factors S are computed by the BIN
 *> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with  *> algorithm (see references) so that the scaled matrix B with elements
 *> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This  *> B(i,j) = S(i)*A(i,j)*S(j) has a condition number within a factor N of
 *> choice of S puts the condition number of B within a factor N of the  *> the smallest possible condition number over all possible diagonal
 *> smallest possible condition number over all possible diagonal  
 *> scalings.  *> scalings.
 *> \endverbatim  *> \endverbatim
 *  *
Line 52 Line 51
 *> \param[in] UPLO  *> \param[in] UPLO
 *> \verbatim  *> \verbatim
 *>          UPLO is CHARACTER*1  *>          UPLO is CHARACTER*1
 *>          = 'U':  Upper triangles of A and B are stored;  *>          = 'U':  Upper triangle of A is stored;
 *>          = 'L':  Lower triangles of A and B are stored.  *>          = 'L':  Lower triangle of A is stored.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] N  *> \param[in] N
 *> \verbatim  *> \verbatim
 *>          N is INTEGER  *>          N is INTEGER
 *>          The order of the matrix A.  N >= 0.  *>          The order of the matrix A. N >= 0.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] A  *> \param[in] A
 *> \verbatim  *> \verbatim
 *>          A is COMPLEX*16 array, dimension (LDA,N)  *>          A is COMPLEX*16 array, dimension (LDA,N)
 *>          The N-by-N Hermitian matrix whose scaling  *>          The N-by-N Hermitian matrix whose scaling factors are to be
 *>          factors are to be computed.  Only the diagonal elements of A  *>          computed.
 *>          are referenced.  
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LDA  *> \param[in] LDA
 *> \verbatim  *> \verbatim
 *>          LDA is INTEGER  *>          LDA is INTEGER
 *>          The leading dimension of the array A.  LDA >= max(1,N).  *>          The leading dimension of the array A. LDA >= max(1,N).
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] S  *> \param[out] S
Line 86 Line 84
 *> \verbatim  *> \verbatim
 *>          SCOND is DOUBLE PRECISION  *>          SCOND is DOUBLE PRECISION
 *>          If INFO = 0, S contains the ratio of the smallest S(i) to  *>          If INFO = 0, S contains the ratio of the smallest S(i) to
 *>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too  *>          the largest S(i). If SCOND >= 0.1 and AMAX is neither too
 *>          large nor too small, it is not worth scaling by S.  *>          large nor too small, it is not worth scaling by S.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] AMAX  *> \param[out] AMAX
 *> \verbatim  *> \verbatim
 *>          AMAX is DOUBLE PRECISION  *>          AMAX is DOUBLE PRECISION
 *>          Absolute value of largest matrix element.  If AMAX is very  *>          Largest absolute value of any matrix element. If AMAX is
 *>          close to overflow or very close to underflow, the matrix  *>          very close to overflow or very close to underflow, the
 *>          should be scaled.  *>          matrix should be scaled.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] WORK  *> \param[out] WORK
 *> \verbatim  *> \verbatim
 *>          WORK is COMPLEX*16 array, dimension (3*N)  *>          WORK is COMPLEX*16 array, dimension (2*N)
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] INFO  *> \param[out] INFO
Line 114 Line 112
 *  Authors:  *  Authors:
 *  ========  *  ========
 *  *
 *> \author Univ. of Tennessee   *> \author Univ. of Tennessee
 *> \author Univ. of California Berkeley   *> \author Univ. of California Berkeley
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.   *> \author NAG Ltd.
 *  *
 *> \date April 2012  *> \date April 2012
 *  *
 *> \ingroup complex16HEcomputational  *> \ingroup complex16HEcomputational
 *  *
   *> \par References:
   *  ================
   *>
   *>  Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
   *>  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
   *>  DOI 10.1023/B:NUMA.0000016606.32820.69 \n
   *>  Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
   *>
 *  =====================================================================  *  =====================================================================
       SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )        SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.4.1) --  *  -- LAPACK computational routine (version 3.7.0) --
 *  -- 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..--
 *     April 2012  *     April 2012
Line 145 Line 151
 *  *
 *     .. Parameters ..  *     .. Parameters ..
       DOUBLE PRECISION   ONE, ZERO        DOUBLE PRECISION   ONE, ZERO
       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )        PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
       INTEGER            MAX_ITER        INTEGER            MAX_ITER
       PARAMETER          ( MAX_ITER = 100 )        PARAMETER          ( MAX_ITER = 100 )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       INTEGER            I, J, ITER        INTEGER            I, J, ITER
       DOUBLE PRECISION   AVG, STD, TOL, C0, C1, C2, T, U, SI, D,        DOUBLE PRECISION   AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
      $                   BASE, SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ       $                   SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
       LOGICAL            UP        LOGICAL            UP
       COMPLEX*16         ZDUM        COMPLEX*16         ZDUM
 *     ..  *     ..
Line 172 Line 178
 *     ..  *     ..
 *     .. Statement Function Definitions ..  *     .. Statement Function Definitions ..
       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )        CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
   *     ..
   *     .. Executable Statements ..
 *  *
 *     Test input parameters.  *     Test the input parameters.
 *  *
       INFO = 0        INFO = 0
       IF (.NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN        IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN
         INFO = -1           INFO = -1
       ELSE IF ( N .LT. 0 ) THEN        ELSE IF ( N .LT. 0 ) THEN
         INFO = -2           INFO = -2
       ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN        ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN
         INFO = -4           INFO = -4
       END IF        END IF
       IF ( INFO .NE. 0 ) THEN        IF ( INFO .NE. 0 ) THEN
         CALL XERBLA( 'ZHEEQUB', -INFO )           CALL XERBLA( 'ZHEEQUB', -INFO )
         RETURN           RETURN
       END IF        END IF
   
       UP = LSAME( UPLO, 'U' )        UP = LSAME( UPLO, 'U' )
Line 194 Line 202
 *     Quick return if possible.  *     Quick return if possible.
 *  *
       IF ( N .EQ. 0 ) THEN        IF ( N .EQ. 0 ) THEN
         SCOND = ONE           SCOND = ONE
         RETURN           RETURN
       END IF        END IF
   
       DO I = 1, N        DO I = 1, N
         S( I ) = ZERO           S( I ) = ZERO
       END DO        END DO
   
       AMAX = ZERO        AMAX = ZERO
Line 220 Line 228
             DO I = J+1, N              DO I = J+1, N
                S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )                 S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )
                S( J ) = MAX( S( J ), CABS1( A( I, J ) ) )                 S( J ) = MAX( S( J ), CABS1( A( I, J ) ) )
                AMAX = MAX( AMAX, CABS1( A(I, J ) ) )                 AMAX = MAX( AMAX, CABS1( A( I, J ) ) )
             END DO              END DO
          END DO           END DO
       END IF        END IF
       DO J = 1, N        DO J = 1, N
          S( J ) = 1.0D+0 / S( J )           S( J ) = 1.0D0 / S( J )
       END DO        END DO
   
       TOL = ONE / SQRT( 2.0D0 * N )        TOL = ONE / SQRT( 2.0D0 * N )
   
       DO ITER = 1, MAX_ITER        DO ITER = 1, MAX_ITER
          SCALE = 0.0D+0           SCALE = 0.0D0
          SUMSQ = 0.0D+0           SUMSQ = 0.0D0
 *       beta = |A|s  *        beta = |A|s
         DO I = 1, N           DO I = 1, N
            WORK( I ) = ZERO              WORK( I ) = ZERO
         END DO           END DO
         IF ( UP ) THEN           IF ( UP ) THEN
            DO J = 1, N              DO J = 1, N
               DO I = 1, J-1                 DO I = 1, J-1
                  T = CABS1( A( I, J ) )                    WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )                    WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )                 END DO
               END DO                 WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )  
            END DO  
         ELSE  
            DO J = 1, N  
               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )  
               DO I = J+1, N  
                  T = CABS1( A( I, J ) )  
                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )  
                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )  
               END DO  
            END DO  
         END IF  
   
 *       avg = s^T beta / n  
         AVG = 0.0D+0  
         DO I = 1, N  
           AVG = AVG + S( I )*WORK( I )  
         END DO  
         AVG = AVG / N  
   
         STD = 0.0D+0  
         DO I = 2*N+1, 3*N  
            WORK( I ) = S( I-2*N ) * WORK( I-2*N ) - AVG  
         END DO  
         CALL ZLASSQ( N, WORK( 2*N+1 ), 1, SCALE, SUMSQ )  
         STD = SCALE * SQRT( SUMSQ / N )  
   
         IF ( STD .LT. TOL * AVG ) GOTO 999  
   
         DO I = 1, N  
           T = CABS1( A( I, I ) )  
           SI = S( I )  
           C2 = ( N-1 ) * T  
           C1 = ( N-2 ) * ( WORK( I ) - T*SI )  
           C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG  
   
           D = C1*C1 - 4*C0*C2  
           IF ( D .LE. 0 ) THEN  
             INFO = -1  
             RETURN  
           END IF  
           SI = -2*C0 / ( C1 + SQRT( D ) )  
   
           D = SI - S(I)  
           U = ZERO  
           IF ( UP ) THEN  
             DO J = 1, I  
               T = CABS1( A( J, I ) )  
               U = U + S( J )*T  
               WORK( J ) = WORK( J ) + D*T  
             END DO  
             DO J = I+1,N  
               T = CABS1( A( I, J ) )  
               U = U + S( J )*T  
               WORK( J ) = WORK( J ) + D*T  
             END DO  
           ELSE  
             DO J = 1, I  
               T = CABS1( A( I, J ) )  
               U = U + S( J )*T  
               WORK( J ) = WORK( J ) + D*T  
             END DO              END DO
             DO J = I+1,N           ELSE
               T = CABS1( A( J, I ) )              DO J = 1, N
               U = U + S( J )*T                 WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
               WORK( J ) = WORK( J ) + D*T                 DO I = J+1, N
                     WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
                     WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
                  END DO
             END DO              END DO
           END IF           END IF
           AVG = AVG + ( U + WORK( I ) ) * D / N  
           S( I ) = SI  
         END DO  
   
   *        avg = s^T beta / n
            AVG = 0.0D0
            DO I = 1, N
               AVG = AVG + S( I )*WORK( I )
            END DO
            AVG = AVG / N
   
            STD = 0.0D0
            DO I = N+1, N
               WORK( I ) = S( I-N ) * WORK( I-N ) - AVG
            END DO
            CALL ZLASSQ( N, WORK( N+1 ), 1, SCALE, SUMSQ )
            STD = SCALE * SQRT( SUMSQ / N )
   
            IF ( STD .LT. TOL * AVG ) GOTO 999
   
            DO I = 1, N
               T = CABS1( A( I, I ) )
               SI = S( I )
               C2 = ( N-1 ) * T
               C1 = ( N-2 ) * ( WORK( I ) - T*SI )
               C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG
               D = C1*C1 - 4*C0*C2
   
               IF ( D .LE. 0 ) THEN
                  INFO = -1
                  RETURN
               END IF
               SI = -2*C0 / ( C1 + SQRT( D ) )
   
               D = SI - S( I )
               U = ZERO
               IF ( UP ) THEN
                  DO J = 1, I
                     T = CABS1( A( J, I ) )
                     U = U + S( J )*T
                     WORK( J ) = WORK( J ) + D*T
                  END DO
                  DO J = I+1,N
                     T = CABS1( A( I, J ) )
                     U = U + S( J )*T
                     WORK( J ) = WORK( J ) + D*T
                  END DO
               ELSE
                  DO J = 1, I
                     T = CABS1( A( I, J ) )
                     U = U + S( J )*T
                     WORK( J ) = WORK( J ) + D*T
                  END DO
                  DO J = I+1,N
                     T = CABS1( A( J, I ) )
                     U = U + S( J )*T
                     WORK( J ) = WORK( J ) + D*T
                  END DO
               END IF
   
               AVG = AVG + ( U + WORK( I ) ) * D / N
               S( I ) = SI
            END DO
       END DO        END DO
   
  999  CONTINUE   999  CONTINUE
Line 328 Line 334
       BASE = DLAMCH( 'B' )        BASE = DLAMCH( 'B' )
       U = ONE / LOG( BASE )        U = ONE / LOG( BASE )
       DO I = 1, N        DO I = 1, N
         S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )           S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )
         SMIN = MIN( SMIN, S( I ) )           SMIN = MIN( SMIN, S( I ) )
         SMAX = MAX( SMAX, S( I ) )           SMAX = MAX( SMAX, S( I ) )
       END DO        END DO
       SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )        SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
   *
       END        END
Removed from v.1.11  
changed lines
  Added in v.1.12

CVSweb interface <joel.bertrand@systella.fr>