1: *> \brief \b DLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian positive-definite 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_PORPVGRW + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_porpvgrw.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_porpvgrw.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_porpvgrw.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * DOUBLE PRECISION FUNCTION DLA_PORPVGRW( UPLO, NCOLS, A, LDA, AF,
22: * LDAF, WORK )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER*1 UPLO
26: * INTEGER NCOLS, LDA, LDAF
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *>
39: *> DLA_PORPVGRW computes the reciprocal pivot growth factor
40: *> norm(A)/norm(U). The "max absolute element" norm is used. If this is
41: *> much less than 1, the stability of the LU factorization of the
42: *> (equilibrated) matrix A could be poor. This also means that the
43: *> solution X, estimated condition numbers, and error bounds could be
44: *> unreliable.
45: *> \endverbatim
46: *
47: * Arguments:
48: * ==========
49: *
50: *> \param[in] UPLO
51: *> \verbatim
52: *> UPLO is CHARACTER*1
53: *> = 'U': Upper triangle of A is stored;
54: *> = 'L': Lower triangle of A is stored.
55: *> \endverbatim
56: *>
57: *> \param[in] NCOLS
58: *> \verbatim
59: *> NCOLS is INTEGER
60: *> The number of columns of the matrix A. NCOLS >= 0.
61: *> \endverbatim
62: *>
63: *> \param[in] A
64: *> \verbatim
65: *> A is DOUBLE PRECISION array, dimension (LDA,N)
66: *> On entry, the N-by-N matrix A.
67: *> \endverbatim
68: *>
69: *> \param[in] LDA
70: *> \verbatim
71: *> LDA is INTEGER
72: *> The leading dimension of the array A. LDA >= max(1,N).
73: *> \endverbatim
74: *>
75: *> \param[in] AF
76: *> \verbatim
77: *> AF is DOUBLE PRECISION array, dimension (LDAF,N)
78: *> The triangular factor U or L from the Cholesky factorization
79: *> A = U**T*U or A = L*L**T, as computed by DPOTRF.
80: *> \endverbatim
81: *>
82: *> \param[in] LDAF
83: *> \verbatim
84: *> LDAF is INTEGER
85: *> The leading dimension of the array AF. LDAF >= max(1,N).
86: *> \endverbatim
87: *>
88: *> \param[out] WORK
89: *> \verbatim
90: *> WORK is DOUBLE PRECISION array, dimension (2*N)
91: *> \endverbatim
92: *
93: * Authors:
94: * ========
95: *
96: *> \author Univ. of Tennessee
97: *> \author Univ. of California Berkeley
98: *> \author Univ. of Colorado Denver
99: *> \author NAG Ltd.
100: *
101: *> \ingroup doublePOcomputational
102: *
103: * =====================================================================
104: DOUBLE PRECISION FUNCTION DLA_PORPVGRW( UPLO, NCOLS, A, LDA, AF,
105: $ LDAF, WORK )
106: *
107: * -- LAPACK computational routine --
108: * -- LAPACK is a software package provided by Univ. of Tennessee, --
109: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
110: *
111: * .. Scalar Arguments ..
112: CHARACTER*1 UPLO
113: INTEGER NCOLS, LDA, LDAF
114: * ..
115: * .. Array Arguments ..
116: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * )
117: * ..
118: *
119: * =====================================================================
120: *
121: * .. Local Scalars ..
122: INTEGER I, J
123: DOUBLE PRECISION AMAX, UMAX, RPVGRW
124: LOGICAL UPPER
125: * ..
126: * .. Intrinsic Functions ..
127: INTRINSIC ABS, MAX, MIN
128: * ..
129: * .. External Functions ..
130: EXTERNAL LSAME
131: LOGICAL LSAME
132: * ..
133: * .. Executable Statements ..
134: *
135: UPPER = LSAME( 'Upper', UPLO )
136: *
137: * DPOTRF will have factored only the NCOLSxNCOLS leading minor, so
138: * we restrict the growth search to that minor and use only the first
139: * 2*NCOLS workspace entries.
140: *
141: RPVGRW = 1.0D+0
142: DO I = 1, 2*NCOLS
143: WORK( I ) = 0.0D+0
144: END DO
145: *
146: * Find the max magnitude entry of each column.
147: *
148: IF ( UPPER ) THEN
149: DO J = 1, NCOLS
150: DO I = 1, J
151: WORK( NCOLS+J ) =
152: $ MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
153: END DO
154: END DO
155: ELSE
156: DO J = 1, NCOLS
157: DO I = J, NCOLS
158: WORK( NCOLS+J ) =
159: $ MAX( ABS( A( I, J ) ), WORK( NCOLS+J ) )
160: END DO
161: END DO
162: END IF
163: *
164: * Now find the max magnitude entry of each column of the factor in
165: * AF. No pivoting, so no permutations.
166: *
167: IF ( LSAME( 'Upper', UPLO ) ) THEN
168: DO J = 1, NCOLS
169: DO I = 1, J
170: WORK( J ) = MAX( ABS( AF( I, J ) ), WORK( J ) )
171: END DO
172: END DO
173: ELSE
174: DO J = 1, NCOLS
175: DO I = J, NCOLS
176: WORK( J ) = MAX( ABS( AF( I, J ) ), WORK( J ) )
177: END DO
178: END DO
179: END IF
180: *
181: * Compute the *inverse* of the max element growth factor. Dividing
182: * by zero would imply the largest entry of the factor's column is
183: * zero. Than can happen when either the column of A is zero or
184: * massive pivots made the factor underflow to zero. Neither counts
185: * as growth in itself, so simply ignore terms with zero
186: * denominators.
187: *
188: IF ( LSAME( 'Upper', UPLO ) ) THEN
189: DO I = 1, NCOLS
190: UMAX = WORK( I )
191: AMAX = WORK( NCOLS+I )
192: IF ( UMAX /= 0.0D+0 ) THEN
193: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
194: END IF
195: END DO
196: ELSE
197: DO I = 1, NCOLS
198: UMAX = WORK( I )
199: AMAX = WORK( NCOLS+I )
200: IF ( UMAX /= 0.0D+0 ) THEN
201: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
202: END IF
203: END DO
204: END IF
205:
206: DLA_PORPVGRW = RPVGRW
207: *
208: * End of DLA_PORPVGRW
209: *
210: END
CVSweb interface <joel.bertrand@systella.fr>