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