![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_syrpvgrw.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_syrpvgrw.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION ZLA_SYRPVGRW( 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: COMPLEX*16 A( LDA, * ), AF( LDAF, * )
20: DOUBLE PRECISION WORK( * )
21: INTEGER IPIV( * )
22: * ..
23: *
24: * Purpose
25: * =======
26: *
27: * ZLA_SYRPVGRW 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 ZSYTRF, .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 ZSYTRF.
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 ZSYTRF.
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
77: COMPLEX*16 ZDUM
78: * ..
79: * .. Intrinsic Functions ..
80: INTRINSIC ABS, REAL, DIMAG, MAX, MIN
81: * ..
82: * .. External Subroutines ..
83: EXTERNAL LSAME, ZLASET
84: LOGICAL LSAME
85: * ..
86: * .. Statement Functions ..
87: DOUBLE PRECISION CABS1
88: * ..
89: * .. Statement Function Definitions ..
90: CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
91: * ..
92: * .. Executable Statements ..
93: *
94: UPPER = LSAME( 'Upper', UPLO )
95: IF ( INFO.EQ.0 ) THEN
96: IF ( UPPER ) THEN
97: NCOLS = 1
98: ELSE
99: 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