![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLA_SYRPVGRW 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLA_SYRPVGRW + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_syrpvgrw.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_syrpvgrw.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_syrpvgrw.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF, 00022 * LDAF, IPIV, WORK ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER*1 UPLO 00026 * INTEGER N, INFO, LDA, LDAF 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IPIV( * ) 00030 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> 00040 *> DLA_SYRPVGRW computes the reciprocal pivot growth factor 00041 *> norm(A)/norm(U). The "max absolute element" norm is used. If this is 00042 *> much less than 1, the stability of the LU factorization of the 00043 *> (equilibrated) matrix A could be poor. This also means that the 00044 *> solution X, estimated condition numbers, and error bounds could be 00045 *> unreliable. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] UPLO 00052 *> \verbatim 00053 *> UPLO is CHARACTER*1 00054 *> = 'U': Upper triangle of A is stored; 00055 *> = 'L': Lower triangle of A is stored. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] N 00059 *> \verbatim 00060 *> N is INTEGER 00061 *> The number of linear equations, i.e., the order of the 00062 *> matrix A. N >= 0. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] INFO 00066 *> \verbatim 00067 *> INFO is INTEGER 00068 *> The value of INFO returned from DSYTRF, .i.e., the pivot in 00069 *> column INFO is exactly 0. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] A 00073 *> \verbatim 00074 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00075 *> On entry, the N-by-N matrix A. 00076 *> \endverbatim 00077 *> 00078 *> \param[in] LDA 00079 *> \verbatim 00080 *> LDA is INTEGER 00081 *> The leading dimension of the array A. LDA >= max(1,N). 00082 *> \endverbatim 00083 *> 00084 *> \param[in] AF 00085 *> \verbatim 00086 *> AF is DOUBLE PRECISION array, dimension (LDAF,N) 00087 *> The block diagonal matrix D and the multipliers used to 00088 *> obtain the factor U or L as computed by DSYTRF. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] LDAF 00092 *> \verbatim 00093 *> LDAF is INTEGER 00094 *> The leading dimension of the array AF. LDAF >= max(1,N). 00095 *> \endverbatim 00096 *> 00097 *> \param[in] IPIV 00098 *> \verbatim 00099 *> IPIV is INTEGER array, dimension (N) 00100 *> Details of the interchanges and the block structure of D 00101 *> as determined by DSYTRF. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] WORK 00105 *> \verbatim 00106 *> WORK is DOUBLE PRECISION array, dimension (2*N) 00107 *> \endverbatim 00108 * 00109 * Authors: 00110 * ======== 00111 * 00112 *> \author Univ. of Tennessee 00113 *> \author Univ. of California Berkeley 00114 *> \author Univ. of Colorado Denver 00115 *> \author NAG Ltd. 00116 * 00117 *> \date November 2011 00118 * 00119 *> \ingroup doubleSYcomputational 00120 * 00121 * ===================================================================== 00122 DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF, 00123 $ LDAF, IPIV, WORK ) 00124 * 00125 * -- LAPACK computational routine (version 3.4.0) -- 00126 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00127 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00128 * November 2011 00129 * 00130 * .. Scalar Arguments .. 00131 CHARACTER*1 UPLO 00132 INTEGER N, INFO, LDA, LDAF 00133 * .. 00134 * .. Array Arguments .. 00135 INTEGER IPIV( * ) 00136 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ) 00137 * .. 00138 * 00139 * ===================================================================== 00140 * 00141 * .. Local Scalars .. 00142 INTEGER NCOLS, I, J, K, KP 00143 DOUBLE PRECISION AMAX, UMAX, RPVGRW, TMP 00144 LOGICAL UPPER 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC ABS, MAX, MIN 00148 * .. 00149 * .. External Functions .. 00150 EXTERNAL LSAME, DLASET 00151 LOGICAL LSAME 00152 * .. 00153 * .. Executable Statements .. 00154 * 00155 UPPER = LSAME( 'Upper', UPLO ) 00156 IF ( INFO.EQ.0 ) THEN 00157 IF ( UPPER ) THEN 00158 NCOLS = 1 00159 ELSE 00160 NCOLS = N 00161 END IF 00162 ELSE 00163 NCOLS = INFO 00164 END IF 00165 00166 RPVGRW = 1.0D+0 00167 DO I = 1, 2*N 00168 WORK( I ) = 0.0D+0 00169 END DO 00170 * 00171 * Find the max magnitude entry of each column of A. Compute the max 00172 * for all N columns so we can apply the pivot permutation while 00173 * looping below. Assume a full factorization is the common case. 00174 * 00175 IF ( UPPER ) THEN 00176 DO J = 1, N 00177 DO I = 1, J 00178 WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) ) 00179 WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) ) 00180 END DO 00181 END DO 00182 ELSE 00183 DO J = 1, N 00184 DO I = J, N 00185 WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) ) 00186 WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) ) 00187 END DO 00188 END DO 00189 END IF 00190 * 00191 * Now find the max magnitude entry of each column of U or L. Also 00192 * permute the magnitudes of A above so they're in the same order as 00193 * the factor. 00194 * 00195 * The iteration orders and permutations were copied from dsytrs. 00196 * Calls to SSWAP would be severe overkill. 00197 * 00198 IF ( UPPER ) THEN 00199 K = N 00200 DO WHILE ( K .LT. NCOLS .AND. K.GT.0 ) 00201 IF ( IPIV( K ).GT.0 ) THEN 00202 ! 1x1 pivot 00203 KP = IPIV( K ) 00204 IF ( KP .NE. K ) THEN 00205 TMP = WORK( N+K ) 00206 WORK( N+K ) = WORK( N+KP ) 00207 WORK( N+KP ) = TMP 00208 END IF 00209 DO I = 1, K 00210 WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) ) 00211 END DO 00212 K = K - 1 00213 ELSE 00214 ! 2x2 pivot 00215 KP = -IPIV( K ) 00216 TMP = WORK( N+K-1 ) 00217 WORK( N+K-1 ) = WORK( N+KP ) 00218 WORK( N+KP ) = TMP 00219 DO I = 1, K-1 00220 WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) ) 00221 WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) ) 00222 END DO 00223 WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) ) 00224 K = K - 2 00225 END IF 00226 END DO 00227 K = NCOLS 00228 DO WHILE ( K .LE. N ) 00229 IF ( IPIV( K ).GT.0 ) THEN 00230 KP = IPIV( K ) 00231 IF ( KP .NE. K ) THEN 00232 TMP = WORK( N+K ) 00233 WORK( N+K ) = WORK( N+KP ) 00234 WORK( N+KP ) = TMP 00235 END IF 00236 K = K + 1 00237 ELSE 00238 KP = -IPIV( K ) 00239 TMP = WORK( N+K ) 00240 WORK( N+K ) = WORK( N+KP ) 00241 WORK( N+KP ) = TMP 00242 K = K + 2 00243 END IF 00244 END DO 00245 ELSE 00246 K = 1 00247 DO WHILE ( K .LE. NCOLS ) 00248 IF ( IPIV( K ).GT.0 ) THEN 00249 ! 1x1 pivot 00250 KP = IPIV( K ) 00251 IF ( KP .NE. K ) THEN 00252 TMP = WORK( N+K ) 00253 WORK( N+K ) = WORK( N+KP ) 00254 WORK( N+KP ) = TMP 00255 END IF 00256 DO I = K, N 00257 WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) ) 00258 END DO 00259 K = K + 1 00260 ELSE 00261 ! 2x2 pivot 00262 KP = -IPIV( K ) 00263 TMP = WORK( N+K+1 ) 00264 WORK( N+K+1 ) = WORK( N+KP ) 00265 WORK( N+KP ) = TMP 00266 DO I = K+1, N 00267 WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) ) 00268 WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) ) 00269 END DO 00270 WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) ) 00271 K = K + 2 00272 END IF 00273 END DO 00274 K = NCOLS 00275 DO WHILE ( K .GE. 1 ) 00276 IF ( IPIV( K ).GT.0 ) THEN 00277 KP = IPIV( K ) 00278 IF ( KP .NE. K ) THEN 00279 TMP = WORK( N+K ) 00280 WORK( N+K ) = WORK( N+KP ) 00281 WORK( N+KP ) = TMP 00282 END IF 00283 K = K - 1 00284 ELSE 00285 KP = -IPIV( K ) 00286 TMP = WORK( N+K ) 00287 WORK( N+K ) = WORK( N+KP ) 00288 WORK( N+KP ) = TMP 00289 K = K - 2 00290 ENDIF 00291 END DO 00292 END IF 00293 * 00294 * Compute the *inverse* of the max element growth factor. Dividing 00295 * by zero would imply the largest entry of the factor's column is 00296 * zero. Than can happen when either the column of A is zero or 00297 * massive pivots made the factor underflow to zero. Neither counts 00298 * as growth in itself, so simply ignore terms with zero 00299 * denominators. 00300 * 00301 IF ( UPPER ) THEN 00302 DO I = NCOLS, N 00303 UMAX = WORK( I ) 00304 AMAX = WORK( N+I ) 00305 IF ( UMAX /= 0.0D+0 ) THEN 00306 RPVGRW = MIN( AMAX / UMAX, RPVGRW ) 00307 END IF 00308 END DO 00309 ELSE 00310 DO I = 1, NCOLS 00311 UMAX = WORK( I ) 00312 AMAX = WORK( N+I ) 00313 IF ( UMAX /= 0.0D+0 ) THEN 00314 RPVGRW = MIN( AMAX / UMAX, RPVGRW ) 00315 END IF 00316 END DO 00317 END IF 00318 00319 DLA_SYRPVGRW = RPVGRW 00320 END