001:       SUBROUTINE CPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          UPLO
009:       INTEGER            INFO, N
010:       REAL               AMAX, SCOND
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               S( * )
014:       COMPLEX            AP( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  CPPEQU computes row and column scalings intended to equilibrate a
021: *  Hermitian positive definite matrix A in packed storage and reduce
022: *  its condition number (with respect to the two-norm).  S contains the
023: *  scale factors, S(i)=1/sqrt(A(i,i)), chosen so that the scaled matrix
024: *  B with elements B(i,j)=S(i)*A(i,j)*S(j) has ones on the diagonal.
025: *  This choice of S puts the condition number of B within a factor N of
026: *  the smallest possible condition number over all possible diagonal
027: *  scalings.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  UPLO    (input) CHARACTER*1
033: *          = 'U':  Upper triangle of A is stored;
034: *          = 'L':  Lower triangle of A is stored.
035: *
036: *  N       (input) INTEGER
037: *          The order of the matrix A.  N >= 0.
038: *
039: *  AP      (input) COMPLEX array, dimension (N*(N+1)/2)
040: *          The upper or lower triangle of the Hermitian matrix A, packed
041: *          columnwise in a linear array.  The j-th column of A is stored
042: *          in the array AP as follows:
043: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
044: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
045: *
046: *  S       (output) REAL array, dimension (N)
047: *          If INFO = 0, S contains the scale factors for A.
048: *
049: *  SCOND   (output) REAL
050: *          If INFO = 0, S contains the ratio of the smallest S(i) to
051: *          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
052: *          large nor too small, it is not worth scaling by S.
053: *
054: *  AMAX    (output) REAL
055: *          Absolute value of largest matrix element.  If AMAX is very
056: *          close to overflow or very close to underflow, the matrix
057: *          should be scaled.
058: *
059: *  INFO    (output) INTEGER
060: *          = 0:  successful exit
061: *          < 0:  if INFO = -i, the i-th argument had an illegal value
062: *          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
063: *
064: *  =====================================================================
065: *
066: *     .. Parameters ..
067:       REAL               ONE, ZERO
068:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
069: *     ..
070: *     .. Local Scalars ..
071:       LOGICAL            UPPER
072:       INTEGER            I, JJ
073:       REAL               SMIN
074: *     ..
075: *     .. External Functions ..
076:       LOGICAL            LSAME
077:       EXTERNAL           LSAME
078: *     ..
079: *     .. External Subroutines ..
080:       EXTERNAL           XERBLA
081: *     ..
082: *     .. Intrinsic Functions ..
083:       INTRINSIC          MAX, MIN, REAL, SQRT
084: *     ..
085: *     .. Executable Statements ..
086: *
087: *     Test the input parameters.
088: *
089:       INFO = 0
090:       UPPER = LSAME( UPLO, 'U' )
091:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
092:          INFO = -1
093:       ELSE IF( N.LT.0 ) THEN
094:          INFO = -2
095:       END IF
096:       IF( INFO.NE.0 ) THEN
097:          CALL XERBLA( 'CPPEQU', -INFO )
098:          RETURN
099:       END IF
100: *
101: *     Quick return if possible
102: *
103:       IF( N.EQ.0 ) THEN
104:          SCOND = ONE
105:          AMAX = ZERO
106:          RETURN
107:       END IF
108: *
109: *     Initialize SMIN and AMAX.
110: *
111:       S( 1 ) = REAL( AP( 1 ) )
112:       SMIN = S( 1 )
113:       AMAX = S( 1 )
114: *
115:       IF( UPPER ) THEN
116: *
117: *        UPLO = 'U':  Upper triangle of A is stored.
118: *        Find the minimum and maximum diagonal elements.
119: *
120:          JJ = 1
121:          DO 10 I = 2, N
122:             JJ = JJ + I
123:             S( I ) = REAL( AP( JJ ) )
124:             SMIN = MIN( SMIN, S( I ) )
125:             AMAX = MAX( AMAX, S( I ) )
126:    10    CONTINUE
127: *
128:       ELSE
129: *
130: *        UPLO = 'L':  Lower triangle of A is stored.
131: *        Find the minimum and maximum diagonal elements.
132: *
133:          JJ = 1
134:          DO 20 I = 2, N
135:             JJ = JJ + N - I + 2
136:             S( I ) = REAL( AP( JJ ) )
137:             SMIN = MIN( SMIN, S( I ) )
138:             AMAX = MAX( AMAX, S( I ) )
139:    20    CONTINUE
140:       END IF
141: *
142:       IF( SMIN.LE.ZERO ) THEN
143: *
144: *        Find the first non-positive diagonal element and return.
145: *
146:          DO 30 I = 1, N
147:             IF( S( I ).LE.ZERO ) THEN
148:                INFO = I
149:                RETURN
150:             END IF
151:    30    CONTINUE
152:       ELSE
153: *
154: *        Set the scale factors to the reciprocals
155: *        of the diagonal elements.
156: *
157:          DO 40 I = 1, N
158:             S( I ) = ONE / SQRT( S( I ) )
159:    40    CONTINUE
160: *
161: *        Compute SCOND = min(S(I)) / max(S(I))
162: *
163:          SCOND = SQRT( SMIN ) / SQRT( AMAX )
164:       END IF
165:       RETURN
166: *
167: *     End of CPPEQU
168: *
169:       END
170: