LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zla_herpvgrw.f
Go to the documentation of this file.
1 *> \brief \b ZLA_HERPVGRW
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLA_HERPVGRW + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zla_herpvgrw.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zla_herpvgrw.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zla_herpvgrw.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * DOUBLE PRECISION FUNCTION ZLA_HERPVGRW( 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 * COMPLEX*16 A( LDA, * ), AF( LDAF, * )
31 * DOUBLE PRECISION WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *>
41 *> ZLA_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 ZHETRF, .i.e., the pivot in
70 *> column INFO is exactly 0.
71 *> \endverbatim
72 *>
73 *> \param[in] A
74 *> \verbatim
75 *> A is COMPLEX*16 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*16 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 ZHETRF.
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 ZHETRF.
103 *> \endverbatim
104 *>
105 *> \param[in] WORK
106 *> \verbatim
107 *> WORK is COMPLEX*16 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 November 2011
119 *
120 *> \ingroup complex16HEcomputational
121 *
122 * =====================================================================
123  DOUBLE PRECISION FUNCTION zla_herpvgrw( UPLO, N, INFO, A, LDA, AF,
124  $ ldaf, ipiv, work )
125 *
126 * -- LAPACK computational routine (version 3.4.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 * November 2011
130 *
131 * .. Scalar Arguments ..
132  CHARACTER*1 uplo
133  INTEGER n, info, lda, ldaf
134 * ..
135 * .. Array Arguments ..
136  INTEGER ipiv( * )
137  COMPLEX*16 a( lda, * ), af( ldaf, * )
138  DOUBLE PRECISION work( * )
139 * ..
140 *
141 * =====================================================================
142 *
143 * .. Local Scalars ..
144  INTEGER ncols, i, j, k, kp
145  DOUBLE PRECISION amax, umax, rpvgrw, tmp
146  LOGICAL upper, lsame
147  COMPLEX*16 zdum
148 * ..
149 * .. External Functions ..
150  EXTERNAL lsame, zlaset
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC abs, REAL, dimag, max, min
154 * ..
155 * .. Statement Functions ..
156  DOUBLE PRECISION cabs1
157 * ..
158 * .. Statement Function Definitions ..
159  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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.0d+0
175  DO i = 1, 2*n
176  work( i ) = 0.0d+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 zsytrs.
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.0d+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.0d+0 ) THEN
324  rpvgrw = min( amax / umax, rpvgrw )
325  END IF
326  END DO
327  END IF
328 
329  zla_herpvgrw = rpvgrw
330  END