001:       SUBROUTINE DPOEQUB( N, A, LDA, S, SCOND, AMAX, INFO )
002: *
003: *     -- LAPACK routine (version 3.2)                                 --
004: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
005: *     -- Jason Riedy of Univ. of California Berkeley.                 --
006: *     -- November 2008                                                --
007: *
008: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
009: *     -- Univ. of California Berkeley and NAG Ltd.                    --
010: *
011:       IMPLICIT NONE
012: *     ..
013: *     .. Scalar Arguments ..
014:       INTEGER            INFO, LDA, N
015:       DOUBLE PRECISION   AMAX, SCOND
016: *     ..
017: *     .. Array Arguments ..
018:       DOUBLE PRECISION   A( LDA, * ), S( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  DPOEQU computes row and column scalings intended to equilibrate a
025: *  symmetric positive definite matrix A and reduce its condition number
026: *  (with respect to the two-norm).  S contains the scale factors,
027: *  S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
028: *  elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
029: *  choice of S puts the condition number of B within a factor N of the
030: *  smallest possible condition number over all possible diagonal
031: *  scalings.
032: *
033: *  Arguments
034: *  =========
035: *
036: *  N       (input) INTEGER
037: *          The order of the matrix A.  N >= 0.
038: *
039: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
040: *          The N-by-N symmetric positive definite matrix whose scaling
041: *          factors are to be computed.  Only the diagonal elements of A
042: *          are referenced.
043: *
044: *  LDA     (input) INTEGER
045: *          The leading dimension of the array A.  LDA >= max(1,N).
046: *
047: *  S       (output) DOUBLE PRECISION array, dimension (N)
048: *          If INFO = 0, S contains the scale factors for A.
049: *
050: *  SCOND   (output) DOUBLE PRECISION
051: *          If INFO = 0, S contains the ratio of the smallest S(i) to
052: *          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
053: *          large nor too small, it is not worth scaling by S.
054: *
055: *  AMAX    (output) DOUBLE PRECISION
056: *          Absolute value of largest matrix element.  If AMAX is very
057: *          close to overflow or very close to underflow, the matrix
058: *          should be scaled.
059: *
060: *  INFO    (output) INTEGER
061: *          = 0:  successful exit
062: *          < 0:  if INFO = -i, the i-th argument had an illegal value
063: *          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
064: *
065: *  =====================================================================
066: *
067: *     .. Parameters ..
068:       DOUBLE PRECISION   ZERO, ONE
069:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
070: *     ..
071: *     .. Local Scalars ..
072:       INTEGER            I
073:       DOUBLE PRECISION   SMIN, BASE, TMP
074: *     ..
075: *     .. External Functions ..
076:       DOUBLE PRECISION   DLAMCH
077:       EXTERNAL           DLAMCH
078: *     ..
079: *     .. External Subroutines ..
080:       EXTERNAL           XERBLA
081: *     ..
082: *     .. Intrinsic Functions ..
083:       INTRINSIC          MAX, MIN, SQRT, LOG, INT
084: *     ..
085: *     .. Executable Statements ..
086: *
087: *     Test the input parameters.
088: *
089: *     Positive definite only performs 1 pass of equilibration.
090: *
091:       INFO = 0
092:       IF( N.LT.0 ) THEN
093:          INFO = -1
094:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
095:          INFO = -3
096:       END IF
097:       IF( INFO.NE.0 ) THEN
098:          CALL XERBLA( 'DPOEQUB', -INFO )
099:          RETURN
100:       END IF
101: *
102: *     Quick return if possible.
103: *
104:       IF( N.EQ.0 ) THEN
105:          SCOND = ONE
106:          AMAX = ZERO
107:          RETURN
108:       END IF
109: 
110:       BASE = DLAMCH( 'B' )
111:       TMP = -0.5D+0 / LOG ( BASE )
112: *
113: *     Find the minimum and maximum diagonal elements.
114: *
115:       S( 1 ) = A( 1, 1 )
116:       SMIN = S( 1 )
117:       AMAX = S( 1 )
118:       DO 10 I = 2, N
119:          S( I ) = A( I, I )
120:          SMIN = MIN( SMIN, S( I ) )
121:          AMAX = MAX( AMAX, S( I ) )
122:    10 CONTINUE
123: *
124:       IF( SMIN.LE.ZERO ) THEN
125: *
126: *        Find the first non-positive diagonal element and return.
127: *
128:          DO 20 I = 1, N
129:             IF( S( I ).LE.ZERO ) THEN
130:                INFO = I
131:                RETURN
132:             END IF
133:    20    CONTINUE
134:       ELSE
135: *
136: *        Set the scale factors to the reciprocals
137: *        of the diagonal elements.
138: *
139:          DO 30 I = 1, N
140:             S( I ) = BASE ** INT( TMP * LOG( S( I ) ) )
141:    30    CONTINUE
142: *
143: *        Compute SCOND = min(S(I)) / max(S(I)).
144: *
145:          SCOND = SQRT( SMIN ) / SQRT( AMAX )
146:       END IF
147: *
148:       RETURN
149: *
150: *     End of DPOEQUB
151: *
152:       END
153: