![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLA_HERPVGRW 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLA_HERPVGRW + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cla_herpvgrw.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_herpvgrw.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_herpvgrw.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * REAL FUNCTION CLA_HERPVGRW( 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 * INTEGER IPIV( * ) 00030 * COMPLEX A( LDA, * ), AF( LDAF, * ) 00031 * REAL WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> 00041 *> CLA_HERPVGRW 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 SSYTRF, .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 CHETRF. 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 CHETRF. 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 complexHEcomputational 00121 * 00122 * ===================================================================== 00123 REAL FUNCTION CLA_HERPVGRW( 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 INTEGER IPIV( * ) 00137 COMPLEX A( LDA, * ), AF( LDAF, * ) 00138 REAL WORK( * ) 00139 * .. 00140 * 00141 * ===================================================================== 00142 * 00143 * .. Local Scalars .. 00144 INTEGER NCOLS, I, J, K, KP 00145 REAL AMAX, UMAX, RPVGRW, TMP 00146 LOGICAL UPPER, LSAME 00147 COMPLEX ZDUM 00148 * .. 00149 * .. External Functions .. 00150 EXTERNAL LSAME, CLASET 00151 * .. 00152 * .. Intrinsic Functions .. 00153 INTRINSIC ABS, REAL, AIMAG, MAX, MIN 00154 * .. 00155 * .. Statement Functions .. 00156 REAL CABS1 00157 * .. 00158 * .. Statement Function Definitions .. 00159 CABS1( ZDUM ) = ABS( REAL ( ZDUM ) ) + ABS( AIMAG ( ZDUM ) ) 00160 * .. 00161 * .. Executable Statements .. 00162 * 00163 UPPER = LSAME( 'Upper', UPLO ) 00164 IF ( INFO.EQ.0 ) THEN 00165 IF (UPPER) THEN 00166 NCOLS = 1 00167 ELSE 00168 NCOLS = N 00169 END IF 00170 ELSE 00171 NCOLS = INFO 00172 END IF 00173 00174 RPVGRW = 1.0 00175 DO I = 1, 2*N 00176 WORK( I ) = 0.0 00177 END DO 00178 * 00179 * Find the max magnitude entry of each column of A. Compute the max 00180 * for all N columns so we can apply the pivot permutation while 00181 * looping below. Assume a full factorization is the common case. 00182 * 00183 IF ( UPPER ) THEN 00184 DO J = 1, N 00185 DO I = 1, J 00186 WORK( N+I ) = MAX( CABS1( A( I,J ) ), WORK( N+I ) ) 00187 WORK( N+J ) = MAX( CABS1( A( I,J ) ), WORK( N+J ) ) 00188 END DO 00189 END DO 00190 ELSE 00191 DO J = 1, N 00192 DO I = J, N 00193 WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) ) 00194 WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) ) 00195 END DO 00196 END DO 00197 END IF 00198 * 00199 * Now find the max magnitude entry of each column of U or L. Also 00200 * permute the magnitudes of A above so they're in the same order as 00201 * the factor. 00202 * 00203 * The iteration orders and permutations were copied from csytrs. 00204 * Calls to SSWAP would be severe overkill. 00205 * 00206 IF ( UPPER ) THEN 00207 K = N 00208 DO WHILE ( K .LT. NCOLS .AND. K.GT.0 ) 00209 IF ( IPIV( K ).GT.0 ) THEN 00210 ! 1x1 pivot 00211 KP = IPIV( K ) 00212 IF ( KP .NE. K ) THEN 00213 TMP = WORK( N+K ) 00214 WORK( N+K ) = WORK( N+KP ) 00215 WORK( N+KP ) = TMP 00216 END IF 00217 DO I = 1, K 00218 WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) ) 00219 END DO 00220 K = K - 1 00221 ELSE 00222 ! 2x2 pivot 00223 KP = -IPIV( K ) 00224 TMP = WORK( N+K-1 ) 00225 WORK( N+K-1 ) = WORK( N+KP ) 00226 WORK( N+KP ) = TMP 00227 DO I = 1, K-1 00228 WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) ) 00229 WORK( K-1 ) = 00230 $ MAX( CABS1( AF( I, K-1 ) ), WORK( K-1 ) ) 00231 END DO 00232 WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) ) 00233 K = K - 2 00234 END IF 00235 END DO 00236 K = NCOLS 00237 DO WHILE ( K .LE. N ) 00238 IF ( IPIV( K ).GT.0 ) THEN 00239 KP = IPIV( K ) 00240 IF ( KP .NE. K ) THEN 00241 TMP = WORK( N+K ) 00242 WORK( N+K ) = WORK( N+KP ) 00243 WORK( N+KP ) = TMP 00244 END IF 00245 K = K + 1 00246 ELSE 00247 KP = -IPIV( K ) 00248 TMP = WORK( N+K ) 00249 WORK( N+K ) = WORK( N+KP ) 00250 WORK( N+KP ) = TMP 00251 K = K + 2 00252 END IF 00253 END DO 00254 ELSE 00255 K = 1 00256 DO WHILE ( K .LE. NCOLS ) 00257 IF ( IPIV( K ).GT.0 ) THEN 00258 ! 1x1 pivot 00259 KP = IPIV( K ) 00260 IF ( KP .NE. K ) THEN 00261 TMP = WORK( N+K ) 00262 WORK( N+K ) = WORK( N+KP ) 00263 WORK( N+KP ) = TMP 00264 END IF 00265 DO I = K, N 00266 WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) ) 00267 END DO 00268 K = K + 1 00269 ELSE 00270 ! 2x2 pivot 00271 KP = -IPIV( K ) 00272 TMP = WORK( N+K+1 ) 00273 WORK( N+K+1 ) = WORK( N+KP ) 00274 WORK( N+KP ) = TMP 00275 DO I = K+1, N 00276 WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) ) 00277 WORK( K+1 ) = 00278 $ MAX( CABS1( AF( I, K+1 ) ) , WORK( K+1 ) ) 00279 END DO 00280 WORK(K) = MAX( CABS1( AF( K, K ) ), WORK( K ) ) 00281 K = K + 2 00282 END IF 00283 END DO 00284 K = NCOLS 00285 DO WHILE ( K .GE. 1 ) 00286 IF ( IPIV( K ).GT.0 ) THEN 00287 KP = IPIV( K ) 00288 IF ( KP .NE. K ) THEN 00289 TMP = WORK( N+K ) 00290 WORK( N+K ) = WORK( N+KP ) 00291 WORK( N+KP ) = TMP 00292 END IF 00293 K = K - 1 00294 ELSE 00295 KP = -IPIV( K ) 00296 TMP = WORK( N+K ) 00297 WORK( N+K ) = WORK( N+KP ) 00298 WORK( N+KP ) = TMP 00299 K = K - 2 00300 ENDIF 00301 END DO 00302 END IF 00303 * 00304 * Compute the *inverse* of the max element growth factor. Dividing 00305 * by zero would imply the largest entry of the factor's column is 00306 * zero. Than can happen when either the column of A is zero or 00307 * massive pivots made the factor underflow to zero. Neither counts 00308 * as growth in itself, so simply ignore terms with zero 00309 * denominators. 00310 * 00311 IF ( UPPER ) THEN 00312 DO I = NCOLS, N 00313 UMAX = WORK( I ) 00314 AMAX = WORK( N+I ) 00315 IF ( UMAX /= 0.0 ) THEN 00316 RPVGRW = MIN( AMAX / UMAX, RPVGRW ) 00317 END IF 00318 END DO 00319 ELSE 00320 DO I = 1, NCOLS 00321 UMAX = WORK( I ) 00322 AMAX = WORK( N+I ) 00323 IF ( UMAX /= 0.0 ) THEN 00324 RPVGRW = MIN( AMAX / UMAX, RPVGRW ) 00325 END IF 00326 END DO 00327 END IF 00328 00329 CLA_HERPVGRW = RPVGRW 00330 END