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