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