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[in] 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: *> \date September 2012
118: *
119: *> \ingroup doubleSYcomputational
120: *
121: * =====================================================================
122: DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
123: $ LDAF, IPIV, WORK )
124: *
125: * -- LAPACK computational routine (version 3.4.2) --
126: * -- LAPACK is a software package provided by Univ. of Tennessee, --
127: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
128: * September 2012
129: *
130: * .. Scalar Arguments ..
131: CHARACTER*1 UPLO
132: INTEGER N, INFO, LDA, LDAF
133: * ..
134: * .. Array Arguments ..
135: INTEGER IPIV( * )
136: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * )
137: * ..
138: *
139: * =====================================================================
140: *
141: * .. Local Scalars ..
142: INTEGER NCOLS, I, J, K, KP
143: DOUBLE PRECISION AMAX, UMAX, RPVGRW, TMP
144: LOGICAL UPPER
145: * ..
146: * .. Intrinsic Functions ..
147: INTRINSIC ABS, MAX, MIN
148: * ..
149: * .. External Functions ..
150: EXTERNAL LSAME, DLASET
151: LOGICAL LSAME
152: * ..
153: * .. Executable Statements ..
154: *
155: UPPER = LSAME( 'Upper', UPLO )
156: IF ( INFO.EQ.0 ) THEN
157: IF ( UPPER ) THEN
158: NCOLS = 1
159: ELSE
160: NCOLS = N
161: END IF
162: ELSE
163: NCOLS = INFO
164: END IF
165:
166: RPVGRW = 1.0D+0
167: DO I = 1, 2*N
168: WORK( I ) = 0.0D+0
169: END DO
170: *
171: * Find the max magnitude entry of each column of A. Compute the max
172: * for all N columns so we can apply the pivot permutation while
173: * looping below. Assume a full factorization is the common case.
174: *
175: IF ( UPPER ) THEN
176: DO J = 1, N
177: DO I = 1, J
178: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
179: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
180: END DO
181: END DO
182: ELSE
183: DO J = 1, N
184: DO I = J, N
185: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
186: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
187: END DO
188: END DO
189: END IF
190: *
191: * Now find the max magnitude entry of each column of U or L. Also
192: * permute the magnitudes of A above so they're in the same order as
193: * the factor.
194: *
195: * The iteration orders and permutations were copied from dsytrs.
196: * Calls to SSWAP would be severe overkill.
197: *
198: IF ( UPPER ) THEN
199: K = N
200: DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
201: IF ( IPIV( K ).GT.0 ) THEN
202: ! 1x1 pivot
203: KP = IPIV( K )
204: IF ( KP .NE. K ) THEN
205: TMP = WORK( N+K )
206: WORK( N+K ) = WORK( N+KP )
207: WORK( N+KP ) = TMP
208: END IF
209: DO I = 1, K
210: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
211: END DO
212: K = K - 1
213: ELSE
214: ! 2x2 pivot
215: KP = -IPIV( K )
216: TMP = WORK( N+K-1 )
217: WORK( N+K-1 ) = WORK( N+KP )
218: WORK( N+KP ) = TMP
219: DO I = 1, K-1
220: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
221: WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) )
222: END DO
223: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
224: K = K - 2
225: END IF
226: END DO
227: K = NCOLS
228: DO WHILE ( K .LE. N )
229: IF ( IPIV( K ).GT.0 ) THEN
230: KP = IPIV( K )
231: IF ( KP .NE. K ) THEN
232: TMP = WORK( N+K )
233: WORK( N+K ) = WORK( N+KP )
234: WORK( N+KP ) = TMP
235: END IF
236: K = K + 1
237: ELSE
238: KP = -IPIV( K )
239: TMP = WORK( N+K )
240: WORK( N+K ) = WORK( N+KP )
241: WORK( N+KP ) = TMP
242: K = K + 2
243: END IF
244: END DO
245: ELSE
246: K = 1
247: DO WHILE ( K .LE. NCOLS )
248: IF ( IPIV( K ).GT.0 ) THEN
249: ! 1x1 pivot
250: KP = IPIV( K )
251: IF ( KP .NE. K ) THEN
252: TMP = WORK( N+K )
253: WORK( N+K ) = WORK( N+KP )
254: WORK( N+KP ) = TMP
255: END IF
256: DO I = K, N
257: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
258: END DO
259: K = K + 1
260: ELSE
261: ! 2x2 pivot
262: KP = -IPIV( K )
263: TMP = WORK( N+K+1 )
264: WORK( N+K+1 ) = WORK( N+KP )
265: WORK( N+KP ) = TMP
266: DO I = K+1, N
267: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
268: WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) )
269: END DO
270: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
271: K = K + 2
272: END IF
273: END DO
274: K = NCOLS
275: DO WHILE ( K .GE. 1 )
276: IF ( IPIV( K ).GT.0 ) THEN
277: KP = IPIV( K )
278: IF ( KP .NE. K ) THEN
279: TMP = WORK( N+K )
280: WORK( N+K ) = WORK( N+KP )
281: WORK( N+KP ) = TMP
282: END IF
283: K = K - 1
284: ELSE
285: KP = -IPIV( K )
286: TMP = WORK( N+K )
287: WORK( N+K ) = WORK( N+KP )
288: WORK( N+KP ) = TMP
289: K = K - 2
290: ENDIF
291: END DO
292: END IF
293: *
294: * Compute the *inverse* of the max element growth factor. Dividing
295: * by zero would imply the largest entry of the factor's column is
296: * zero. Than can happen when either the column of A is zero or
297: * massive pivots made the factor underflow to zero. Neither counts
298: * as growth in itself, so simply ignore terms with zero
299: * denominators.
300: *
301: IF ( UPPER ) THEN
302: DO I = NCOLS, N
303: UMAX = WORK( I )
304: AMAX = WORK( N+I )
305: IF ( UMAX /= 0.0D+0 ) THEN
306: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
307: END IF
308: END DO
309: ELSE
310: DO I = 1, NCOLS
311: UMAX = WORK( I )
312: AMAX = WORK( N+I )
313: IF ( UMAX /= 0.0D+0 ) THEN
314: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
315: END IF
316: END DO
317: END IF
318:
319: DLA_SYRPVGRW = RPVGRW
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