![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DPTSVX 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DPTSVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dptsvx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dptsvx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dptsvx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, 00022 * RCOND, FERR, BERR, WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER FACT 00026 * INTEGER INFO, LDB, LDX, N, NRHS 00027 * DOUBLE PRECISION RCOND 00028 * .. 00029 * .. Array Arguments .. 00030 * DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ), 00031 * $ E( * ), EF( * ), FERR( * ), WORK( * ), 00032 * $ X( LDX, * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> DPTSVX uses the factorization A = L*D*L**T to compute the solution 00042 *> to a real system of linear equations A*X = B, where A is an N-by-N 00043 *> symmetric positive definite tridiagonal matrix and X and B are 00044 *> N-by-NRHS matrices. 00045 *> 00046 *> Error bounds on the solution and a condition estimate are also 00047 *> provided. 00048 *> \endverbatim 00049 * 00050 *> \par Description: 00051 * ================= 00052 *> 00053 *> \verbatim 00054 *> 00055 *> The following steps are performed: 00056 *> 00057 *> 1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L 00058 *> is a unit lower bidiagonal matrix and D is diagonal. The 00059 *> factorization can also be regarded as having the form 00060 *> A = U**T*D*U. 00061 *> 00062 *> 2. If the leading i-by-i principal minor is not positive definite, 00063 *> then the routine returns with INFO = i. Otherwise, the factored 00064 *> form of A is used to estimate the condition number of the matrix 00065 *> A. If the reciprocal of the condition number is less than machine 00066 *> precision, INFO = N+1 is returned as a warning, but the routine 00067 *> still goes on to solve for X and compute error bounds as 00068 *> described below. 00069 *> 00070 *> 3. The system of equations is solved for X using the factored form 00071 *> of A. 00072 *> 00073 *> 4. Iterative refinement is applied to improve the computed solution 00074 *> matrix and calculate error bounds and backward error estimates 00075 *> for it. 00076 *> \endverbatim 00077 * 00078 * Arguments: 00079 * ========== 00080 * 00081 *> \param[in] FACT 00082 *> \verbatim 00083 *> FACT is CHARACTER*1 00084 *> Specifies whether or not the factored form of A has been 00085 *> supplied on entry. 00086 *> = 'F': On entry, DF and EF contain the factored form of A. 00087 *> D, E, DF, and EF will not be modified. 00088 *> = 'N': The matrix A will be copied to DF and EF and 00089 *> factored. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] N 00093 *> \verbatim 00094 *> N is INTEGER 00095 *> The order of the matrix A. N >= 0. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] NRHS 00099 *> \verbatim 00100 *> NRHS is INTEGER 00101 *> The number of right hand sides, i.e., the number of columns 00102 *> of the matrices B and X. NRHS >= 0. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] D 00106 *> \verbatim 00107 *> D is DOUBLE PRECISION array, dimension (N) 00108 *> The n diagonal elements of the tridiagonal matrix A. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] E 00112 *> \verbatim 00113 *> E is DOUBLE PRECISION array, dimension (N-1) 00114 *> The (n-1) subdiagonal elements of the tridiagonal matrix A. 00115 *> \endverbatim 00116 *> 00117 *> \param[in,out] DF 00118 *> \verbatim 00119 *> DF is DOUBLE PRECISION array, dimension (N) 00120 *> If FACT = 'F', then DF is an input argument and on entry 00121 *> contains the n diagonal elements of the diagonal matrix D 00122 *> from the L*D*L**T factorization of A. 00123 *> If FACT = 'N', then DF is an output argument and on exit 00124 *> contains the n diagonal elements of the diagonal matrix D 00125 *> from the L*D*L**T factorization of A. 00126 *> \endverbatim 00127 *> 00128 *> \param[in,out] EF 00129 *> \verbatim 00130 *> EF is DOUBLE PRECISION array, dimension (N-1) 00131 *> If FACT = 'F', then EF is an input argument and on entry 00132 *> contains the (n-1) subdiagonal elements of the unit 00133 *> bidiagonal factor L from the L*D*L**T factorization of A. 00134 *> If FACT = 'N', then EF is an output argument and on exit 00135 *> contains the (n-1) subdiagonal elements of the unit 00136 *> bidiagonal factor L from the L*D*L**T factorization of A. 00137 *> \endverbatim 00138 *> 00139 *> \param[in] B 00140 *> \verbatim 00141 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00142 *> The N-by-NRHS right hand side matrix B. 00143 *> \endverbatim 00144 *> 00145 *> \param[in] LDB 00146 *> \verbatim 00147 *> LDB is INTEGER 00148 *> The leading dimension of the array B. LDB >= max(1,N). 00149 *> \endverbatim 00150 *> 00151 *> \param[out] X 00152 *> \verbatim 00153 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00154 *> If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X. 00155 *> \endverbatim 00156 *> 00157 *> \param[in] LDX 00158 *> \verbatim 00159 *> LDX is INTEGER 00160 *> The leading dimension of the array X. LDX >= max(1,N). 00161 *> \endverbatim 00162 *> 00163 *> \param[out] RCOND 00164 *> \verbatim 00165 *> RCOND is DOUBLE PRECISION 00166 *> The reciprocal condition number of the matrix A. If RCOND 00167 *> is less than the machine precision (in particular, if 00168 *> RCOND = 0), the matrix is singular to working precision. 00169 *> This condition is indicated by a return code of INFO > 0. 00170 *> \endverbatim 00171 *> 00172 *> \param[out] FERR 00173 *> \verbatim 00174 *> FERR is DOUBLE PRECISION array, dimension (NRHS) 00175 *> The forward error bound for each solution vector 00176 *> X(j) (the j-th column of the solution matrix X). 00177 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00178 *> is an estimated upper bound for the magnitude of the largest 00179 *> element in (X(j) - XTRUE) divided by the magnitude of the 00180 *> largest element in X(j). 00181 *> \endverbatim 00182 *> 00183 *> \param[out] BERR 00184 *> \verbatim 00185 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00186 *> The componentwise relative backward error of each solution 00187 *> vector X(j) (i.e., the smallest relative change in any 00188 *> element of A or B that makes X(j) an exact solution). 00189 *> \endverbatim 00190 *> 00191 *> \param[out] WORK 00192 *> \verbatim 00193 *> WORK is DOUBLE PRECISION array, dimension (2*N) 00194 *> \endverbatim 00195 *> 00196 *> \param[out] INFO 00197 *> \verbatim 00198 *> INFO is INTEGER 00199 *> = 0: successful exit 00200 *> < 0: if INFO = -i, the i-th argument had an illegal value 00201 *> > 0: if INFO = i, and i is 00202 *> <= N: the leading minor of order i of A is 00203 *> not positive definite, so the factorization 00204 *> could not be completed, and the solution has not 00205 *> been computed. RCOND = 0 is returned. 00206 *> = N+1: U is nonsingular, but RCOND is less than machine 00207 *> precision, meaning that the matrix is singular 00208 *> to working precision. Nevertheless, the 00209 *> solution and error bounds are computed because 00210 *> there are a number of situations where the 00211 *> computed solution can be more accurate than the 00212 *> value of RCOND would suggest. 00213 *> \endverbatim 00214 * 00215 * Authors: 00216 * ======== 00217 * 00218 *> \author Univ. of Tennessee 00219 *> \author Univ. of California Berkeley 00220 *> \author Univ. of Colorado Denver 00221 *> \author NAG Ltd. 00222 * 00223 *> \date April 2012 00224 * 00225 *> \ingroup doubleOTHERcomputational 00226 * 00227 * ===================================================================== 00228 SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, 00229 $ RCOND, FERR, BERR, WORK, INFO ) 00230 * 00231 * -- LAPACK computational routine (version 3.4.1) -- 00232 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00233 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00234 * April 2012 00235 * 00236 * .. Scalar Arguments .. 00237 CHARACTER FACT 00238 INTEGER INFO, LDB, LDX, N, NRHS 00239 DOUBLE PRECISION RCOND 00240 * .. 00241 * .. Array Arguments .. 00242 DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ), 00243 $ E( * ), EF( * ), FERR( * ), WORK( * ), 00244 $ X( LDX, * ) 00245 * .. 00246 * 00247 * ===================================================================== 00248 * 00249 * .. Parameters .. 00250 DOUBLE PRECISION ZERO 00251 PARAMETER ( ZERO = 0.0D+0 ) 00252 * .. 00253 * .. Local Scalars .. 00254 LOGICAL NOFACT 00255 DOUBLE PRECISION ANORM 00256 * .. 00257 * .. External Functions .. 00258 LOGICAL LSAME 00259 DOUBLE PRECISION DLAMCH, DLANST 00260 EXTERNAL LSAME, DLAMCH, DLANST 00261 * .. 00262 * .. External Subroutines .. 00263 EXTERNAL DCOPY, DLACPY, DPTCON, DPTRFS, DPTTRF, DPTTRS, 00264 $ XERBLA 00265 * .. 00266 * .. Intrinsic Functions .. 00267 INTRINSIC MAX 00268 * .. 00269 * .. Executable Statements .. 00270 * 00271 * Test the input parameters. 00272 * 00273 INFO = 0 00274 NOFACT = LSAME( FACT, 'N' ) 00275 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 00276 INFO = -1 00277 ELSE IF( N.LT.0 ) THEN 00278 INFO = -2 00279 ELSE IF( NRHS.LT.0 ) THEN 00280 INFO = -3 00281 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00282 INFO = -9 00283 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00284 INFO = -11 00285 END IF 00286 IF( INFO.NE.0 ) THEN 00287 CALL XERBLA( 'DPTSVX', -INFO ) 00288 RETURN 00289 END IF 00290 * 00291 IF( NOFACT ) THEN 00292 * 00293 * Compute the L*D*L**T (or U**T*D*U) factorization of A. 00294 * 00295 CALL DCOPY( N, D, 1, DF, 1 ) 00296 IF( N.GT.1 ) 00297 $ CALL DCOPY( N-1, E, 1, EF, 1 ) 00298 CALL DPTTRF( N, DF, EF, INFO ) 00299 * 00300 * Return if INFO is non-zero. 00301 * 00302 IF( INFO.GT.0 )THEN 00303 RCOND = ZERO 00304 RETURN 00305 END IF 00306 END IF 00307 * 00308 * Compute the norm of the matrix A. 00309 * 00310 ANORM = DLANST( '1', N, D, E ) 00311 * 00312 * Compute the reciprocal of the condition number of A. 00313 * 00314 CALL DPTCON( N, DF, EF, ANORM, RCOND, WORK, INFO ) 00315 * 00316 * Compute the solution vectors X. 00317 * 00318 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00319 CALL DPTTRS( N, NRHS, DF, EF, X, LDX, INFO ) 00320 * 00321 * Use iterative refinement to improve the computed solutions and 00322 * compute error bounds and backward error estimates for them. 00323 * 00324 CALL DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, 00325 $ WORK, INFO ) 00326 * 00327 * Set INFO = N+1 if the matrix is singular to working precision. 00328 * 00329 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00330 $ INFO = N + 1 00331 * 00332 RETURN 00333 * 00334 * End of DPTSVX 00335 * 00336 END