LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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
9 *> Download CLA_HERPVGRW + dependencies
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 complexHEcomputational
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function cla_herpvgrw(UPLO, N, INFO, A, LDA, AF, LDAF, IPIV, WORK)
CLA_HERPVGRW
Definition: cla_herpvgrw.f:123