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