00001 /* dptsvx.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 00020 /* Subroutine */ int dptsvx_(char *fact, integer *n, integer *nrhs, 00021 doublereal *d__, doublereal *e, doublereal *df, doublereal *ef, 00022 doublereal *b, integer *ldb, doublereal *x, integer *ldx, doublereal * 00023 rcond, doublereal *ferr, doublereal *berr, doublereal *work, integer * 00024 info) 00025 { 00026 /* System generated locals */ 00027 integer b_dim1, b_offset, x_dim1, x_offset, i__1; 00028 00029 /* Local variables */ 00030 extern logical lsame_(char *, char *); 00031 doublereal anorm; 00032 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 00033 doublereal *, integer *); 00034 extern doublereal dlamch_(char *); 00035 logical nofact; 00036 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 00037 doublereal *, integer *, doublereal *, integer *), 00038 xerbla_(char *, integer *); 00039 extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *); 00040 extern /* Subroutine */ int dptcon_(integer *, doublereal *, doublereal *, 00041 doublereal *, doublereal *, doublereal *, integer *), dptrfs_( 00042 integer *, integer *, doublereal *, doublereal *, doublereal *, 00043 doublereal *, doublereal *, integer *, doublereal *, integer *, 00044 doublereal *, doublereal *, doublereal *, integer *), dpttrf_( 00045 integer *, doublereal *, doublereal *, integer *), dpttrs_( 00046 integer *, integer *, doublereal *, doublereal *, doublereal *, 00047 integer *, integer *); 00048 00049 00050 /* -- LAPACK routine (version 3.2) -- */ 00051 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00052 /* November 2006 */ 00053 00054 /* .. Scalar Arguments .. */ 00055 /* .. */ 00056 /* .. Array Arguments .. */ 00057 /* .. */ 00058 00059 /* Purpose */ 00060 /* ======= */ 00061 00062 /* DPTSVX uses the factorization A = L*D*L**T to compute the solution */ 00063 /* to a real system of linear equations A*X = B, where A is an N-by-N */ 00064 /* symmetric positive definite tridiagonal matrix and X and B are */ 00065 /* N-by-NRHS matrices. */ 00066 00067 /* Error bounds on the solution and a condition estimate are also */ 00068 /* provided. */ 00069 00070 /* Description */ 00071 /* =========== */ 00072 00073 /* The following steps are performed: */ 00074 00075 /* 1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L */ 00076 /* is a unit lower bidiagonal matrix and D is diagonal. The */ 00077 /* factorization can also be regarded as having the form */ 00078 /* A = U**T*D*U. */ 00079 00080 /* 2. If the leading i-by-i principal minor is not positive definite, */ 00081 /* then the routine returns with INFO = i. Otherwise, the factored */ 00082 /* form of A is used to estimate the condition number of the matrix */ 00083 /* A. If the reciprocal of the condition number is less than machine */ 00084 /* precision, INFO = N+1 is returned as a warning, but the routine */ 00085 /* still goes on to solve for X and compute error bounds as */ 00086 /* described below. */ 00087 00088 /* 3. The system of equations is solved for X using the factored form */ 00089 /* of A. */ 00090 00091 /* 4. Iterative refinement is applied to improve the computed solution */ 00092 /* matrix and calculate error bounds and backward error estimates */ 00093 /* for it. */ 00094 00095 /* Arguments */ 00096 /* ========= */ 00097 00098 /* FACT (input) CHARACTER*1 */ 00099 /* Specifies whether or not the factored form of A has been */ 00100 /* supplied on entry. */ 00101 /* = 'F': On entry, DF and EF contain the factored form of A. */ 00102 /* D, E, DF, and EF will not be modified. */ 00103 /* = 'N': The matrix A will be copied to DF and EF and */ 00104 /* factored. */ 00105 00106 /* N (input) INTEGER */ 00107 /* The order of the matrix A. N >= 0. */ 00108 00109 /* NRHS (input) INTEGER */ 00110 /* The number of right hand sides, i.e., the number of columns */ 00111 /* of the matrices B and X. NRHS >= 0. */ 00112 00113 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00114 /* The n diagonal elements of the tridiagonal matrix A. */ 00115 00116 /* E (input) DOUBLE PRECISION array, dimension (N-1) */ 00117 /* The (n-1) subdiagonal elements of the tridiagonal matrix A. */ 00118 00119 /* DF (input or output) 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 00127 /* EF (input or output) DOUBLE PRECISION array, dimension (N-1) */ 00128 /* If FACT = 'F', then EF is an input argument and on entry */ 00129 /* contains the (n-1) subdiagonal elements of the unit */ 00130 /* bidiagonal factor L from the L*D*L**T factorization of A. */ 00131 /* If FACT = 'N', then EF is an output argument and on exit */ 00132 /* contains the (n-1) subdiagonal elements of the unit */ 00133 /* bidiagonal factor L from the L*D*L**T factorization of A. */ 00134 00135 /* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00136 /* The N-by-NRHS right hand side matrix B. */ 00137 00138 /* LDB (input) INTEGER */ 00139 /* The leading dimension of the array B. LDB >= max(1,N). */ 00140 00141 /* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) */ 00142 /* If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X. */ 00143 00144 /* LDX (input) INTEGER */ 00145 /* The leading dimension of the array X. LDX >= max(1,N). */ 00146 00147 /* RCOND (output) DOUBLE PRECISION */ 00148 /* The reciprocal condition number of the matrix A. If RCOND */ 00149 /* is less than the machine precision (in particular, if */ 00150 /* RCOND = 0), the matrix is singular to working precision. */ 00151 /* This condition is indicated by a return code of INFO > 0. */ 00152 00153 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00154 /* The forward error bound for each solution vector */ 00155 /* X(j) (the j-th column of the solution matrix X). */ 00156 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00157 /* is an estimated upper bound for the magnitude of the largest */ 00158 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00159 /* largest element in X(j). */ 00160 00161 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00162 /* The componentwise relative backward error of each solution */ 00163 /* vector X(j) (i.e., the smallest relative change in any */ 00164 /* element of A or B that makes X(j) an exact solution). */ 00165 00166 /* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) */ 00167 00168 /* INFO (output) INTEGER */ 00169 /* = 0: successful exit */ 00170 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00171 /* > 0: if INFO = i, and i is */ 00172 /* <= N: the leading minor of order i of A is */ 00173 /* not positive definite, so the factorization */ 00174 /* could not be completed, and the solution has not */ 00175 /* been computed. RCOND = 0 is returned. */ 00176 /* = N+1: U is nonsingular, but RCOND is less than machine */ 00177 /* precision, meaning that the matrix is singular */ 00178 /* to working precision. Nevertheless, the */ 00179 /* solution and error bounds are computed because */ 00180 /* there are a number of situations where the */ 00181 /* computed solution can be more accurate than the */ 00182 /* value of RCOND would suggest. */ 00183 00184 /* ===================================================================== */ 00185 00186 /* .. Parameters .. */ 00187 /* .. */ 00188 /* .. Local Scalars .. */ 00189 /* .. */ 00190 /* .. External Functions .. */ 00191 /* .. */ 00192 /* .. External Subroutines .. */ 00193 /* .. */ 00194 /* .. Intrinsic Functions .. */ 00195 /* .. */ 00196 /* .. Executable Statements .. */ 00197 00198 /* Test the input parameters. */ 00199 00200 /* Parameter adjustments */ 00201 --d__; 00202 --e; 00203 --df; 00204 --ef; 00205 b_dim1 = *ldb; 00206 b_offset = 1 + b_dim1; 00207 b -= b_offset; 00208 x_dim1 = *ldx; 00209 x_offset = 1 + x_dim1; 00210 x -= x_offset; 00211 --ferr; 00212 --berr; 00213 --work; 00214 00215 /* Function Body */ 00216 *info = 0; 00217 nofact = lsame_(fact, "N"); 00218 if (! nofact && ! lsame_(fact, "F")) { 00219 *info = -1; 00220 } else if (*n < 0) { 00221 *info = -2; 00222 } else if (*nrhs < 0) { 00223 *info = -3; 00224 } else if (*ldb < max(1,*n)) { 00225 *info = -9; 00226 } else if (*ldx < max(1,*n)) { 00227 *info = -11; 00228 } 00229 if (*info != 0) { 00230 i__1 = -(*info); 00231 xerbla_("DPTSVX", &i__1); 00232 return 0; 00233 } 00234 00235 if (nofact) { 00236 00237 /* Compute the L*D*L' (or U'*D*U) factorization of A. */ 00238 00239 dcopy_(n, &d__[1], &c__1, &df[1], &c__1); 00240 if (*n > 1) { 00241 i__1 = *n - 1; 00242 dcopy_(&i__1, &e[1], &c__1, &ef[1], &c__1); 00243 } 00244 dpttrf_(n, &df[1], &ef[1], info); 00245 00246 /* Return if INFO is non-zero. */ 00247 00248 if (*info > 0) { 00249 *rcond = 0.; 00250 return 0; 00251 } 00252 } 00253 00254 /* Compute the norm of the matrix A. */ 00255 00256 anorm = dlanst_("1", n, &d__[1], &e[1]); 00257 00258 /* Compute the reciprocal of the condition number of A. */ 00259 00260 dptcon_(n, &df[1], &ef[1], &anorm, rcond, &work[1], info); 00261 00262 /* Compute the solution vectors X. */ 00263 00264 dlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00265 dpttrs_(n, nrhs, &df[1], &ef[1], &x[x_offset], ldx, info); 00266 00267 /* Use iterative refinement to improve the computed solutions and */ 00268 /* compute error bounds and backward error estimates for them. */ 00269 00270 dptrfs_(n, nrhs, &d__[1], &e[1], &df[1], &ef[1], &b[b_offset], ldb, &x[ 00271 x_offset], ldx, &ferr[1], &berr[1], &work[1], info); 00272 00273 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00274 00275 if (*rcond < dlamch_("Epsilon")) { 00276 *info = *n + 1; 00277 } 00278 00279 return 0; 00280 00281 /* End of DPTSVX */ 00282 00283 } /* dptsvx_ */