1: *> \brief \b DLA_SYRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric indefinite matrix.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLA_SYRPVGRW + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_syrpvgrw.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_syrpvgrw.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_syrpvgrw.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
22: * LDAF, IPIV, WORK )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER*1 UPLO
26: * INTEGER N, INFO, LDA, LDAF
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IPIV( * )
30: * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *>
40: *> DLA_SYRPVGRW computes the reciprocal pivot growth factor
41: *> norm(A)/norm(U). The "max absolute element" norm is used. If this is
42: *> much less than 1, the stability of the LU factorization of the
43: *> (equilibrated) matrix A could be poor. This also means that the
44: *> solution X, estimated condition numbers, and error bounds could be
45: *> unreliable.
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] UPLO
52: *> \verbatim
53: *> UPLO is CHARACTER*1
54: *> = 'U': Upper triangle of A is stored;
55: *> = 'L': Lower triangle of A is stored.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The number of linear equations, i.e., the order of the
62: *> matrix A. N >= 0.
63: *> \endverbatim
64: *>
65: *> \param[in] INFO
66: *> \verbatim
67: *> INFO is INTEGER
68: *> The value of INFO returned from DSYTRF, .i.e., the pivot in
69: *> column INFO is exactly 0.
70: *> \endverbatim
71: *>
72: *> \param[in] A
73: *> \verbatim
74: *> A is DOUBLE PRECISION array, dimension (LDA,N)
75: *> On entry, the N-by-N matrix A.
76: *> \endverbatim
77: *>
78: *> \param[in] LDA
79: *> \verbatim
80: *> LDA is INTEGER
81: *> The leading dimension of the array A. LDA >= max(1,N).
82: *> \endverbatim
83: *>
84: *> \param[in] AF
85: *> \verbatim
86: *> AF is DOUBLE PRECISION array, dimension (LDAF,N)
87: *> The block diagonal matrix D and the multipliers used to
88: *> obtain the factor U or L as computed by DSYTRF.
89: *> \endverbatim
90: *>
91: *> \param[in] LDAF
92: *> \verbatim
93: *> LDAF is INTEGER
94: *> The leading dimension of the array AF. LDAF >= max(1,N).
95: *> \endverbatim
96: *>
97: *> \param[in] IPIV
98: *> \verbatim
99: *> IPIV is INTEGER array, dimension (N)
100: *> Details of the interchanges and the block structure of D
101: *> as determined by DSYTRF.
102: *> \endverbatim
103: *>
104: *> \param[out] WORK
105: *> \verbatim
106: *> WORK is DOUBLE PRECISION array, dimension (2*N)
107: *> \endverbatim
108: *
109: * Authors:
110: * ========
111: *
112: *> \author Univ. of Tennessee
113: *> \author Univ. of California Berkeley
114: *> \author Univ. of Colorado Denver
115: *> \author NAG Ltd.
116: *
117: *> \ingroup doubleSYcomputational
118: *
119: * =====================================================================
120: DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
121: $ LDAF, IPIV, WORK )
122: *
123: * -- LAPACK computational routine --
124: * -- LAPACK is a software package provided by Univ. of Tennessee, --
125: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126: *
127: * .. Scalar Arguments ..
128: CHARACTER*1 UPLO
129: INTEGER N, INFO, LDA, LDAF
130: * ..
131: * .. Array Arguments ..
132: INTEGER IPIV( * )
133: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * )
134: * ..
135: *
136: * =====================================================================
137: *
138: * .. Local Scalars ..
139: INTEGER NCOLS, I, J, K, KP
140: DOUBLE PRECISION AMAX, UMAX, RPVGRW, TMP
141: LOGICAL UPPER
142: * ..
143: * .. Intrinsic Functions ..
144: INTRINSIC ABS, MAX, MIN
145: * ..
146: * .. External Functions ..
147: EXTERNAL LSAME
148: LOGICAL LSAME
149: * ..
150: * .. Executable Statements ..
151: *
152: UPPER = LSAME( 'Upper', UPLO )
153: IF ( INFO.EQ.0 ) THEN
154: IF ( UPPER ) THEN
155: NCOLS = 1
156: ELSE
157: NCOLS = N
158: END IF
159: ELSE
160: NCOLS = INFO
161: END IF
162:
163: RPVGRW = 1.0D+0
164: DO I = 1, 2*N
165: WORK( I ) = 0.0D+0
166: END DO
167: *
168: * Find the max magnitude entry of each column of A. Compute the max
169: * for all N columns so we can apply the pivot permutation while
170: * looping below. Assume a full factorization is the common case.
171: *
172: IF ( UPPER ) THEN
173: DO J = 1, N
174: DO I = 1, J
175: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
176: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
177: END DO
178: END DO
179: ELSE
180: DO J = 1, N
181: DO I = J, N
182: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
183: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
184: END DO
185: END DO
186: END IF
187: *
188: * Now find the max magnitude entry of each column of U or L. Also
189: * permute the magnitudes of A above so they're in the same order as
190: * the factor.
191: *
192: * The iteration orders and permutations were copied from dsytrs.
193: * Calls to SSWAP would be severe overkill.
194: *
195: IF ( UPPER ) THEN
196: K = N
197: DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
198: IF ( IPIV( K ).GT.0 ) THEN
199: ! 1x1 pivot
200: KP = IPIV( K )
201: IF ( KP .NE. K ) THEN
202: TMP = WORK( N+K )
203: WORK( N+K ) = WORK( N+KP )
204: WORK( N+KP ) = TMP
205: END IF
206: DO I = 1, K
207: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
208: END DO
209: K = K - 1
210: ELSE
211: ! 2x2 pivot
212: KP = -IPIV( K )
213: TMP = WORK( N+K-1 )
214: WORK( N+K-1 ) = WORK( N+KP )
215: WORK( N+KP ) = TMP
216: DO I = 1, K-1
217: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
218: WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) )
219: END DO
220: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
221: K = K - 2
222: END IF
223: END DO
224: K = NCOLS
225: DO WHILE ( K .LE. N )
226: IF ( IPIV( K ).GT.0 ) THEN
227: KP = IPIV( K )
228: IF ( KP .NE. K ) THEN
229: TMP = WORK( N+K )
230: WORK( N+K ) = WORK( N+KP )
231: WORK( N+KP ) = TMP
232: END IF
233: K = K + 1
234: ELSE
235: KP = -IPIV( K )
236: TMP = WORK( N+K )
237: WORK( N+K ) = WORK( N+KP )
238: WORK( N+KP ) = TMP
239: K = K + 2
240: END IF
241: END DO
242: ELSE
243: K = 1
244: DO WHILE ( K .LE. NCOLS )
245: IF ( IPIV( K ).GT.0 ) THEN
246: ! 1x1 pivot
247: KP = IPIV( K )
248: IF ( KP .NE. K ) THEN
249: TMP = WORK( N+K )
250: WORK( N+K ) = WORK( N+KP )
251: WORK( N+KP ) = TMP
252: END IF
253: DO I = K, N
254: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
255: END DO
256: K = K + 1
257: ELSE
258: ! 2x2 pivot
259: KP = -IPIV( K )
260: TMP = WORK( N+K+1 )
261: WORK( N+K+1 ) = WORK( N+KP )
262: WORK( N+KP ) = TMP
263: DO I = K+1, N
264: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
265: WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) )
266: END DO
267: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
268: K = K + 2
269: END IF
270: END DO
271: K = NCOLS
272: DO WHILE ( K .GE. 1 )
273: IF ( IPIV( K ).GT.0 ) THEN
274: KP = IPIV( K )
275: IF ( KP .NE. K ) THEN
276: TMP = WORK( N+K )
277: WORK( N+K ) = WORK( N+KP )
278: WORK( N+KP ) = TMP
279: END IF
280: K = K - 1
281: ELSE
282: KP = -IPIV( K )
283: TMP = WORK( N+K )
284: WORK( N+K ) = WORK( N+KP )
285: WORK( N+KP ) = TMP
286: K = K - 2
287: ENDIF
288: END DO
289: END IF
290: *
291: * Compute the *inverse* of the max element growth factor. Dividing
292: * by zero would imply the largest entry of the factor's column is
293: * zero. Than can happen when either the column of A is zero or
294: * massive pivots made the factor underflow to zero. Neither counts
295: * as growth in itself, so simply ignore terms with zero
296: * denominators.
297: *
298: IF ( UPPER ) THEN
299: DO I = NCOLS, N
300: UMAX = WORK( I )
301: AMAX = WORK( N+I )
302: IF ( UMAX /= 0.0D+0 ) THEN
303: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
304: END IF
305: END DO
306: ELSE
307: DO I = 1, NCOLS
308: UMAX = WORK( I )
309: AMAX = WORK( N+I )
310: IF ( UMAX /= 0.0D+0 ) THEN
311: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
312: END IF
313: END DO
314: END IF
315:
316: DLA_SYRPVGRW = RPVGRW
317: *
318: * End of DLA_SYRPVGRW
319: *
320: END
CVSweb interface <joel.bertrand@systella.fr>
![[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