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