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

Diff for /rpl/lapack/lapack/zsyequb.f between versions 1.10 and 1.11

version 1.10, 2016/08/27 15:35:07 version 1.11, 2017/06/17 10:54:28
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 ZSYEQUB + dependencies   *> Download ZSYEQUB + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsyequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsyequb.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsyequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsyequb.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsyequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsyequb.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )  *       SUBROUTINE ZSYEQUB( 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
 *>  *>
 *> ZSYEQUB computes row and column scalings intended to equilibrate a  *> ZSYEQUB computes row and column scalings intended to equilibrate a
 *> symmetric matrix A and reduce its condition number  *> symmetric 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
 *>          Specifies whether the details of the factorization are stored  *>          = 'U':  Upper triangle of A is stored;
 *>          as an upper or lower triangular matrix.  *>          = 'L':  Lower triangle of A is stored.
 *>          = 'U':  Upper triangular, form is A = U*D*U**T;  
 *>          = 'L':  Lower triangular, form is A = L*D*L**T.  
 *> \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 symmetric matrix whose scaling  *>          The N-by-N symmetric 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 88 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 116 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 November 2011  *> \date December 2016
 *  *
 *> \ingroup complex16SYcomputational  *> \ingroup complex16SYcomputational
 *  *
Line 130 Line 126
 *>  *>
 *>  Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n  *>  Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
 *>  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n  *>  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
 *>  DOI 10.1023/B:NUMA.0000016606.32820.69 \n   *>  DOI 10.1023/B:NUMA.0000016606.32820.69 \n
 *>  Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf  *>  Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
 *>  *>
 *  =====================================================================  *  =====================================================================
       SUBROUTINE ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )        SUBROUTINE ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.4.0) --  *  -- 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..--
 *     November 2011  *     December 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, N        INTEGER            INFO, LDA, N
Line 180 Line 176
 *     .. Statement Functions ..  *     .. Statement Functions ..
       DOUBLE PRECISION   CABS1        DOUBLE PRECISION   CABS1
 *     ..  *     ..
 *     Statement Function Definitions  *     .. Statement Function Definitions ..
       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )        CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 *     ..  *     ..
 *     .. Executable Statements ..  *     .. Executable Statements ..
Line 189 Line 185
 *  *
       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( 'ZSYEQUB', -INFO )           CALL XERBLA( 'ZSYEQUB', -INFO )
         RETURN           RETURN
       END IF        END IF
   
       UP = LSAME( UPLO, 'U' )        UP = LSAME( UPLO, 'U' )
Line 206 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 222 Line 218
                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
             S( J ) = MAX( S( J ), CABS1( A( J, J) ) )              S( J ) = MAX( S( J ), CABS1( A( J, J ) ) )
             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )              AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
          END DO           END DO
       ELSE        ELSE
Line 231 Line 227
             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )              AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
             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 = N+1, 2*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              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  *        avg = s^T beta / n
         END DO           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, 2*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 339 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 )
 *  *
Removed from v.1.10  
changed lines
  Added in v.1.11

CVSweb interface <joel.bertrand@systella.fr>