Annotation of rpl/lapack/lapack/dsyequb.f, revision 1.11

1.5       bertrand    1: *> \brief \b DSYEQUB
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.11    ! bertrand    5: * Online html documentation available at
        !             6: *            http://www.netlib.org/lapack/explore-html/
1.5       bertrand    7: *
                      8: *> \htmlonly
1.11    ! bertrand    9: *> Download DSYEQUB + dependencies
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyequb.f">
        !            11: *> [TGZ]</a>
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyequb.f">
        !            13: *> [ZIP]</a>
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyequb.f">
1.5       bertrand   15: *> [TXT]</a>
1.11    ! bertrand   16: *> \endhtmlonly
1.5       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
1.11    ! bertrand   22: *
1.5       bertrand   23: *       .. Scalar Arguments ..
                     24: *       INTEGER            INFO, LDA, N
                     25: *       DOUBLE PRECISION   AMAX, SCOND
                     26: *       CHARACTER          UPLO
                     27: *       ..
                     28: *       .. Array Arguments ..
                     29: *       DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( * )
                     30: *       ..
1.11    ! bertrand   31: *
1.5       bertrand   32: *
                     33: *> \par Purpose:
                     34: *  =============
                     35: *>
                     36: *> \verbatim
                     37: *>
                     38: *> DSYEQUB computes row and column scalings intended to equilibrate a
1.11    ! bertrand   39: *> symmetric matrix A (with respect to the Euclidean norm) and reduce
        !            40: *> its condition number. The scale factors S are computed by the BIN
        !            41: *> algorithm (see references) so that the scaled matrix B with elements
        !            42: *> B(i,j) = S(i)*A(i,j)*S(j) has a condition number within a factor N of
        !            43: *> the smallest possible condition number over all possible diagonal
1.5       bertrand   44: *> scalings.
                     45: *> \endverbatim
                     46: *
                     47: *  Arguments:
                     48: *  ==========
                     49: *
                     50: *> \param[in] UPLO
                     51: *> \verbatim
                     52: *>          UPLO is CHARACTER*1
1.11    ! bertrand   53: *>          = 'U':  Upper triangle of A is stored;
        !            54: *>          = 'L':  Lower triangle of A is stored.
1.5       bertrand   55: *> \endverbatim
                     56: *>
                     57: *> \param[in] N
                     58: *> \verbatim
                     59: *>          N is INTEGER
1.11    ! bertrand   60: *>          The order of the matrix A. N >= 0.
1.5       bertrand   61: *> \endverbatim
                     62: *>
                     63: *> \param[in] A
                     64: *> \verbatim
                     65: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
1.11    ! bertrand   66: *>          The N-by-N symmetric matrix whose scaling factors are to be
        !            67: *>          computed.
1.5       bertrand   68: *> \endverbatim
                     69: *>
                     70: *> \param[in] LDA
                     71: *> \verbatim
                     72: *>          LDA is INTEGER
1.11    ! bertrand   73: *>          The leading dimension of the array A. LDA >= max(1,N).
1.5       bertrand   74: *> \endverbatim
                     75: *>
                     76: *> \param[out] S
                     77: *> \verbatim
                     78: *>          S is DOUBLE PRECISION array, dimension (N)
                     79: *>          If INFO = 0, S contains the scale factors for A.
                     80: *> \endverbatim
                     81: *>
                     82: *> \param[out] SCOND
                     83: *> \verbatim
                     84: *>          SCOND is DOUBLE PRECISION
                     85: *>          If INFO = 0, S contains the ratio of the smallest S(i) to
1.11    ! bertrand   86: *>          the largest S(i). If SCOND >= 0.1 and AMAX is neither too
1.5       bertrand   87: *>          large nor too small, it is not worth scaling by S.
                     88: *> \endverbatim
                     89: *>
                     90: *> \param[out] AMAX
                     91: *> \verbatim
                     92: *>          AMAX is DOUBLE PRECISION
1.11    ! bertrand   93: *>          Largest absolute value of any matrix element. If AMAX is
        !            94: *>          very close to overflow or very close to underflow, the
        !            95: *>          matrix should be scaled.
1.5       bertrand   96: *> \endverbatim
                     97: *>
                     98: *> \param[out] WORK
                     99: *> \verbatim
1.11    ! bertrand  100: *>          WORK is DOUBLE PRECISION array, dimension (2*N)
1.5       bertrand  101: *> \endverbatim
                    102: *>
                    103: *> \param[out] INFO
                    104: *> \verbatim
                    105: *>          INFO is INTEGER
                    106: *>          = 0:  successful exit
                    107: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
                    108: *>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
                    109: *> \endverbatim
                    110: *
                    111: *  Authors:
                    112: *  ========
                    113: *
1.11    ! bertrand  114: *> \author Univ. of Tennessee
        !           115: *> \author Univ. of California Berkeley
        !           116: *> \author Univ. of Colorado Denver
        !           117: *> \author NAG Ltd.
