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