[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.3 and 1.15

version 1.3, 2010/08/13 21:04:14 version 1.15, 2023/08/07 08:39:38
Line 1 Line 1
       SUBROUTINE ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )  *> \brief \b ZSYEQUB
 *  *
 *     -- LAPACK routine (version 3.2.2)                                 --  *  =========== DOCUMENTATION ===========
 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --  
 *     -- Jason Riedy of Univ. of California Berkeley.                 --  
 *     -- June 2010                                                    --  
 *  *
 *     -- LAPACK is a software package provided by Univ. of Tennessee, --  * Online html documentation available at
 *     -- Univ. of California Berkeley and NAG Ltd.                    --  *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download ZSYEQUB + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsyequb.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsyequb.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsyequb.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       INTEGER            INFO, LDA, N
   *       DOUBLE PRECISION   AMAX, SCOND
   *       CHARACTER          UPLO
   *       ..
   *       .. Array Arguments ..
   *       COMPLEX*16         A( LDA, * ), WORK( * )
   *       DOUBLE PRECISION   S( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> ZSYEQUB computes row and column scalings intended to equilibrate a
   *> symmetric matrix A (with respect to the Euclidean norm) and reduce
   *> its condition number. The scale factors S are computed by the BIN
   *> algorithm (see references) so that the scaled matrix B with elements
   *> B(i,j) = S(i)*A(i,j)*S(j) has a condition number within a factor N of
   *> the smallest possible condition number over all possible diagonal
   *> scalings.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] UPLO
   *> \verbatim
   *>          UPLO is CHARACTER*1
   *>          = 'U':  Upper triangle of A is stored;
   *>          = 'L':  Lower triangle of A is stored.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrix A. N >= 0.
   *> \endverbatim
   *>
   *> \param[in] A
   *> \verbatim
   *>          A is COMPLEX*16 array, dimension (LDA,N)
   *>          The N-by-N symmetric matrix whose scaling factors are to be
   *>          computed.
   *> \endverbatim
   *>
   *> \param[in] LDA
   *> \verbatim
   *>          LDA is INTEGER
   *>          The leading dimension of the array A. LDA >= max(1,N).
   *> \endverbatim
   *>
   *> \param[out] S
   *> \verbatim
   *>          S is DOUBLE PRECISION array, dimension (N)
   *>          If INFO = 0, S contains the scale factors for A.
   *> \endverbatim
   *>
   *> \param[out] SCOND
   *> \verbatim
   *>          SCOND is DOUBLE PRECISION
   *>          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
   *>          large nor too small, it is not worth scaling by S.
   *> \endverbatim
   *>
   *> \param[out] AMAX
   *> \verbatim
   *>          AMAX is DOUBLE PRECISION
   *>          Largest absolute value of any matrix element. If AMAX is
   *>          very close to overflow or very close to underflow, the
   *>          matrix should be scaled.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is COMPLEX*16 array, dimension (2*N)
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit
   *>          < 0:  if INFO = -i, the i-th argument had an illegal value
   *>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \ingroup complex16SYcomputational
   *
   *> \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 ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
   *
   *  -- LAPACK computational routine --
   *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *  *
       IMPLICIT NONE  
 *     ..  
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, N        INTEGER            INFO, LDA, N
       DOUBLE PRECISION   AMAX, SCOND        DOUBLE PRECISION   AMAX, SCOND
Line 20 Line 144
       DOUBLE PRECISION   S( * )        DOUBLE PRECISION   S( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  ZSYEQUB computes row and column scalings intended to equilibrate a  
 *  symmetric matrix A and reduce its condition number  
 *  (with respect to the two-norm).  S contains the scale factors,  
 *  S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with  
 *  elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This  
 *  choice of S puts the condition number of B within a factor N of the  
 *  smallest possible condition number over all possible diagonal  
 *  scalings.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  UPLO    (input) CHARACTER*1  
 *          Specifies whether the details of the factorization are stored  
 *          as an upper or lower triangular matrix.  
 *          = 'U':  Upper triangular, form is A = U*D*U**T;  
 *          = 'L':  Lower triangular, form is A = L*D*L**T.  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrix A.  N >= 0.  
 *  
 *  A       (input) COMPLEX*16 array, dimension (LDA,N)  
 *          The N-by-N symmetric matrix whose scaling  
 *          factors are to be computed.  Only the diagonal elements of A  
 *          are referenced.  
 *  
 *  LDA     (input) INTEGER  
 *          The leading dimension of the array A.  LDA >= max(1,N).  
 *  
 *  S       (output) DOUBLE PRECISION array, dimension (N)  
 *          If INFO = 0, S contains the scale factors for A.  
 *  
 *  SCOND   (output) DOUBLE PRECISION  
 *          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  
 *          large nor too small, it is not worth scaling by S.  
 *  
 *  AMAX    (output) DOUBLE PRECISION  
 *          Absolute value of largest matrix element.  If AMAX is very  
 *          close to overflow or very close to underflow, the matrix  
 *          should be scaled.  
 *  
 *  WORK    (workspace) COMPLEX*16 array, dimension (3*N)  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value  
 *          > 0:  if INFO = i, the i-th diagonal element is nonpositive.  
 *  
 *  Further Details  
 *  ======= =======  
 *  
 *  Reference: Livne, O.E. and Golub, G.H., "Scaling by Binormalization",  
 *  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004.  
 *  DOI 10.1023/B:NUMA.0000016606.32820.69  
 *  Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 101 Line 165
       EXTERNAL           DLAMCH, LSAME        EXTERNAL           DLAMCH, LSAME
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           ZLASSQ        EXTERNAL           ZLASSQ, XERBLA
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          ABS, DBLE, DIMAG, INT, LOG, MAX, MIN, SQRT        INTRINSIC          ABS, DBLE, DIMAG, INT, LOG, MAX, MIN, SQRT
Line 109 Line 173
 *     .. 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 118 Line 182
 *  *
       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 135 Line 199
 *     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 151 Line 215
                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 160 Line 224
             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 ) * DBLE( 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 ) * ( DBLE( WORK( I ) ) - T*SI )
               C0 = -(T*SI)*SI + 2 * DBLE( 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 + DBLE( WORK( I ) ) ) * D / N
               S( I ) = SI
            END DO
       END DO        END DO
   
  999  CONTINUE   999  CONTINUE
Line 268 Line 331
       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.3  
changed lines
  Added in v.1.15

CVSweb interface <joel.bertrand@systella.fr>