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