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