1.5       bertrand  118: *
1.11    ! bertrand  119: *> \date December 2016
1.5       bertrand  120: *
                    121: *> \ingroup doubleSYcomputational
                    122: *
                    123: *> \par References:
                    124: *  ================
                    125: *>
                    126: *>  Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
                    127: *>  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
                    128: *>  DOI 10.1023/B:NUMA.0000016606.32820.69 \n
1.11    ! bertrand  129: *>  Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679
1.5       bertrand  130: *>
                    131: *  =====================================================================
1.1       bertrand  132:       SUBROUTINE DSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
                    133: *
1.11    ! bertrand  134: *  -- LAPACK computational routine (version 3.7.0) --
1.5       bertrand  135: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    136: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11    ! bertrand  137: *     December 2016
1.1       bertrand  138: *
                    139: *     .. Scalar Arguments ..
                    140:       INTEGER            INFO, LDA, N
                    141:       DOUBLE PRECISION   AMAX, SCOND
                    142:       CHARACTER          UPLO
                    143: *     ..
                    144: *     .. Array Arguments ..
                    145:       DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( * )
                    146: *     ..
                    147: *
                    148: *  =====================================================================
                    149: *
                    150: *     .. Parameters ..
                    151:       DOUBLE PRECISION   ONE, ZERO
1.11    ! bertrand  152:       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
1.1       bertrand  153:       INTEGER            MAX_ITER
                    154:       PARAMETER          ( MAX_ITER = 100 )
                    155: *     ..
                    156: *     .. Local Scalars ..
                    157:       INTEGER            I, J, ITER
                    158:       DOUBLE PRECISION   AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
                    159:      $                   SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
                    160:       LOGICAL            UP
                    161: *     ..
                    162: *     .. External Functions ..
                    163:       DOUBLE PRECISION   DLAMCH
                    164:       LOGICAL            LSAME
                    165:       EXTERNAL           DLAMCH, LSAME
                    166: *     ..
                    167: *     .. External Subroutines ..
                    168:       EXTERNAL           DLASSQ
                    169: *     ..
                    170: *     .. Intrinsic Functions ..
                    171:       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
                    172: *     ..
                    173: *     .. Executable Statements ..
                    174: *
1.11    ! bertrand  175: *     Test the input parameters.
1.1       bertrand  176: *
                    177:       INFO = 0
                    178:       IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN
1.11    ! bertrand  179:          INFO = -1
1.1       bertrand  180:       ELSE IF ( N .LT. 0 ) THEN
1.11    ! bertrand  181:          INFO = -2
1.1       bertrand  182:       ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN
1.11    ! bertrand  183:          INFO = -4
1.1       bertrand  184:       END IF
                    185:       IF ( INFO .NE. 0 ) THEN
1.11    ! bertrand  186:          CALL XERBLA( 'DSYEQUB', -INFO )
        !           187:          RETURN
1.1       bertrand  188:       END IF
                    189: 
                    190:       UP = LSAME( UPLO, 'U' )
                    191:       AMAX = ZERO
                    192: *
                    193: *     Quick return if possible.
                    194: *
                    195:       IF ( N .EQ. 0 ) THEN
1.11    ! bertrand  196:          SCOND = ONE
        !           197:          RETURN
1.1       bertrand  198:       END IF
                    199: 
                    200:       DO I = 1, N
1.11    ! bertrand  201:          S( I ) = ZERO
1.1       bertrand  202:       END DO
                    203: 
                    204:       AMAX = ZERO
                    205:       IF ( UP ) THEN
                    206:          DO J = 1, N
                    207:             DO I = 1, J-1
                    208:                S( I ) = MAX( S( I ), ABS( A( I, J ) ) )
                    209:                S( J ) = MAX( S( J ), ABS( A( I, J ) ) )
1.11    ! bertrand  210:                AMAX = MAX( AMAX, ABS( A( I, J ) ) )
1.1       bertrand  211:             END DO
                    212:             S( J ) = MAX( S( J ), ABS( A( J, J ) ) )
                    213:             AMAX = MAX( AMAX, ABS( A( J, J ) ) )
                    214:          END DO
                    215:       ELSE
                    216:          DO J = 1, N
                    217:             S( J ) = MAX( S( J ), ABS( A( J, J ) ) )
                    218:             AMAX = MAX( AMAX, ABS( A( J, J ) ) )
                    219:             DO I = J+1, N
                    220:                S( I ) = MAX( S( I ), ABS( A( I, J ) ) )
                    221:                S( J ) = MAX( S( J ), ABS( A( I, J ) ) )
                    222:                AMAX = MAX( AMAX, ABS( A( I, J ) ) )
                    223:             END DO
                    224:          END DO
                    225:       END IF
                    226:       DO J = 1, N
1.11    ! bertrand  227:          S( J ) = 1.0D0 / S( J )
1.1       bertrand  228:       END DO
                    229: 
