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