![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_herpvgrw.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zla_herpvgrw.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 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