1.11    ! bertrand  230:       TOL = ONE / SQRT( 2.0D0 * N )
1.1       bertrand  231: 
                    232:       DO ITER = 1, MAX_ITER
1.11    ! bertrand  233:          SCALE = 0.0D0
        !           234:          SUMSQ = 0.0D0
        !           235: *        beta = |A|s
        !           236:          DO I = 1, N
        !           237:             WORK( I ) = ZERO
        !           238:          END DO
        !           239:          IF ( UP ) THEN
        !           240:             DO J = 1, N
        !           241:                DO I = 1, J-1
        !           242:                   WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
        !           243:                   WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
        !           244:                END DO
        !           245:                WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
1.1       bertrand  246:             END DO
1.11    ! bertrand  247:          ELSE
        !           248:             DO J = 1, N
        !           249:                WORK( J ) = WORK( J ) + ABS( A( J, J ) ) * S( J )
        !           250:                DO I = J+1, N
        !           251:                   WORK( I ) = WORK( I ) + ABS( A( I, J ) ) * S( J )
        !           252:                   WORK( J ) = WORK( J ) + ABS( A( I, J ) ) * S( I )
        !           253:                END DO
1.1       bertrand  254:             END DO
1.11    ! bertrand  255:          END IF
        !           256: 
        !           257: *        avg = s^T beta / n
        !           258:          AVG = 0.0D0
        !           259:          DO I = 1, N
        !           260:             AVG = AVG + S( I )*WORK( I )
        !           261:          END DO
        !           262:          AVG = AVG / N
        !           263: 
        !           264:          STD = 0.0D0
        !           265:          DO I = N+1, 2*N
        !           266:             WORK( I ) = S( I-N ) * WORK( I-N ) - AVG
        !           267:          END DO
        !           268:          CALL DLASSQ( N, WORK( N+1 ), 1, SCALE, SUMSQ )
        !           269:          STD = SCALE * SQRT( SUMSQ / N )
1.1       bertrand  270: 
1.11    ! bertrand  271:          IF ( STD .LT. TOL * AVG ) GOTO 999
1.1       bertrand  272: 
1.11    ! bertrand  273:          DO I = 1, N
        !           274:             T = ABS( A( I, I ) )
        !           275:             SI = S( I )
        !           276:             C2 = ( N-1 ) * T
        !           277:             C1 = ( N-2 ) * ( WORK( I ) - T*SI )
        !           278:             C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG
        !           279:             D = C1*C1 - 4*C0*C2
        !           280: 
        !           281:             IF ( D .LE. 0 ) THEN
        !           282:                INFO = -1
        !           283:                RETURN
        !           284:             END IF
        !           285:             SI = -2*C0 / ( C1 + SQRT( D ) )
        !           286: 
        !           287:             D = SI - S( I )
        !           288:             U = ZERO
        !           289:             IF ( UP ) THEN
        !           290:                DO J = 1, I
        !           291:                   T = ABS( A( J, I ) )
        !           292:                   U = U + S( J )*T
        !           293:                   WORK( J ) = WORK( J ) + D*T
        !           294:                END DO
        !           295:                DO J = I+1,N
        !           296:                   T = ABS( A( I, J ) )
        !           297:                   U = U + S( J )*T
        !           298:                   WORK( J ) = WORK( J ) + D*T
        !           299:                END DO
        !           300:             ELSE
        !           301:                DO J = 1, I
        !           302:                   T = ABS( A( I, J ) )
        !           303:                   U = U + S( J )*T
        !           304:                   WORK( J ) = WORK( J ) + D*T
        !           305:                END DO
        !           306:                DO J = I+1,N
        !           307:                   T = ABS( A( J, I ) )
        !           308:                   U = U + S( J )*T
        !           309:                   WORK( J ) = WORK( J ) + D*T
        !           310:                END DO
        !           311:             END IF
1.1       bertrand  312: 
1.11    ! bertrand  313:             AVG = AVG + ( U + WORK( I ) ) * D / N
        !           314:             S( I ) = SI
        !           315:          END DO
1.1       bertrand  316:       END DO
                    317: 
                    318:  999  CONTINUE
                    319: 
                    320:       SMLNUM = DLAMCH( 'SAFEMIN' )
                    321:       BIGNUM = ONE / SMLNUM
                    322:       SMIN = BIGNUM
                    323:       SMAX = ZERO
1.11    ! bertrand  324:       T = ONE / SQRT( AVG )
1.1       bertrand  325:       BASE = DLAMCH( 'B' )
                    326:       U = ONE / LOG( BASE )
                    327:       DO I = 1, N
1.11    ! bertrand  328:          S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )
        !           329:          SMIN = MIN( SMIN, S( I ) )
        !           330:          SMAX = MAX( SMAX, S( I ) )
1.1       bertrand  331:       END DO
                    332:       SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
                    333: *
                    334:       END

CVSweb interface <joel.bertrand@systella.fr>