1 DOUBLE PRECISION FUNCTION ZLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF,
2 $ LDAF, IPIV, WORK )
3 *
4 * -- LAPACK routine (version 3.2.2) --
5 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
6 * -- Jason Riedy of Univ. of California Berkeley. --
7 * -- June 2010 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley and NAG Ltd. --
11 *
12 IMPLICIT NONE
13 * ..
14 * .. Scalar Arguments ..
15 CHARACTER*1 UPLO
16 INTEGER N, INFO, LDA, LDAF
17 * ..
18 * .. Array Arguments ..
19 INTEGER IPIV( * )
20 COMPLEX*16 A( LDA, * ), AF( LDAF, * )
21 DOUBLE PRECISION WORK( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * ZLA_HERPVGRW computes the reciprocal pivot growth factor
28 * norm(A)/norm(U). The "max absolute element" norm is used. If this is
29 * much less than 1, the stability of the LU factorization of the
30 * (equilibrated) matrix A could be poor. This also means that the
31 * solution X, estimated condition numbers, and error bounds could be
32 * unreliable.
33 *
34 * Arguments
35 * =========
36 *
37 * UPLO (input) CHARACTER*1
38 * = 'U': Upper triangle of A is stored;
39 * = 'L': Lower triangle of A is stored.
40 *
41 * N (input) INTEGER
42 * The number of linear equations, i.e., the order of the
43 * matrix A. N >= 0.
44 *
45 * INFO (input) INTEGER
46 * The value of INFO returned from ZHETRF, .i.e., the pivot in
47 * column INFO is exactly 0.
48 *
49 * NCOLS (input) INTEGER
50 * The number of columns of the matrix A. NCOLS >= 0.
51 *
52 * A (input) COMPLEX*16 array, dimension (LDA,N)
53 * On entry, the N-by-N matrix A.
54 *
55 * LDA (input) INTEGER
56 * The leading dimension of the array A. LDA >= max(1,N).
57 *
58 * AF (input) COMPLEX*16 array, dimension (LDAF,N)
59 * The block diagonal matrix D and the multipliers used to
60 * obtain the factor U or L as computed by ZHETRF.
61 *
62 * LDAF (input) INTEGER
63 * The leading dimension of the array AF. LDAF >= max(1,N).
64 *
65 * IPIV (input) INTEGER array, dimension (N)
66 * Details of the interchanges and the block structure of D
67 * as determined by ZHETRF.
68 *
69 * WORK (input) COMPLEX*16 array, dimension (2*N)
70 *
71 * =====================================================================
72 *
73 * .. Local Scalars ..
74 INTEGER NCOLS, I, J, K, KP
75 DOUBLE PRECISION AMAX, UMAX, RPVGRW, TMP
76 LOGICAL UPPER, LSAME
77 COMPLEX*16 ZDUM
78 * ..
79 * .. External Functions ..
80 EXTERNAL LSAME, ZLASET
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS, REAL, DIMAG, MAX, MIN
84 * ..
85 * .. Statement Functions ..
86 DOUBLE PRECISION CABS1
87 * ..
88 * .. Statement Function Definitions ..
89 CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
90 * ..
91 * .. Executable Statements ..
92 *
93 UPPER = LSAME( 'Upper', UPLO )
94 IF ( INFO.EQ.0 ) THEN
95 IF (UPPER) THEN
96 NCOLS = 1
97 ELSE
98 NCOLS = N
99 END IF
100 ELSE
101 NCOLS = INFO
102 END IF
103
104 RPVGRW = 1.0D+0
105 DO I = 1, 2*N
106 WORK( I ) = 0.0D+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 zsytrs.
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.0D+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.0D+0 ) THEN
254 RPVGRW = MIN( AMAX / UMAX, RPVGRW )
255 END IF
256 END DO
257 END IF
258
259 ZLA_HERPVGRW = RPVGRW
260 END
2 $ LDAF, IPIV, WORK )
3 *
4 * -- LAPACK routine (version 3.2.2) --
5 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
6 * -- Jason Riedy of Univ. of California Berkeley. --
7 * -- June 2010 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley and NAG Ltd. --
11 *
12 IMPLICIT NONE
13 * ..
14 * .. Scalar Arguments ..
15 CHARACTER*1 UPLO
16 INTEGER N, INFO, LDA, LDAF
17 * ..
18 * .. Array Arguments ..
19 INTEGER IPIV( * )
20 COMPLEX*16 A( LDA, * ), AF( LDAF, * )
21 DOUBLE PRECISION WORK( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * ZLA_HERPVGRW computes the reciprocal pivot growth factor
28 * norm(A)/norm(U). The "max absolute element" norm is used. If this is
29 * much less than 1, the stability of the LU factorization of the
30 * (equilibrated) matrix A could be poor. This also means that the
31 * solution X, estimated condition numbers, and error bounds could be
32 * unreliable.
33 *
34 * Arguments
35 * =========
36 *
37 * UPLO (input) CHARACTER*1
38 * = 'U': Upper triangle of A is stored;
39 * = 'L': Lower triangle of A is stored.
40 *
41 * N (input) INTEGER
42 * The number of linear equations, i.e., the order of the
43 * matrix A. N >= 0.
44 *
45 * INFO (input) INTEGER
46 * The value of INFO returned from ZHETRF, .i.e., the pivot in
47 * column INFO is exactly 0.
48 *
49 * NCOLS (input) INTEGER
50 * The number of columns of the matrix A. NCOLS >= 0.
51 *
52 * A (input) COMPLEX*16 array, dimension (LDA,N)
53 * On entry, the N-by-N matrix A.
54 *
55 * LDA (input) INTEGER
56 * The leading dimension of the array A. LDA >= max(1,N).
57 *
58 * AF (input) COMPLEX*16 array, dimension (LDAF,N)
59 * The block diagonal matrix D and the multipliers used to
60 * obtain the factor U or L as computed by ZHETRF.
61 *
62 * LDAF (input) INTEGER
63 * The leading dimension of the array AF. LDAF >= max(1,N).
64 *
65 * IPIV (input) INTEGER array, dimension (N)
66 * Details of the interchanges and the block structure of D
67 * as determined by ZHETRF.
68 *
69 * WORK (input) COMPLEX*16 array, dimension (2*N)
70 *
71 * =====================================================================
72 *
73 * .. Local Scalars ..
74 INTEGER NCOLS, I, J, K, KP
75 DOUBLE PRECISION AMAX, UMAX, RPVGRW, TMP
76 LOGICAL UPPER, LSAME
77 COMPLEX*16 ZDUM
78 * ..
79 * .. External Functions ..
80 EXTERNAL LSAME, ZLASET
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS, REAL, DIMAG, MAX, MIN
84 * ..
85 * .. Statement Functions ..
86 DOUBLE PRECISION CABS1
87 * ..
88 * .. Statement Function Definitions ..
89 CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
90 * ..
91 * .. Executable Statements ..
92 *
93 UPPER = LSAME( 'Upper', UPLO )
94 IF ( INFO.EQ.0 ) THEN
95 IF (UPPER) THEN
96 NCOLS = 1
97 ELSE
98 NCOLS = N
99 END IF
100 ELSE
101 NCOLS = INFO
102 END IF
103
104 RPVGRW = 1.0D+0
105 DO I = 1, 2*N
106 WORK( I ) = 0.0D+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 zsytrs.
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.0D+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.0D+0 ) THEN
254 RPVGRW = MIN( AMAX / UMAX, RPVGRW )
255 END IF
256 END DO
257 END IF
258
259 ZLA_HERPVGRW = RPVGRW
260 END