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