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