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