001:       REAL FUNCTION CLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, IPIV,
002:      $                            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:       COMPLEX            A( LDA, * ), AF( LDAF, * )
021:       REAL               WORK( * )
022: *     ..
023: *
024: *  Purpose
025: *  =======
026: * 
027: *  CLA_HERPVGRW 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 SSYTRF, .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 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 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 CHETRF.
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 CHETRF.
068: *
069: *     WORK    (input) COMPLEX array, dimension (2*N)
070: *
071: *  =====================================================================
072: *
073: *     .. Local Scalars ..
074:       INTEGER            NCOLS, I, J, K, KP
075:       REAL               AMAX, UMAX, RPVGRW, TMP
076:       LOGICAL            UPPER, LSAME
077:       COMPLEX            ZDUM
078: *     ..
079: *     .. External Functions ..
080:       EXTERNAL           LSAME, CLASET
081: *     ..
082: *     .. Intrinsic Functions ..
083:       INTRINSIC          ABS, REAL, AIMAG, MAX, MIN
084: *     ..
085: *     .. Statement Functions ..
086:       REAL               CABS1
087: *     ..
088: *     .. Statement Function Definitions ..
089:       CABS1( ZDUM ) = ABS( REAL ( ZDUM ) ) + ABS( AIMAG ( ZDUM ) )
090: *     ..
091: *     .. Executable Statements ..
092: *
093:       UPPER = LSAME( 'Upper', UPLO )
094:       IF ( INFO.EQ.0 ) THEN
095:          IF (UPPER) THEN
096:             NCOLS = 1
097:          ELSE
098:             NCOLS = N
099:          END IF
100:       ELSE
101:          NCOLS = INFO
102:       END IF
103: 
104:       RPVGRW = 1.0
105:       DO I = 1, 2*N
106:          WORK( I ) = 0.0
107:       END DO
108: *
109: *     Find the max magnitude entry of each column of A.  Compute the max
110: *     for all N columns so we can apply the pivot permutation while
111: *     looping below.  Assume a full factorization is the common case.
112: *
113:       IF ( UPPER ) THEN
114:          DO J = 1, N
115:             DO I = 1, J
116:                WORK( N+I ) = MAX( CABS1( A( I,J ) ), WORK( N+I ) )
117:                WORK( N+J ) = MAX( CABS1( A( I,J ) ), WORK( N+J ) )
118:             END DO
119:          END DO
120:       ELSE
121:          DO J = 1, N
122:             DO I = J, N
123:                WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
124:                WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
125:             END DO
126:          END DO
127:       END IF
128: *
129: *     Now find the max magnitude entry of each column of U or L.  Also
130: *     permute the magnitudes of A above so they're in the same order as
131: *     the factor.
132: *
133: *     The iteration orders and permutations were copied from csytrs.
134: *     Calls to SSWAP would be severe overkill.
135: *
136:       IF ( UPPER ) THEN
137:          K = N
138:          DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
139:             IF ( IPIV( K ).GT.0 ) THEN
140: !              1x1 pivot
141:                KP = IPIV( K )
142:                IF ( KP .NE. K ) THEN
143:                   TMP = WORK( N+K )
144:                   WORK( N+K ) = WORK( N+KP )
145:                   WORK( N+KP ) = TMP
146:                END IF
147:                DO I = 1, K
148:                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
149:                END DO
150:                K = K - 1
151:             ELSE
152: !              2x2 pivot
153:                KP = -IPIV( K )
154:                TMP = WORK( N+K-1 )
155:                WORK( N+K-1 ) = WORK( N+KP )
156:                WORK( N+KP ) = TMP
157:                DO I = 1, K-1
158:                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
159:                   WORK( K-1 ) =
160:      $                 MAX( CABS1( AF( I, K-1 ) ), WORK( K-1 ) )
161:                END DO
162:                WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
163:                K = K - 2
164:             END IF
165:          END DO
166:          K = NCOLS
167:          DO WHILE ( K .LE. N )
168:             IF ( IPIV( K ).GT.0 ) THEN
169:                KP = IPIV( K )
170:                IF ( KP .NE. K ) THEN
171:                   TMP = WORK( N+K )
172:                   WORK( N+K ) = WORK( N+KP )
173:                   WORK( N+KP ) = TMP
174:                END IF
175:                K = K + 1
176:             ELSE
177:                KP = -IPIV( K )
178:                TMP = WORK( N+K )
179:                WORK( N+K ) = WORK( N+KP )
180:                WORK( N+KP ) = TMP
181:                K = K + 2
182:             END IF
183:          END DO
184:       ELSE
185:          K = 1
186:          DO WHILE ( K .LE. NCOLS )
187:             IF ( IPIV( K ).GT.0 ) THEN
188: !              1x1 pivot
189:                KP = IPIV( K )
190:                IF ( KP .NE. K ) THEN
191:                   TMP = WORK( N+K )
192:                   WORK( N+K ) = WORK( N+KP )
193:                   WORK( N+KP ) = TMP
194:                END IF
195:                DO I = K, N
196:                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
197:                END DO
198:                K = K + 1
199:             ELSE
200: !              2x2 pivot
201:                KP = -IPIV( K )
202:                TMP = WORK( N+K+1 )
203:                WORK( N+K+1 ) = WORK( N+KP )
204:                WORK( N+KP ) = TMP
205:                DO I = K+1, N
206:                   WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
207:                   WORK( K+1 ) =
208:      $                 MAX( CABS1( AF( I, K+1 ) ) , WORK( K+1 ) )
209:                END DO
210:                WORK(K) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
211:                K = K + 2
212:             END IF
213:          END DO
214:          K = NCOLS
215:          DO WHILE ( K .GE. 1 )
216:             IF ( IPIV( K ).GT.0 ) THEN
217:                KP = IPIV( K )
218:                IF ( KP .NE. K ) THEN
219:                   TMP = WORK( N+K )
220:                   WORK( N+K ) = WORK( N+KP )
221:                   WORK( N+KP ) = TMP
222:                END IF
223:                K = K - 1
224:             ELSE
225:                KP = -IPIV( K )
226:                TMP = WORK( N+K )
227:                WORK( N+K ) = WORK( N+KP )
228:                WORK( N+KP ) = TMP
229:                K = K - 2
230:             ENDIF
231:          END DO
232:       END IF
233: *
234: *     Compute the *inverse* of the max element growth factor.  Dividing
235: *     by zero would imply the largest entry of the factor's column is
236: *     zero.  Than can happen when either the column of A is zero or
237: *     massive pivots made the factor underflow to zero.  Neither counts
238: *     as growth in itself, so simply ignore terms with zero
239: *     denominators.
240: *
241:       IF ( UPPER ) THEN
242:          DO I = NCOLS, N
243:             UMAX = WORK( I )
244:             AMAX = WORK( N+I )
245:             IF ( UMAX /= 0.0 ) THEN
246:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
247:             END IF
248:          END DO
249:       ELSE
250:          DO I = 1, NCOLS
251:             UMAX = WORK( I )
252:             AMAX = WORK( N+I )
253:             IF ( UMAX /= 0.0 ) THEN
254:                RPVGRW = MIN( AMAX / UMAX, RPVGRW )
255:             END IF
256:          END DO
257:       END IF
258: 
259:       CLA_HERPVGRW = RPVGRW
260:       END FUNCTION
261: