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

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

version 1.10, 2016/08/27 15:34:39 version 1.11, 2017/06/17 10:54:04
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 DSYEQUB + dependencies   *> Download DSYEQUB + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyequb.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyequb.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyequb.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyequb.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )  *       SUBROUTINE DSYEQUB( 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 28 Line 28
 *       .. Array Arguments ..  *       .. Array Arguments ..
 *       DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( * )  *       DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 36 Line 36
 *> \verbatim  *> \verbatim
 *>  *>
 *> DSYEQUB computes row and column scalings intended to equilibrate a  *> DSYEQUB 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 51 Line 50
 *> \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 DOUBLE PRECISION array, dimension (LDA,N)  *>          A is DOUBLE PRECISION 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 87 Line 83
 *> \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 DOUBLE PRECISION array, dimension (3*N)  *>          WORK is DOUBLE PRECISION array, dimension (2*N)
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] INFO  *> \param[out] INFO
Line 115 Line 111
 *  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 doubleSYcomputational  *> \ingroup doubleSYcomputational
 *  *
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 DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )        SUBROUTINE DSYEQUB( 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 153 Line 149
 *  *
 *     .. 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 )
 *     ..  *     ..
Line 176 Line 172
 *     ..  *     ..
 *     .. Executable Statements ..  *     .. 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( 'DSYEQUB', -INFO )           CALL XERBLA( 'DSYEQUB', -INFO )
         RETURN           RETURN
       END IF        END IF
   
       UP = LSAME( UPLO, 'U' )        UP = LSAME( UPLO, 'U' )
Line 197 Line 193
 *     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 211 Line 207
             DO I = 1, J-1              DO I = 1, J-1
                S( I ) = MAX( S( I ), ABS( A( I, J ) ) )                 S( I ) = MAX( S( I ), ABS( A( I, J ) ) )
                S( J ) = MAX( S( J ), ABS( A( I, J ) ) )                 S( J ) = MAX( S( J ), ABS( A( I, J ) ) )
                AMAX = MAX( AMAX, ABS( A(I, J) ) )                 AMAX = MAX( AMAX, ABS( A( I, J ) ) )
             END DO              END DO
             S( J ) = MAX( S( J ), ABS( A( J, J ) ) )              S( J ) = MAX( S( J ), ABS( A( J, J ) ) )
             AMAX = MAX( AMAX, ABS( A( J, J ) ) )              AMAX = MAX( AMAX, ABS( A( J, J ) ) )
Line 228 Line 224
          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 = ABS( A( I, J ) )                    WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
                  WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )                    WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
                  WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )                 END DO
               END DO                 WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
               WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )  
            END DO  
         ELSE  
            DO J = 1, N  
               WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )  
               DO I = J+1, N  
                  T = ABS( A( I, J ) )  
                  WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )  
                  WORK( J ) = WORK( J ) + ABS( 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 DLASSQ( 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 = ABS( 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 = ABS( A( J, I ) )  
               U = U + S( J )*T  
               WORK( J ) = WORK( J ) + D*T  
             END DO  
             DO J = I+1,N  
               T = ABS( A( I, J ) )  
               U = U + S( J )*T  
               WORK( J ) = WORK( J ) + D*T  
             END DO  
           ELSE  
             DO J = 1, I  
               T = ABS( 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 = ABS( A( J, I ) )              DO J = 1, N
               U = U + S( J )*T                 WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
               WORK( J ) = WORK( J ) + D*T                 DO I = J+1, N
                     WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
                     WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
                  END DO
             END DO              END DO
           END IF           END IF
   
           AVG = AVG + ( U + WORK( I ) ) * D / N  *        avg = s^T beta / n
           S( I ) = SI           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 DLASSQ( N, WORK( N+1 ), 1, SCALE, SUMSQ )
            STD = SCALE * SQRT( SUMSQ / N )
   
         END DO           IF ( STD .LT. TOL * AVG ) GOTO 999
   
            DO I = 1, N
               T = ABS( 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 = ABS( A( J, I ) )
                     U = U + S( J )*T
                     WORK( J ) = WORK( J ) + D*T
                  END DO
                  DO J = I+1,N
                     T = ABS( A( I, J ) )
                     U = U + S( J )*T
                     WORK( J ) = WORK( J ) + D*T
                  END DO
               ELSE
                  DO J = 1, I
                     T = ABS( A( I, J ) )
                     U = U + S( J )*T
                     WORK( J ) = WORK( J ) + D*T
                  END DO
                  DO J = I+1,N
                     T = ABS( 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 329 Line 321
       BIGNUM = ONE / SMLNUM        BIGNUM = ONE / SMLNUM
       SMIN = BIGNUM        SMIN = BIGNUM
       SMAX = ZERO        SMAX = ZERO
       T = ONE / SQRT(AVG)        T = ONE / SQRT( AVG )
       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>