![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGBCON 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZGBCON + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbcon.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbcon.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbcon.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND, 00022 * WORK, RWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER NORM 00026 * INTEGER INFO, KL, KU, LDAB, N 00027 * DOUBLE PRECISION ANORM, RCOND 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER IPIV( * ) 00031 * DOUBLE PRECISION RWORK( * ) 00032 * COMPLEX*16 AB( LDAB, * ), WORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> ZGBCON estimates the reciprocal of the condition number of a complex 00042 *> general band matrix A, in either the 1-norm or the infinity-norm, 00043 *> using the LU factorization computed by ZGBTRF. 00044 *> 00045 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the 00046 *> condition number is computed as 00047 *> RCOND = 1 / ( norm(A) * norm(inv(A)) ). 00048 *> \endverbatim 00049 * 00050 * Arguments: 00051 * ========== 00052 * 00053 *> \param[in] NORM 00054 *> \verbatim 00055 *> NORM is CHARACTER*1 00056 *> Specifies whether the 1-norm condition number or the 00057 *> infinity-norm condition number is required: 00058 *> = '1' or 'O': 1-norm; 00059 *> = 'I': Infinity-norm. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] N 00063 *> \verbatim 00064 *> N is INTEGER 00065 *> The order of the matrix A. N >= 0. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] KL 00069 *> \verbatim 00070 *> KL is INTEGER 00071 *> The number of subdiagonals within the band of A. KL >= 0. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] KU 00075 *> \verbatim 00076 *> KU is INTEGER 00077 *> The number of superdiagonals within the band of A. KU >= 0. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] AB 00081 *> \verbatim 00082 *> AB is COMPLEX*16 array, dimension (LDAB,N) 00083 *> Details of the LU factorization of the band matrix A, as 00084 *> computed by ZGBTRF. U is stored as an upper triangular band 00085 *> matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and 00086 *> the multipliers used during the factorization are stored in 00087 *> rows KL+KU+2 to 2*KL+KU+1. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] LDAB 00091 *> \verbatim 00092 *> LDAB is INTEGER 00093 *> The leading dimension of the array AB. LDAB >= 2*KL+KU+1. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] IPIV 00097 *> \verbatim 00098 *> IPIV is INTEGER array, dimension (N) 00099 *> The pivot indices; for 1 <= i <= N, row i of the matrix was 00100 *> interchanged with row IPIV(i). 00101 *> \endverbatim 00102 *> 00103 *> \param[in] ANORM 00104 *> \verbatim 00105 *> ANORM is DOUBLE PRECISION 00106 *> If NORM = '1' or 'O', the 1-norm of the original matrix A. 00107 *> If NORM = 'I', the infinity-norm of the original matrix A. 00108 *> \endverbatim 00109 *> 00110 *> \param[out] RCOND 00111 *> \verbatim 00112 *> RCOND is DOUBLE PRECISION 00113 *> The reciprocal of the condition number of the matrix A, 00114 *> computed as RCOND = 1/(norm(A) * norm(inv(A))). 00115 *> \endverbatim 00116 *> 00117 *> \param[out] WORK 00118 *> \verbatim 00119 *> WORK is COMPLEX*16 array, dimension (2*N) 00120 *> \endverbatim 00121 *> 00122 *> \param[out] RWORK 00123 *> \verbatim 00124 *> RWORK is DOUBLE PRECISION array, dimension (N) 00125 *> \endverbatim 00126 *> 00127 *> \param[out] INFO 00128 *> \verbatim 00129 *> INFO is INTEGER 00130 *> = 0: successful exit 00131 *> < 0: if INFO = -i, the i-th argument had an illegal value 00132 *> \endverbatim 00133 * 00134 * Authors: 00135 * ======== 00136 * 00137 *> \author Univ. of Tennessee 00138 *> \author Univ. of California Berkeley 00139 *> \author Univ. of Colorado Denver 00140 *> \author NAG Ltd. 00141 * 00142 *> \date November 2011 00143 * 00144 *> \ingroup complex16GBcomputational 00145 * 00146 * ===================================================================== 00147 SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND, 00148 $ WORK, RWORK, INFO ) 00149 * 00150 * -- LAPACK computational routine (version 3.4.0) -- 00151 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00152 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00153 * November 2011 00154 * 00155 * .. Scalar Arguments .. 00156 CHARACTER NORM 00157 INTEGER INFO, KL, KU, LDAB, N 00158 DOUBLE PRECISION ANORM, RCOND 00159 * .. 00160 * .. Array Arguments .. 00161 INTEGER IPIV( * ) 00162 DOUBLE PRECISION RWORK( * ) 00163 COMPLEX*16 AB( LDAB, * ), WORK( * ) 00164 * .. 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 DOUBLE PRECISION ONE, ZERO 00170 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00171 * .. 00172 * .. Local Scalars .. 00173 LOGICAL LNOTI, ONENRM 00174 CHARACTER NORMIN 00175 INTEGER IX, J, JP, KASE, KASE1, KD, LM 00176 DOUBLE PRECISION AINVNM, SCALE, SMLNUM 00177 COMPLEX*16 T, ZDUM 00178 * .. 00179 * .. Local Arrays .. 00180 INTEGER ISAVE( 3 ) 00181 * .. 00182 * .. External Functions .. 00183 LOGICAL LSAME 00184 INTEGER IZAMAX 00185 DOUBLE PRECISION DLAMCH 00186 COMPLEX*16 ZDOTC 00187 EXTERNAL LSAME, IZAMAX, DLAMCH, ZDOTC 00188 * .. 00189 * .. External Subroutines .. 00190 EXTERNAL XERBLA, ZAXPY, ZDRSCL, ZLACN2, ZLATBS 00191 * .. 00192 * .. Intrinsic Functions .. 00193 INTRINSIC ABS, DBLE, DIMAG, MIN 00194 * .. 00195 * .. Statement Functions .. 00196 DOUBLE PRECISION CABS1 00197 * .. 00198 * .. Statement Function definitions .. 00199 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00200 * .. 00201 * .. Executable Statements .. 00202 * 00203 * Test the input parameters. 00204 * 00205 INFO = 0 00206 ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) 00207 IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN 00208 INFO = -1 00209 ELSE IF( N.LT.0 ) THEN 00210 INFO = -2 00211 ELSE IF( KL.LT.0 ) THEN 00212 INFO = -3 00213 ELSE IF( KU.LT.0 ) THEN 00214 INFO = -4 00215 ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN 00216 INFO = -6 00217 ELSE IF( ANORM.LT.ZERO ) THEN 00218 INFO = -8 00219 END IF 00220 IF( INFO.NE.0 ) THEN 00221 CALL XERBLA( 'ZGBCON', -INFO ) 00222 RETURN 00223 END IF 00224 * 00225 * Quick return if possible 00226 * 00227 RCOND = ZERO 00228 IF( N.EQ.0 ) THEN 00229 RCOND = ONE 00230 RETURN 00231 ELSE IF( ANORM.EQ.ZERO ) THEN 00232 RETURN 00233 END IF 00234 * 00235 SMLNUM = DLAMCH( 'Safe minimum' ) 00236 * 00237 * Estimate the norm of inv(A). 00238 * 00239 AINVNM = ZERO 00240 NORMIN = 'N' 00241 IF( ONENRM ) THEN 00242 KASE1 = 1 00243 ELSE 00244 KASE1 = 2 00245 END IF 00246 KD = KL + KU + 1 00247 LNOTI = KL.GT.0 00248 KASE = 0 00249 10 CONTINUE 00250 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00251 IF( KASE.NE.0 ) THEN 00252 IF( KASE.EQ.KASE1 ) THEN 00253 * 00254 * Multiply by inv(L). 00255 * 00256 IF( LNOTI ) THEN 00257 DO 20 J = 1, N - 1 00258 LM = MIN( KL, N-J ) 00259 JP = IPIV( J ) 00260 T = WORK( JP ) 00261 IF( JP.NE.J ) THEN 00262 WORK( JP ) = WORK( J ) 00263 WORK( J ) = T 00264 END IF 00265 CALL ZAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 ) 00266 20 CONTINUE 00267 END IF 00268 * 00269 * Multiply by inv(U). 00270 * 00271 CALL ZLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 00272 $ KL+KU, AB, LDAB, WORK, SCALE, RWORK, INFO ) 00273 ELSE 00274 * 00275 * Multiply by inv(U**H). 00276 * 00277 CALL ZLATBS( 'Upper', 'Conjugate transpose', 'Non-unit', 00278 $ NORMIN, N, KL+KU, AB, LDAB, WORK, SCALE, RWORK, 00279 $ INFO ) 00280 * 00281 * Multiply by inv(L**H). 00282 * 00283 IF( LNOTI ) THEN 00284 DO 30 J = N - 1, 1, -1 00285 LM = MIN( KL, N-J ) 00286 WORK( J ) = WORK( J ) - ZDOTC( LM, AB( KD+1, J ), 1, 00287 $ WORK( J+1 ), 1 ) 00288 JP = IPIV( J ) 00289 IF( JP.NE.J ) THEN 00290 T = WORK( JP ) 00291 WORK( JP ) = WORK( J ) 00292 WORK( J ) = T 00293 END IF 00294 30 CONTINUE 00295 END IF 00296 END IF 00297 * 00298 * Divide X by 1/SCALE if doing so will not cause overflow. 00299 * 00300 NORMIN = 'Y' 00301 IF( SCALE.NE.ONE ) THEN 00302 IX = IZAMAX( N, WORK, 1 ) 00303 IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 00304 $ GO TO 40 00305 CALL ZDRSCL( N, SCALE, WORK, 1 ) 00306 END IF 00307 GO TO 10 00308 END IF 00309 * 00310 * Compute the estimate of the reciprocal condition number. 00311 * 00312 IF( AINVNM.NE.ZERO ) 00313 $ RCOND = ( ONE / AINVNM ) / ANORM 00314 * 00315 40 CONTINUE 00316 RETURN 00317 * 00318 * End of ZGBCON 00319 * 00320 END