001:       DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
002:      $                                        LDAF, IPIV, WORK )
003: *
004: *     -- LAPACK routine (version 3.2.1)                                 --
005: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
006: *     -- Jason Riedy of Univ. of California Berkeley.                 --
007: *     -- April 2009                                                   --
008: *
009: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
010: *     -- Univ. of California Berkeley and NAG Ltd.                    --
011: *
012:       IMPLICIT NONE
013: *     ..
014: *     .. Scalar Arguments ..
015:       CHARACTER*1        UPLO
016:       INTEGER            N, INFO, LDA, LDAF
017: *     ..
018: *     .. Array Arguments ..
019:       INTEGER            IPIV( * )
020:       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: * 
026: *  DLA_SYRPVGRW computes the reciprocal pivot growth factor
027: *  norm(A)/norm(U). The "max absolute element" norm is used. If this is
028: *  much less than 1, the stability of the LU factorization of the
029: *  (equilibrated) matrix A could be poor. This also means that the
030: *  solution X, estimated condition numbers, and error bounds could be
031: *  unreliable.
032: *
033: *  Arguments
034: *  =========
035: *
036: *     UPLO    (input) CHARACTER*1
037: *       = 'U':  Upper triangle of A is stored;
038: *       = 'L':  Lower triangle of A is stored.
039: *
040: *     N       (input) INTEGER
041: *     The number of linear equations, i.e., the order of the
042: *     matrix A.  N >= 0.
043: *
044: *     INFO    (input) INTEGER
045: *     The value of INFO returned from DSYTRF, .i.e., the pivot in
046: *     column INFO is exactly 0.
047: *
048: *     NCOLS   (input) INTEGER
049: *     The number of columns of the matrix A. NCOLS >= 0.
050: *
051: *     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
052: *     On entry, the N-by-N matrix A.
053: *
054: *     LDA     (input) INTEGER
055: *     The leading dimension of the array A.  LDA >= max(1,N).
056: *
057: *     AF      (input) DOUBLE PRECISION array, dimension (LDAF,N)
058: *     The block diagonal matrix D and the multipliers used to
059: *     obtain the factor U or L as computed by DSYTRF.
060: *
061: *     LDAF    (input) INTEGER
062: *     The leading dimension of the array AF.  LDAF >= max(1,N).
063: *
064: *     IPIV    (input) INTEGER array, dimension (N)
065: *     Details of the interchanges and the block structure of D
066: *     as determined by DSYTRF.
067: *
068: *     WORK    (input) DOUBLE PRECISION array, dimension (2*N)
069: *
070: *  =====================================================================
071: *
072: *     .. Local Scalars ..
073:       INTEGER            NCOLS, I, J, K, KP
074:       DOUBLE PRECISION   AMAX, UMAX, RPVGRW, TMP
075:       LOGICAL            UPPER
076: *     ..
077: *     .. Intrinsic Functions ..
078:       INTRINSIC          ABS, MAX, MIN
079: *     ..
080: *     .. External Functions ..
081:       EXTERNAL           LSAME, DLASET
082:       LOGICAL            LSAME
083: *     ..
084: *     .. Executable Statements ..
085: *
086:       UPPER = LSAME( 'Upper', UPLO )
087:       IF ( INFO.EQ.0 ) THEN
088:          IF ( UPPER ) THEN
089:             NCOLS = 1
090:          ELSE
091:             NCOLS = N
092:          END IF
093:       ELSE
094:          NCOLS = INFO
095:       END IF
096: 
097:       RPVGRW = 1.0D+0
098:       DO I = 1, 2*N
099:          WORK( I ) = 0.0D+0
100:       END DO
101: *
102: *     Find the max magnitude entry of each column of A.  Compute the max
103: *     for all N columns so we can apply the pivot permutation while
104: *     looping below.  Assume a full factorization is the common case.
105: *
106:       IF ( UPPER ) THEN
107:          DO J = 1, N
108:             DO I = 1, J
109:                WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
110:                WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
111:             END DO
112:          END DO
113:       ELSE
114:          DO J = 1, N
115:             DO I = J, N
116:                WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
117:                WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
118:             END DO
119:          END DO
120:       END IF
121: *
122: *     Now find the max magnitude entry of each column of U or L.  Also
123: *     permute the magnitudes of A above so they're in the same order as
124: *     the factor.
125: *
126: *     The iteration orders and permutations were copied from dsytrs.
127: *     Calls to SSWAP would be severe overkill.
128: *
129:       IF ( UPPER ) THEN
130:          K = N
131:          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
132:             IF ( IPIV( K ).GT.0 ) THEN
133: !              1x1 pivot
134:                KP = IPIV( K )
135:                IF ( KP .NE. K ) THEN
136:                   TMP = WORK( N+K )
137:                   WORK( N+K ) = WORK( N+KP )
138:                   WORK( N+KP ) = TMP
139:                END IF
140:                DO I = 1, K
141:                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
142:                END DO
143:                K = K - 1
144:             ELSE
145: !              2x2 pivot
146:                KP = -IPIV( K )
147:                TMP = WORK( N+K-1 )
148:                WORK( N+K-1 ) = WORK( N+KP )
149:                WORK( N+KP ) = TMP
150:                DO I = 1, K-1
151:                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
152:                   WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) )
153:                END DO
154:                WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
155:                K = K - 2
156:             END IF
157:          END DO
158:          K = NCOLS
159:          DO WHILE ( K .LE. N )
160:             IF ( IPIV( K ).GT.0 ) THEN
161:                KP = IPIV( K )
162:                IF ( KP .NE. K ) THEN
163:                   TMP = WORK( N+K )
164:                   WORK( N+K ) = WORK( N+KP )
165:                   WORK( N+KP ) = TMP
166:                END IF
167:                K = K + 1
168:             ELSE
169:                KP = -IPIV( K )
170:                TMP = WORK( N+K )
171:                WORK( N+K ) = WORK( N+KP )
172:                WORK( N+KP ) = TMP
173:                K = K + 2
174:             END IF
175:          END DO
176:       ELSE
177:          K = 1
178:          DO WHILE ( K .LE. NCOLS )
179:             IF ( IPIV( K ).GT.0 ) THEN
180: !              1x1 pivot
181:                KP = IPIV( K )
182:                IF ( KP .NE. K ) THEN
183:                   TMP = WORK( N+K )
184:                   WORK( N+K ) = WORK( N+KP )
185:                   WORK( N+KP ) = TMP
186:                END IF
187:                DO I = K, N
188:                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
189:                END DO
190:                K = K + 1
191:             ELSE
192: !              2x2 pivot
193:                KP = -IPIV( K )
194:                TMP = WORK( N+K+1 )
195:                WORK( N+K+1 ) = WORK( N+KP )
196:                WORK( N+KP ) = TMP
197:                DO I = K+1, N
198:                   WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
199:                   WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) )
200:                END DO
201:                WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
202:                K = K + 2
203:             END IF
204:          END DO
205:          K = NCOLS
206:          DO WHILE ( K .GE. 1 )
207:             IF ( IPIV( K ).GT.0 ) THEN
208:                KP = IPIV( K )
209:                IF ( KP .NE. K ) THEN
210:                   TMP = WORK( N+K )
211:                   WORK( N+K ) = WORK( N+KP )
212:                   WORK( N+KP ) = TMP
213:                END IF
214:                K = K - 1
215:             ELSE
216:                KP = -IPIV( K )
217:                TMP = WORK( N+K )
218:                WORK( N+K ) = WORK( N+KP )
219:                WORK( N+KP ) = TMP
220:                K = K - 2
221:             ENDIF
222:          END DO
223:       END IF
224: *
225: *     Compute the *inverse* of the max element growth factor.  Dividing
226: *     by zero would imply the largest entry of the factor's column is
227: *     zero.  Than can happen when either the column of A is zero or
228: *     massive pivots made the factor underflow to zero.  Neither counts
229: *     as growth in itself, so simply ignore terms with zero
230: *     denominators.
231: *
232:       IF ( UPPER ) THEN
233:          DO I = NCOLS, N
234:             UMAX = WORK( I )
235:             AMAX = WORK( N+I )
236:             IF ( UMAX /= 0.0D+0 ) THEN
237:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
238:             END IF
239:          END DO
240:       ELSE
241:          DO I = 1, NCOLS
242:             UMAX = WORK( I )
243:             AMAX = WORK( N+I )
244:             IF ( UMAX /= 0.0D+0 ) THEN
245:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
246:             END IF
247:          END DO
248:       END IF
249: 
250:       DLA_SYRPVGRW = RPVGRW
251:       END FUNCTION
252: