![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLA_GBRCOND 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLA_GBRCOND + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_gbrcond.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_gbrcond.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_gbrcond.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * REAL FUNCTION SLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB, 00022 * IPIV, CMODE, C, INFO, WORK, IWORK ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER TRANS 00026 * INTEGER N, LDAB, LDAFB, INFO, KL, KU, CMODE 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ), IPIV( * ) 00030 * REAL AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ), 00031 * $ C( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> SLA_GBRCOND Estimates the Skeel condition number of op(A) * op2(C) 00041 *> where op2 is determined by CMODE as follows 00042 *> CMODE = 1 op2(C) = C 00043 *> CMODE = 0 op2(C) = I 00044 *> CMODE = -1 op2(C) = inv(C) 00045 *> The Skeel condition number cond(A) = norminf( |inv(A)||A| ) 00046 *> is computed by computing scaling factors R such that 00047 *> diag(R)*A*op2(C) is row equilibrated and computing the standard 00048 *> infinity-norm condition number. 00049 *> \endverbatim 00050 * 00051 * Arguments: 00052 * ========== 00053 * 00054 *> \param[in] TRANS 00055 *> \verbatim 00056 *> TRANS is CHARACTER*1 00057 *> Specifies the form of the system of equations: 00058 *> = 'N': A * X = B (No transpose) 00059 *> = 'T': A**T * X = B (Transpose) 00060 *> = 'C': A**H * X = B (Conjugate Transpose = Transpose) 00061 *> \endverbatim 00062 *> 00063 *> \param[in] N 00064 *> \verbatim 00065 *> N is INTEGER 00066 *> The number of linear equations, i.e., the order of the 00067 *> matrix A. N >= 0. 00068 *> \endverbatim 00069 *> 00070 *> \param[in] KL 00071 *> \verbatim 00072 *> KL is INTEGER 00073 *> The number of subdiagonals within the band of A. KL >= 0. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] KU 00077 *> \verbatim 00078 *> KU is INTEGER 00079 *> The number of superdiagonals within the band of A. KU >= 0. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] AB 00083 *> \verbatim 00084 *> AB is REAL array, dimension (LDAB,N) 00085 *> On entry, the matrix A in band storage, in rows 1 to KL+KU+1. 00086 *> The j-th column of A is stored in the j-th column of the 00087 *> array AB as follows: 00088 *> AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) 00089 *> \endverbatim 00090 *> 00091 *> \param[in] LDAB 00092 *> \verbatim 00093 *> LDAB is INTEGER 00094 *> The leading dimension of the array AB. LDAB >= KL+KU+1. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] AFB 00098 *> \verbatim 00099 *> AFB is REAL array, dimension (LDAFB,N) 00100 *> Details of the LU factorization of the band matrix A, as 00101 *> computed by SGBTRF. U is stored as an upper triangular 00102 *> band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, 00103 *> and the multipliers used during the factorization are stored 00104 *> in rows KL+KU+2 to 2*KL+KU+1. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] LDAFB 00108 *> \verbatim 00109 *> LDAFB is INTEGER 00110 *> The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] IPIV 00114 *> \verbatim 00115 *> IPIV is INTEGER array, dimension (N) 00116 *> The pivot indices from the factorization A = P*L*U 00117 *> as computed by SGBTRF; row i of the matrix was interchanged 00118 *> with row IPIV(i). 00119 *> \endverbatim 00120 *> 00121 *> \param[in] CMODE 00122 *> \verbatim 00123 *> CMODE is INTEGER 00124 *> Determines op2(C) in the formula op(A) * op2(C) as follows: 00125 *> CMODE = 1 op2(C) = C 00126 *> CMODE = 0 op2(C) = I 00127 *> CMODE = -1 op2(C) = inv(C) 00128 *> \endverbatim 00129 *> 00130 *> \param[in] C 00131 *> \verbatim 00132 *> C is REAL array, dimension (N) 00133 *> The vector C in the formula op(A) * op2(C). 00134 *> \endverbatim 00135 *> 00136 *> \param[out] INFO 00137 *> \verbatim 00138 *> INFO is INTEGER 00139 *> = 0: Successful exit. 00140 *> i > 0: The ith argument is invalid. 00141 *> \endverbatim 00142 *> 00143 *> \param[in] WORK 00144 *> \verbatim 00145 *> WORK is REAL array, dimension (5*N). 00146 *> Workspace. 00147 *> \endverbatim 00148 *> 00149 *> \param[in] IWORK 00150 *> \verbatim 00151 *> IWORK is INTEGER array, dimension (N). 00152 *> Workspace. 00153 *> \endverbatim 00154 * 00155 * Authors: 00156 * ======== 00157 * 00158 *> \author Univ. of Tennessee 00159 *> \author Univ. of California Berkeley 00160 *> \author Univ. of Colorado Denver 00161 *> \author NAG Ltd. 00162 * 00163 *> \date November 2011 00164 * 00165 *> \ingroup realGBcomputational 00166 * 00167 * ===================================================================== 00168 REAL FUNCTION SLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB, 00169 $ IPIV, CMODE, C, INFO, WORK, IWORK ) 00170 * 00171 * -- LAPACK computational routine (version 3.4.0) -- 00172 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00174 * November 2011 00175 * 00176 * .. Scalar Arguments .. 00177 CHARACTER TRANS 00178 INTEGER N, LDAB, LDAFB, INFO, KL, KU, CMODE 00179 * .. 00180 * .. Array Arguments .. 00181 INTEGER IWORK( * ), IPIV( * ) 00182 REAL AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ), 00183 $ C( * ) 00184 * .. 00185 * 00186 * ===================================================================== 00187 * 00188 * .. Local Scalars .. 00189 LOGICAL NOTRANS 00190 INTEGER KASE, I, J, KD, KE 00191 REAL AINVNM, TMP 00192 * .. 00193 * .. Local Arrays .. 00194 INTEGER ISAVE( 3 ) 00195 * .. 00196 * .. External Functions .. 00197 LOGICAL LSAME 00198 EXTERNAL LSAME 00199 * .. 00200 * .. External Subroutines .. 00201 EXTERNAL SLACN2, SGBTRS, XERBLA 00202 * .. 00203 * .. Intrinsic Functions .. 00204 INTRINSIC ABS, MAX 00205 * .. 00206 * .. Executable Statements .. 00207 * 00208 SLA_GBRCOND = 0.0 00209 * 00210 INFO = 0 00211 NOTRANS = LSAME( TRANS, 'N' ) 00212 IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') 00213 $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN 00214 INFO = -1 00215 ELSE IF( N.LT.0 ) THEN 00216 INFO = -2 00217 ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN 00218 INFO = -3 00219 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 00220 INFO = -4 00221 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 00222 INFO = -6 00223 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 00224 INFO = -8 00225 END IF 00226 IF( INFO.NE.0 ) THEN 00227 CALL XERBLA( 'SLA_GBRCOND', -INFO ) 00228 RETURN 00229 END IF 00230 IF( N.EQ.0 ) THEN 00231 SLA_GBRCOND = 1.0 00232 RETURN 00233 END IF 00234 * 00235 * Compute the equilibration matrix R such that 00236 * inv(R)*A*C has unit 1-norm. 00237 * 00238 KD = KU + 1 00239 KE = KL + 1 00240 IF ( NOTRANS ) THEN 00241 DO I = 1, N 00242 TMP = 0.0 00243 IF ( CMODE .EQ. 1 ) THEN 00244 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00245 TMP = TMP + ABS( AB( KD+I-J, J ) * C( J ) ) 00246 END DO 00247 ELSE IF ( CMODE .EQ. 0 ) THEN 00248 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00249 TMP = TMP + ABS( AB( KD+I-J, J ) ) 00250 END DO 00251 ELSE 00252 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00253 TMP = TMP + ABS( AB( KD+I-J, J ) / C( J ) ) 00254 END DO 00255 END IF 00256 WORK( 2*N+I ) = TMP 00257 END DO 00258 ELSE 00259 DO I = 1, N 00260 TMP = 0.0 00261 IF ( CMODE .EQ. 1 ) THEN 00262 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00263 TMP = TMP + ABS( AB( KE-I+J, I ) * C( J ) ) 00264 END DO 00265 ELSE IF ( CMODE .EQ. 0 ) THEN 00266 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00267 TMP = TMP + ABS( AB( KE-I+J, I ) ) 00268 END DO 00269 ELSE 00270 DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 00271 TMP = TMP + ABS( AB( KE-I+J, I ) / C( J ) ) 00272 END DO 00273 END IF 00274 WORK( 2*N+I ) = TMP 00275 END DO 00276 END IF 00277 * 00278 * Estimate the norm of inv(op(A)). 00279 * 00280 AINVNM = 0.0 00281 00282 KASE = 0 00283 10 CONTINUE 00284 CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 00285 IF( KASE.NE.0 ) THEN 00286 IF( KASE.EQ.2 ) THEN 00287 * 00288 * Multiply by R. 00289 * 00290 DO I = 1, N 00291 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 00292 END DO 00293 00294 IF ( NOTRANS ) THEN 00295 CALL SGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00296 $ IPIV, WORK, N, INFO ) 00297 ELSE 00298 CALL SGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV, 00299 $ WORK, N, INFO ) 00300 END IF 00301 * 00302 * Multiply by inv(C). 00303 * 00304 IF ( CMODE .EQ. 1 ) THEN 00305 DO I = 1, N 00306 WORK( I ) = WORK( I ) / C( I ) 00307 END DO 00308 ELSE IF ( CMODE .EQ. -1 ) THEN 00309 DO I = 1, N 00310 WORK( I ) = WORK( I ) * C( I ) 00311 END DO 00312 END IF 00313 ELSE 00314 * 00315 * Multiply by inv(C**T). 00316 * 00317 IF ( CMODE .EQ. 1 ) THEN 00318 DO I = 1, N 00319 WORK( I ) = WORK( I ) / C( I ) 00320 END DO 00321 ELSE IF ( CMODE .EQ. -1 ) THEN 00322 DO I = 1, N 00323 WORK( I ) = WORK( I ) * C( I ) 00324 END DO 00325 END IF 00326 00327 IF ( NOTRANS ) THEN 00328 CALL SGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV, 00329 $ WORK, N, INFO ) 00330 ELSE 00331 CALL SGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 00332 $ IPIV, WORK, N, INFO ) 00333 END IF 00334 * 00335 * Multiply by R. 00336 * 00337 DO I = 1, N 00338 WORK( I ) = WORK( I ) * WORK( 2*N+I ) 00339 END DO 00340 END IF 00341 GO TO 10 00342 END IF 00343 * 00344 * Compute the estimate of the reciprocal condition number. 00345 * 00346 IF( AINVNM .NE. 0.0 ) 00347 $ SLA_GBRCOND = ( 1.0 / AINVNM ) 00348 * 00349 RETURN 00350 * 00351 END