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