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