LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
cla_syrpvgrw.f
Go to the documentation of this file.
1 *> \brief \b CLA_SYRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric indefinite matrix.
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_syrpvgrw.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_syrpvgrw.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_syrpvgrw.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * REAL FUNCTION CLA_SYRPVGRW( 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 * COMPLEX A( LDA, * ), AF( LDAF, * )
30 * REAL WORK( * )
31 * INTEGER IPIV( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *>
41 *> CLA_SYRPVGRW 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 CSYTRF, .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 CSYTRF.
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 CSYTRF.
103 *> \endverbatim
104 *>
105 *> \param[in] WORK
106 *> \verbatim
107 *> WORK is COMPLEX 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 September 2012
119 *
120 *> \ingroup complexSYcomputational
121 *
122 * =====================================================================
123  REAL FUNCTION cla_syrpvgrw( UPLO, N, INFO, A, LDA, AF, LDAF, IPIV,
124  \$ work )
125 *
126 * -- LAPACK computational routine (version 3.4.2) --
127 * -- LAPACK is a software package provided by Univ. of Tennessee, --
128 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
129 * September 2012
130 *
131 * .. Scalar Arguments ..
132  CHARACTER*1 uplo
133  INTEGER n, info, lda, ldaf
134 * ..
135 * .. Array Arguments ..
136  COMPLEX a( lda, * ), af( ldaf, * )
137  REAL work( * )
138  INTEGER ipiv( * )
139 * ..
140 *
141 * =====================================================================
142 *
143 * .. Local Scalars ..
144  INTEGER ncols, i, j, k, kp
145  REAL amax, umax, rpvgrw, tmp
146  LOGICAL upper
147  COMPLEX zdum
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC abs, REAL, aimag, max, min
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL lsame, claset
154  LOGICAL lsame
155 * ..
156 * .. Statement Functions ..
157  REAL cabs1
158 * ..
159 * .. Statement Function Definitions ..
160  cabs1( zdum ) = abs( REAL ( ZDUM ) ) + abs( aimag ( zdum ) )
161 * ..
162 * .. Executable Statements ..
163 *
164  upper = lsame( 'Upper', uplo )
165  IF ( info.EQ.0 ) THEN
166  IF ( upper ) THEN
167  ncols = 1
168  ELSE
169  ncols = n
170  END IF
171  ELSE
172  ncols = info
173  END IF
174
175  rpvgrw = 1.0
176  DO i = 1, 2*n
177  work( i ) = 0.0
178  END DO
179 *
180 * Find the max magnitude entry of each column of A. Compute the max
181 * for all N columns so we can apply the pivot permutation while
182 * looping below. Assume a full factorization is the common case.
183 *
184  IF ( upper ) THEN
185  DO j = 1, n
186  DO i = 1, j
187  work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
188  work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
189  END DO
190  END DO
191  ELSE
192  DO j = 1, n
193  DO i = j, n
194  work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
195  work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
196  END DO
197  END DO
198  END IF
199 *
200 * Now find the max magnitude entry of each column of U or L. Also
201 * permute the magnitudes of A above so they're in the same order as
202 * the factor.
203 *
204 * The iteration orders and permutations were copied from csytrs.
205 * Calls to SSWAP would be severe overkill.
206 *
207  IF ( upper ) THEN
208  k = n
209  DO WHILE ( k .LT. ncols .AND. k.GT.0 )
210  IF ( ipiv( k ).GT.0 ) THEN
211 ! 1x1 pivot
212  kp = ipiv( k )
213  IF ( kp .NE. k ) THEN
214  tmp = work( n+k )
215  work( n+k ) = work( n+kp )
216  work( n+kp ) = tmp
217  END IF
218  DO i = 1, k
219  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
220  END DO
221  k = k - 1
222  ELSE
223 ! 2x2 pivot
224  kp = -ipiv( k )
225  tmp = work( n+k-1 )
226  work( n+k-1 ) = work( n+kp )
227  work( n+kp ) = tmp
228  DO i = 1, k-1
229  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
230  work( k-1 ) =
231  \$ max( cabs1( af( i, k-1 ) ), work( k-1 ) )
232  END DO
233  work( k ) = max( cabs1( af( k, k ) ), work( k ) )
234  k = k - 2
235  END IF
236  END DO
237  k = ncols
238  DO WHILE ( k .LE. n )
239  IF ( ipiv( k ).GT.0 ) THEN
240  kp = ipiv( k )
241  IF ( kp .NE. k ) THEN
242  tmp = work( n+k )
243  work( n+k ) = work( n+kp )
244  work( n+kp ) = tmp
245  END IF
246  k = k + 1
247  ELSE
248  kp = -ipiv( k )
249  tmp = work( n+k )
250  work( n+k ) = work( n+kp )
251  work( n+kp ) = tmp
252  k = k + 2
253  END IF
254  END DO
255  ELSE
256  k = 1
257  DO WHILE ( k .LE. ncols )
258  IF ( ipiv( k ).GT.0 ) THEN
259 ! 1x1 pivot
260  kp = ipiv( k )
261  IF ( kp .NE. k ) THEN
262  tmp = work( n+k )
263  work( n+k ) = work( n+kp )
264  work( n+kp ) = tmp
265  END IF
266  DO i = k, n
267  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
268  END DO
269  k = k + 1
270  ELSE
271 ! 2x2 pivot
272  kp = -ipiv( k )
273  tmp = work( n+k+1 )
274  work( n+k+1 ) = work( n+kp )
275  work( n+kp ) = tmp
276  DO i = k+1, n
277  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
278  work( k+1 ) =
279  \$ max( cabs1( af( i, k+1 ) ), work( k+1 ) )
280  END DO
281  work( k ) = max( cabs1( af( k, k ) ), work( k ) )
282  k = k + 2
283  END IF
284  END DO
285  k = ncols
286  DO WHILE ( k .GE. 1 )
287  IF ( ipiv( k ).GT.0 ) THEN
288  kp = ipiv( k )
289  IF ( kp .NE. k ) THEN
290  tmp = work( n+k )
291  work( n+k ) = work( n+kp )
292  work( n+kp ) = tmp
293  END IF
294  k = k - 1
295  ELSE
296  kp = -ipiv( k )
297  tmp = work( n+k )
298  work( n+k ) = work( n+kp )
299  work( n+kp ) = tmp
300  k = k - 2
301  ENDIF
302  END DO
303  END IF
304 *
305 * Compute the *inverse* of the max element growth factor. Dividing
306 * by zero would imply the largest entry of the factor's column is
307 * zero. Than can happen when either the column of A is zero or
308 * massive pivots made the factor underflow to zero. Neither counts
309 * as growth in itself, so simply ignore terms with zero
310 * denominators.
311 *
312  IF ( upper ) THEN
313  DO i = ncols, n
314  umax = work( i )
315  amax = work( n+i )
316  IF ( umax /= 0.0 ) THEN
317  rpvgrw = min( amax / umax, rpvgrw )
318  END IF
319  END DO
320  ELSE
321  DO i = 1, ncols
322  umax = work( i )
323  amax = work( n+i )
324  IF ( umax /= 0.0 ) THEN
325  rpvgrw = min( amax / umax, rpvgrw )
326  END IF
327  END DO
328  END IF
329
330  cla_syrpvgrw = rpvgrw
331  END