![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_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 dla_syrpvgrw.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION DLA_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: INTEGER IPIV( * )
20: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * DLA_SYRPVGRW computes the reciprocal pivot growth factor
27: * norm(A)/norm(U). The "max absolute element" norm is used. If this is
28: * much less than 1, the stability of the LU factorization of the
29: * (equilibrated) matrix A could be poor. This also means that the
30: * solution X, estimated condition numbers, and error bounds could be
31: * unreliable.
32: *
33: * Arguments
34: * =========
35: *
36: * UPLO (input) CHARACTER*1
37: * = 'U': Upper triangle of A is stored;
38: * = 'L': Lower triangle of A is stored.
39: *
40: * N (input) INTEGER
41: * The number of linear equations, i.e., the order of the
42: * matrix A. N >= 0.
43: *
44: * INFO (input) INTEGER
45: * The value of INFO returned from DSYTRF, .i.e., the pivot in
46: * column INFO is exactly 0.
47: *
48: * NCOLS (input) INTEGER
49: * The number of columns of the matrix A. NCOLS >= 0.
50: *
51: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
52: * On entry, the N-by-N matrix A.
53: *
54: * LDA (input) INTEGER
55: * The leading dimension of the array A. LDA >= max(1,N).
56: *
57: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
58: * The block diagonal matrix D and the multipliers used to
59: * obtain the factor U or L as computed by DSYTRF.
60: *
61: * LDAF (input) INTEGER
62: * The leading dimension of the array AF. LDAF >= max(1,N).
63: *
64: * IPIV (input) INTEGER array, dimension (N)
65: * Details of the interchanges and the block structure of D
66: * as determined by DSYTRF.
67: *
68: * WORK (input) DOUBLE PRECISION array, dimension (2*N)
69: *
70: * =====================================================================
71: *
72: * .. Local Scalars ..
73: INTEGER NCOLS, I, J, K, KP
74: DOUBLE PRECISION AMAX, UMAX, RPVGRW, TMP
75: LOGICAL UPPER
76: * ..
77: * .. Intrinsic Functions ..
78: INTRINSIC ABS, MAX, MIN
79: * ..
80: * .. External Functions ..
81: EXTERNAL LSAME, DLASET
82: LOGICAL LSAME
83: * ..
84: * .. Executable Statements ..
85: *
86: UPPER = LSAME( 'Upper', UPLO )
87: IF ( INFO.EQ.0 ) THEN
88: IF ( UPPER ) THEN
89: NCOLS = 1
90: ELSE
91: NCOLS = N
92: END IF
93: ELSE
94: NCOLS = INFO
95: END IF
96:
97: RPVGRW = 1.0D+0
98: DO I = 1, 2*N
99: WORK( I ) = 0.0D+0
100: END DO
101: *
102: * Find the max magnitude entry of each column of A. Compute the max
103: * for all N columns so we can apply the pivot permutation while
104: * looping below. Assume a full factorization is the common case.
105: *
106: IF ( UPPER ) THEN
107: DO J = 1, N
108: DO I = 1, J
109: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
110: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
111: END DO
112: END DO
113: ELSE
114: DO J = 1, N
115: DO I = J, N
116: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
117: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
118: END DO
119: END DO
120: END IF
121: *
122: * Now find the max magnitude entry of each column of U or L. Also
123: * permute the magnitudes of A above so they're in the same order as
124: * the factor.
125: *
126: * The iteration orders and permutations were copied from dsytrs.
127: * Calls to SSWAP would be severe overkill.
128: *
129: IF ( UPPER ) THEN
130: K = N
131: DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
132: IF ( IPIV( K ).GT.0 ) THEN
133: ! 1x1 pivot
134: KP = IPIV( K )
135: IF ( KP .NE. K ) THEN
136: TMP = WORK( N+K )
137: WORK( N+K ) = WORK( N+KP )
138: WORK( N+KP ) = TMP
139: END IF
140: DO I = 1, K
141: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
142: END DO
143: K = K - 1
144: ELSE
145: ! 2x2 pivot
146: KP = -IPIV( K )
147: TMP = WORK( N+K-1 )
148: WORK( N+K-1 ) = WORK( N+KP )
149: WORK( N+KP ) = TMP
150: DO I = 1, K-1
151: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
152: WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) )
153: END DO
154: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
155: K = K - 2
156: END IF
157: END DO
158: K = NCOLS
159: DO WHILE ( K .LE. N )
160: IF ( IPIV( K ).GT.0 ) THEN
161: KP = IPIV( K )
162: IF ( KP .NE. K ) THEN
163: TMP = WORK( N+K )
164: WORK( N+K ) = WORK( N+KP )
165: WORK( N+KP ) = TMP
166: END IF
167: K = K + 1
168: ELSE
169: KP = -IPIV( K )
170: TMP = WORK( N+K )
171: WORK( N+K ) = WORK( N+KP )
172: WORK( N+KP ) = TMP
173: K = K + 2
174: END IF
175: END DO
176: ELSE
177: K = 1
178: DO WHILE ( K .LE. NCOLS )
179: IF ( IPIV( K ).GT.0 ) THEN
180: ! 1x1 pivot
181: KP = IPIV( K )
182: IF ( KP .NE. K ) THEN
183: TMP = WORK( N+K )
184: WORK( N+K ) = WORK( N+KP )
185: WORK( N+KP ) = TMP
186: END IF
187: DO I = K, N
188: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
189: END DO
190: K = K + 1
191: ELSE
192: ! 2x2 pivot
193: KP = -IPIV( K )
194: TMP = WORK( N+K+1 )
195: WORK( N+K+1 ) = WORK( N+KP )
196: WORK( N+KP ) = TMP
197: DO I = K+1, N
198: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
199: WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) )
200: END DO
201: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
202: K = K + 2
203: END IF
204: END DO
205: K = NCOLS
206: DO WHILE ( K .GE. 1 )
207: IF ( IPIV( K ).GT.0 ) THEN
208: KP = IPIV( K )
209: IF ( KP .NE. K ) THEN
210: TMP = WORK( N+K )
211: WORK( N+K ) = WORK( N+KP )
212: WORK( N+KP ) = TMP
213: END IF
214: K = K - 1
215: ELSE
216: KP = -IPIV( K )
217: TMP = WORK( N+K )
218: WORK( N+K ) = WORK( N+KP )
219: WORK( N+KP ) = TMP
220: K = K - 2
221: ENDIF
222: END DO
223: END IF
224: *
225: * Compute the *inverse* of the max element growth factor. Dividing
226: * by zero would imply the largest entry of the factor's column is
227: * zero. Than can happen when either the column of A is zero or
228: * massive pivots made the factor underflow to zero. Neither counts
229: * as growth in itself, so simply ignore terms with zero
230: * denominators.
231: *
232: IF ( UPPER ) THEN
233: DO I = NCOLS, N
234: UMAX = WORK( I )
235: AMAX = WORK( N+I )
236: IF ( UMAX /= 0.0D+0 ) THEN
237: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
238: END IF
239: END DO
240: ELSE
241: DO I = 1, NCOLS
242: UMAX = WORK( I )
243: AMAX = WORK( N+I )
244: IF ( UMAX /= 0.0D+0 ) THEN
245: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
246: END IF
247: END DO
248: END IF
249:
250: DLA_SYRPVGRW = RPVGRW
251: END