|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, 00002 $ RCOND, FERR, BERR, WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER FACT 00011 INTEGER INFO, LDB, LDX, N, NRHS 00012 DOUBLE PRECISION RCOND 00013 * .. 00014 * .. Array Arguments .. 00015 DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ), 00016 $ E( * ), EF( * ), FERR( * ), WORK( * ), 00017 $ X( LDX, * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * DPTSVX uses the factorization A = L*D*L**T to compute the solution 00024 * to a real system of linear equations A*X = B, where A is an N-by-N 00025 * symmetric positive definite tridiagonal matrix and X and B are 00026 * N-by-NRHS matrices. 00027 * 00028 * Error bounds on the solution and a condition estimate are also 00029 * provided. 00030 * 00031 * Description 00032 * =========== 00033 * 00034 * The following steps are performed: 00035 * 00036 * 1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L 00037 * is a unit lower bidiagonal matrix and D is diagonal. The 00038 * factorization can also be regarded as having the form 00039 * A = U**T*D*U. 00040 * 00041 * 2. If the leading i-by-i principal minor is not positive definite, 00042 * then the routine returns with INFO = i. Otherwise, the factored 00043 * form of A is used to estimate the condition number of the matrix 00044 * A. If the reciprocal of the condition number is less than machine 00045 * precision, INFO = N+1 is returned as a warning, but the routine 00046 * still goes on to solve for X and compute error bounds as 00047 * described below. 00048 * 00049 * 3. The system of equations is solved for X using the factored form 00050 * of A. 00051 * 00052 * 4. Iterative refinement is applied to improve the computed solution 00053 * matrix and calculate error bounds and backward error estimates 00054 * for it. 00055 * 00056 * Arguments 00057 * ========= 00058 * 00059 * FACT (input) CHARACTER*1 00060 * Specifies whether or not the factored form of A has been 00061 * supplied on entry. 00062 * = 'F': On entry, DF and EF contain the factored form of A. 00063 * D, E, DF, and EF will not be modified. 00064 * = 'N': The matrix A will be copied to DF and EF and 00065 * factored. 00066 * 00067 * N (input) INTEGER 00068 * The order of the matrix A. N >= 0. 00069 * 00070 * NRHS (input) INTEGER 00071 * The number of right hand sides, i.e., the number of columns 00072 * of the matrices B and X. NRHS >= 0. 00073 * 00074 * D (input) DOUBLE PRECISION array, dimension (N) 00075 * The n diagonal elements of the tridiagonal matrix A. 00076 * 00077 * E (input) DOUBLE PRECISION array, dimension (N-1) 00078 * The (n-1) subdiagonal elements of the tridiagonal matrix A. 00079 * 00080 * DF (input or output) DOUBLE PRECISION array, dimension (N) 00081 * If FACT = 'F', then DF is an input argument and on entry 00082 * contains the n diagonal elements of the diagonal matrix D 00083 * from the L*D*L**T factorization of A. 00084 * If FACT = 'N', then DF is an output argument and on exit 00085 * contains the n diagonal elements of the diagonal matrix D 00086 * from the L*D*L**T factorization of A. 00087 * 00088 * EF (input or output) DOUBLE PRECISION array, dimension (N-1) 00089 * If FACT = 'F', then EF is an input argument and on entry 00090 * contains the (n-1) subdiagonal elements of the unit 00091 * bidiagonal factor L from the L*D*L**T factorization of A. 00092 * If FACT = 'N', then EF is an output argument and on exit 00093 * contains the (n-1) subdiagonal elements of the unit 00094 * bidiagonal factor L from the L*D*L**T factorization of A. 00095 * 00096 * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) 00097 * The N-by-NRHS right hand side matrix B. 00098 * 00099 * LDB (input) INTEGER 00100 * The leading dimension of the array B. LDB >= max(1,N). 00101 * 00102 * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) 00103 * If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X. 00104 * 00105 * LDX (input) INTEGER 00106 * The leading dimension of the array X. LDX >= max(1,N). 00107 * 00108 * RCOND (output) DOUBLE PRECISION 00109 * The reciprocal condition number of the matrix A. If RCOND 00110 * is less than the machine precision (in particular, if 00111 * RCOND = 0), the matrix is singular to working precision. 00112 * This condition is indicated by a return code of INFO > 0. 00113 * 00114 * FERR (output) DOUBLE PRECISION array, dimension (NRHS) 00115 * The forward error bound for each solution vector 00116 * X(j) (the j-th column of the solution matrix X). 00117 * If XTRUE is the true solution corresponding to X(j), FERR(j) 00118 * is an estimated upper bound for the magnitude of the largest 00119 * element in (X(j) - XTRUE) divided by the magnitude of the 00120 * largest element in X(j). 00121 * 00122 * BERR (output) DOUBLE PRECISION array, dimension (NRHS) 00123 * The componentwise relative backward error of each solution 00124 * vector X(j) (i.e., the smallest relative change in any 00125 * element of A or B that makes X(j) an exact solution). 00126 * 00127 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) 00128 * 00129 * INFO (output) INTEGER 00130 * = 0: successful exit 00131 * < 0: if INFO = -i, the i-th argument had an illegal value 00132 * > 0: if INFO = i, and i is 00133 * <= N: the leading minor of order i of A is 00134 * not positive definite, so the factorization 00135 * could not be completed, and the solution has not 00136 * been computed. RCOND = 0 is returned. 00137 * = N+1: U is nonsingular, but RCOND is less than machine 00138 * precision, meaning that the matrix is singular 00139 * to working precision. Nevertheless, the 00140 * solution and error bounds are computed because 00141 * there are a number of situations where the 00142 * computed solution can be more accurate than the 00143 * value of RCOND would suggest. 00144 * 00145 * ===================================================================== 00146 * 00147 * .. Parameters .. 00148 DOUBLE PRECISION ZERO 00149 PARAMETER ( ZERO = 0.0D+0 ) 00150 * .. 00151 * .. Local Scalars .. 00152 LOGICAL NOFACT 00153 DOUBLE PRECISION ANORM 00154 * .. 00155 * .. External Functions .. 00156 LOGICAL LSAME 00157 DOUBLE PRECISION DLAMCH, DLANST 00158 EXTERNAL LSAME, DLAMCH, DLANST 00159 * .. 00160 * .. External Subroutines .. 00161 EXTERNAL DCOPY, DLACPY, DPTCON, DPTRFS, DPTTRF, DPTTRS, 00162 $ XERBLA 00163 * .. 00164 * .. Intrinsic Functions .. 00165 INTRINSIC MAX 00166 * .. 00167 * .. Executable Statements .. 00168 * 00169 * Test the input parameters. 00170 * 00171 INFO = 0 00172 NOFACT = LSAME( FACT, 'N' ) 00173 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 00174 INFO = -1 00175 ELSE IF( N.LT.0 ) THEN 00176 INFO = -2 00177 ELSE IF( NRHS.LT.0 ) THEN 00178 INFO = -3 00179 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00180 INFO = -9 00181 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00182 INFO = -11 00183 END IF 00184 IF( INFO.NE.0 ) THEN 00185 CALL XERBLA( 'DPTSVX', -INFO ) 00186 RETURN 00187 END IF 00188 * 00189 IF( NOFACT ) THEN 00190 * 00191 * Compute the L*D*L' (or U'*D*U) factorization of A. 00192 * 00193 CALL DCOPY( N, D, 1, DF, 1 ) 00194 IF( N.GT.1 ) 00195 $ CALL DCOPY( N-1, E, 1, EF, 1 ) 00196 CALL DPTTRF( N, DF, EF, INFO ) 00197 * 00198 * Return if INFO is non-zero. 00199 * 00200 IF( INFO.GT.0 )THEN 00201 RCOND = ZERO 00202 RETURN 00203 END IF 00204 END IF 00205 * 00206 * Compute the norm of the matrix A. 00207 * 00208 ANORM = DLANST( '1', N, D, E ) 00209 * 00210 * Compute the reciprocal of the condition number of A. 00211 * 00212 CALL DPTCON( N, DF, EF, ANORM, RCOND, WORK, INFO ) 00213 * 00214 * Compute the solution vectors X. 00215 * 00216 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00217 CALL DPTTRS( N, NRHS, DF, EF, X, LDX, INFO ) 00218 * 00219 * Use iterative refinement to improve the computed solutions and 00220 * compute error bounds and backward error estimates for them. 00221 * 00222 CALL DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, 00223 $ WORK, INFO ) 00224 * 00225 * Set INFO = N+1 if the matrix is singular to working precision. 00226 * 00227 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00228 $ INFO = N + 1 00229 * 00230 RETURN 00231 * 00232 * End of DPTSVX 00233 * 00234 END
1.7.3 
