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