![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SPBCON 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SPBCON + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spbcon.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spbcon.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spbcon.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK, 00022 * IWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO 00026 * INTEGER INFO, KD, LDAB, N 00027 * REAL ANORM, RCOND 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER IWORK( * ) 00031 * REAL AB( LDAB, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> SPBCON estimates the reciprocal of the condition number (in the 00041 *> 1-norm) of a real symmetric positive definite band matrix using the 00042 *> Cholesky factorization A = U**T*U or A = L*L**T computed by SPBTRF. 00043 *> 00044 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the 00045 *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] UPLO 00052 *> \verbatim 00053 *> UPLO is CHARACTER*1 00054 *> = 'U': Upper triangular factor stored in AB; 00055 *> = 'L': Lower triangular factor stored in AB. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] N 00059 *> \verbatim 00060 *> N is INTEGER 00061 *> The order of the matrix A. N >= 0. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] KD 00065 *> \verbatim 00066 *> KD is INTEGER 00067 *> The number of superdiagonals of the matrix A if UPLO = 'U', 00068 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] AB 00072 *> \verbatim 00073 *> AB is REAL array, dimension (LDAB,N) 00074 *> The triangular factor U or L from the Cholesky factorization 00075 *> A = U**T*U or A = L*L**T of the band matrix A, stored in the 00076 *> first KD+1 rows of the array. The j-th column of U or L is 00077 *> stored in the j-th column of the array AB as follows: 00078 *> if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j; 00079 *> if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd). 00080 *> \endverbatim 00081 *> 00082 *> \param[in] LDAB 00083 *> \verbatim 00084 *> LDAB is INTEGER 00085 *> The leading dimension of the array AB. LDAB >= KD+1. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] ANORM 00089 *> \verbatim 00090 *> ANORM is REAL 00091 *> The 1-norm (or infinity-norm) of the symmetric band matrix A. 00092 *> \endverbatim 00093 *> 00094 *> \param[out] RCOND 00095 *> \verbatim 00096 *> RCOND is REAL 00097 *> The reciprocal of the condition number of the matrix A, 00098 *> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 00099 *> estimate of the 1-norm of inv(A) computed in this routine. 00100 *> \endverbatim 00101 *> 00102 *> \param[out] WORK 00103 *> \verbatim 00104 *> WORK is REAL array, dimension (3*N) 00105 *> \endverbatim 00106 *> 00107 *> \param[out] IWORK 00108 *> \verbatim 00109 *> IWORK is INTEGER array, dimension (N) 00110 *> \endverbatim 00111 *> 00112 *> \param[out] INFO 00113 *> \verbatim 00114 *> INFO is INTEGER 00115 *> = 0: successful exit 00116 *> < 0: if INFO = -i, the i-th argument had an illegal value 00117 *> \endverbatim 00118 * 00119 * Authors: 00120 * ======== 00121 * 00122 *> \author Univ. of Tennessee 00123 *> \author Univ. of California Berkeley 00124 *> \author Univ. of Colorado Denver 00125 *> \author NAG Ltd. 00126 * 00127 *> \date November 2011 00128 * 00129 *> \ingroup realOTHERcomputational 00130 * 00131 * ===================================================================== 00132 SUBROUTINE SPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK, 00133 $ IWORK, INFO ) 00134 * 00135 * -- LAPACK computational routine (version 3.4.0) -- 00136 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00138 * November 2011 00139 * 00140 * .. Scalar Arguments .. 00141 CHARACTER UPLO 00142 INTEGER INFO, KD, LDAB, N 00143 REAL ANORM, RCOND 00144 * .. 00145 * .. Array Arguments .. 00146 INTEGER IWORK( * ) 00147 REAL AB( LDAB, * ), WORK( * ) 00148 * .. 00149 * 00150 * ===================================================================== 00151 * 00152 * .. Parameters .. 00153 REAL ONE, ZERO 00154 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00155 * .. 00156 * .. Local Scalars .. 00157 LOGICAL UPPER 00158 CHARACTER NORMIN 00159 INTEGER IX, KASE 00160 REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 00161 * .. 00162 * .. Local Arrays .. 00163 INTEGER ISAVE( 3 ) 00164 * .. 00165 * .. External Functions .. 00166 LOGICAL LSAME 00167 INTEGER ISAMAX 00168 REAL SLAMCH 00169 EXTERNAL LSAME, ISAMAX, SLAMCH 00170 * .. 00171 * .. External Subroutines .. 00172 EXTERNAL SLACN2, SLATBS, SRSCL, XERBLA 00173 * .. 00174 * .. Intrinsic Functions .. 00175 INTRINSIC ABS 00176 * .. 00177 * .. Executable Statements .. 00178 * 00179 * Test the input parameters. 00180 * 00181 INFO = 0 00182 UPPER = LSAME( UPLO, 'U' ) 00183 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00184 INFO = -1 00185 ELSE IF( N.LT.0 ) THEN 00186 INFO = -2 00187 ELSE IF( KD.LT.0 ) THEN 00188 INFO = -3 00189 ELSE IF( LDAB.LT.KD+1 ) THEN 00190 INFO = -5 00191 ELSE IF( ANORM.LT.ZERO ) THEN 00192 INFO = -6 00193 END IF 00194 IF( INFO.NE.0 ) THEN 00195 CALL XERBLA( 'SPBCON', -INFO ) 00196 RETURN 00197 END IF 00198 * 00199 * Quick return if possible 00200 * 00201 RCOND = ZERO 00202 IF( N.EQ.0 ) THEN 00203 RCOND = ONE 00204 RETURN 00205 ELSE IF( ANORM.EQ.ZERO ) THEN 00206 RETURN 00207 END IF 00208 * 00209 SMLNUM = SLAMCH( 'Safe minimum' ) 00210 * 00211 * Estimate the 1-norm of the inverse. 00212 * 00213 KASE = 0 00214 NORMIN = 'N' 00215 10 CONTINUE 00216 CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 00217 IF( KASE.NE.0 ) THEN 00218 IF( UPPER ) THEN 00219 * 00220 * Multiply by inv(U**T). 00221 * 00222 CALL SLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, 00223 $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ), 00224 $ INFO ) 00225 NORMIN = 'Y' 00226 * 00227 * Multiply by inv(U). 00228 * 00229 CALL SLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 00230 $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ), 00231 $ INFO ) 00232 ELSE 00233 * 00234 * Multiply by inv(L). 00235 * 00236 CALL SLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 00237 $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ), 00238 $ INFO ) 00239 NORMIN = 'Y' 00240 * 00241 * Multiply by inv(L**T). 00242 * 00243 CALL SLATBS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, 00244 $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ), 00245 $ INFO ) 00246 END IF 00247 * 00248 * Multiply by 1/SCALE if doing so will not cause overflow. 00249 * 00250 SCALE = SCALEL*SCALEU 00251 IF( SCALE.NE.ONE ) THEN 00252 IX = ISAMAX( N, WORK, 1 ) 00253 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 00254 $ GO TO 20 00255 CALL SRSCL( N, SCALE, WORK, 1 ) 00256 END IF 00257 GO TO 10 00258 END IF 00259 * 00260 * Compute the estimate of the reciprocal condition number. 00261 * 00262 IF( AINVNM.NE.ZERO ) 00263 $ RCOND = ( ONE / AINVNM ) / ANORM 00264 * 00265 20 CONTINUE 00266 * 00267 RETURN 00268 * 00269 * End of SPBCON 00270 * 00271 END