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