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 321