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