00001 /* dspsvx.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 dspsvx_(char *fact, char *uplo, integer *n, integer * 00021 nrhs, doublereal *ap, doublereal *afp, integer *ipiv, doublereal *b, 00022 integer *ldb, doublereal *x, integer *ldx, doublereal *rcond, 00023 doublereal *ferr, doublereal *berr, doublereal *work, integer *iwork, 00024 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 *); 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 dlansp_(char *, char *, integer *, doublereal *, 00040 doublereal *); 00041 extern /* Subroutine */ int dspcon_(char *, integer *, doublereal *, 00042 integer *, doublereal *, doublereal *, doublereal *, integer *, 00043 integer *), dsprfs_(char *, integer *, integer *, 00044 doublereal *, doublereal *, integer *, doublereal *, integer *, 00045 doublereal *, integer *, doublereal *, doublereal *, doublereal *, 00046 integer *, integer *), dsptrf_(char *, integer *, 00047 doublereal *, integer *, integer *), dsptrs_(char *, 00048 integer *, integer *, doublereal *, integer *, doublereal *, 00049 integer *, integer *); 00050 00051 00052 /* -- LAPACK driver routine (version 3.2) -- */ 00053 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00054 /* November 2006 */ 00055 00056 /* .. Scalar Arguments .. */ 00057 /* .. */ 00058 /* .. Array Arguments .. */ 00059 /* .. */ 00060 00061 /* Purpose */ 00062 /* ======= */ 00063 00064 /* DSPSVX uses the diagonal pivoting factorization A = U*D*U**T or */ 00065 /* A = L*D*L**T to compute the solution to a real system of linear */ 00066 /* equations A * X = B, where A is an N-by-N symmetric matrix stored */ 00067 /* in packed format and X and B are N-by-NRHS matrices. */ 00068 00069 /* Error bounds on the solution and a condition estimate are also */ 00070 /* provided. */ 00071 00072 /* Description */ 00073 /* =========== */ 00074 00075 /* The following steps are performed: */ 00076 00077 /* 1. If FACT = 'N', the diagonal pivoting method is used to factor A as */ 00078 /* A = U * D * U**T, if UPLO = 'U', or */ 00079 /* A = L * D * L**T, if UPLO = 'L', */ 00080 /* where U (or L) is a product of permutation and unit upper (lower) */ 00081 /* triangular matrices and D is symmetric and block diagonal with */ 00082 /* 1-by-1 and 2-by-2 diagonal blocks. */ 00083 00084 /* 2. If some D(i,i)=0, so that D is exactly singular, then the routine */ 00085 /* returns with INFO = i. Otherwise, the factored form of A is used */ 00086 /* to estimate the condition number of the matrix A. If the */ 00087 /* reciprocal of the condition number is less than machine precision, */ 00088 /* INFO = N+1 is returned as a warning, but the routine still goes on */ 00089 /* to solve for X and compute error bounds as 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 A has been */ 00103 /* supplied on entry. */ 00104 /* = 'F': On entry, AFP and IPIV contain the factored form of */ 00105 /* A. AP, AFP and IPIV will not be modified. */ 00106 /* = 'N': The matrix A will be copied to AFP and factored. */ 00107 00108 /* UPLO (input) CHARACTER*1 */ 00109 /* = 'U': Upper triangle of A is stored; */ 00110 /* = 'L': Lower triangle of A is stored. */ 00111 00112 /* N (input) INTEGER */ 00113 /* The number of linear equations, i.e., the order of the */ 00114 /* matrix A. N >= 0. */ 00115 00116 /* NRHS (input) INTEGER */ 00117 /* The number of right hand sides, i.e., the number of columns */ 00118 /* of the matrices B and X. NRHS >= 0. */ 00119 00120 /* AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2) */ 00121 /* The upper or lower triangle of the symmetric matrix A, packed */ 00122 /* columnwise in a linear array. The j-th column of A is stored */ 00123 /* in the array AP as follows: */ 00124 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */ 00125 /* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */ 00126 /* See below for further details. */ 00127 00128 /* AFP (input or output) DOUBLE PRECISION array, dimension */ 00129 /* (N*(N+1)/2) */ 00130 /* If FACT = 'F', then AFP is an input argument and on entry */ 00131 /* contains the block diagonal matrix D and the multipliers used */ 00132 /* to obtain the factor U or L from the factorization */ 00133 /* A = U*D*U**T or A = L*D*L**T as computed by DSPTRF, stored as */ 00134 /* a packed triangular matrix in the same storage format as A. */ 00135 00136 /* If FACT = 'N', then AFP is an output argument and on exit */ 00137 /* contains the block diagonal matrix D and the multipliers used */ 00138 /* to obtain the factor U or L from the factorization */ 00139 /* A = U*D*U**T or A = L*D*L**T as computed by DSPTRF, stored as */ 00140 /* a packed triangular matrix in the same storage format as A. */ 00141 00142 /* IPIV (input or output) INTEGER array, dimension (N) */ 00143 /* If FACT = 'F', then IPIV is an input argument and on entry */ 00144 /* contains details of the interchanges and the block structure */ 00145 /* of D, as determined by DSPTRF. */ 00146 /* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */ 00147 /* interchanged and D(k,k) is a 1-by-1 diagonal block. */ 00148 /* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */ 00149 /* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */ 00150 /* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */ 00151 /* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */ 00152 /* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */ 00153 00154 /* If FACT = 'N', then IPIV is an output argument and on exit */ 00155 /* contains details of the interchanges and the block structure */ 00156 /* of D, as determined by DSPTRF. */ 00157 00158 /* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00159 /* The N-by-NRHS right hand side matrix B. */ 00160 00161 /* LDB (input) INTEGER */ 00162 /* The leading dimension of the array B. LDB >= max(1,N). */ 00163 00164 /* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) */ 00165 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. */ 00166 00167 /* LDX (input) INTEGER */ 00168 /* The leading dimension of the array X. LDX >= max(1,N). */ 00169 00170 /* RCOND (output) DOUBLE PRECISION */ 00171 /* The estimate of the reciprocal condition number of the matrix */ 00172 /* A. If RCOND is less than the machine precision (in */ 00173 /* particular, if RCOND = 0), the matrix is singular to working */ 00174 /* precision. This condition is indicated by a return code of */ 00175 /* INFO > 0. */ 00176 00177 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00178 /* The estimated forward error bound for each solution vector */ 00179 /* X(j) (the j-th column of the solution matrix X). */ 00180 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00181 /* is an estimated upper bound for the magnitude of the largest */ 00182 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00183 /* largest element in X(j). The estimate is as reliable as */ 00184 /* the estimate for RCOND, and is almost always a slight */ 00185 /* overestimate of the true error. */ 00186 00187 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00188 /* The componentwise relative backward error of each solution */ 00189 /* vector X(j) (i.e., the smallest relative change in */ 00190 /* any element of A or B that makes X(j) an exact solution). */ 00191 00192 /* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */ 00193 00194 /* IWORK (workspace) INTEGER array, dimension (N) */ 00195 00196 /* INFO (output) INTEGER */ 00197 /* = 0: successful exit */ 00198 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00199 /* > 0: if INFO = i, and i is */ 00200 /* <= N: D(i,i) is exactly zero. The factorization */ 00201 /* has been completed but the factor D is exactly */ 00202 /* singular, so the solution and error bounds could */ 00203 /* not be computed. RCOND = 0 is returned. */ 00204 /* = N+1: D is nonsingular, but RCOND is less than machine */ 00205 /* precision, meaning that the matrix is singular */ 00206 /* to working precision. Nevertheless, the */ 00207 /* solution and error bounds are computed because */ 00208 /* there are a number of situations where the */ 00209 /* computed solution can be more accurate than the */ 00210 /* value of RCOND would suggest. */ 00211 00212 /* Further Details */ 00213 /* =============== */ 00214 00215 /* The packed storage scheme is illustrated by the following example */ 00216 /* when N = 4, UPLO = 'U': */ 00217 00218 /* Two-dimensional storage of the symmetric matrix A: */ 00219 00220 /* a11 a12 a13 a14 */ 00221 /* a22 a23 a24 */ 00222 /* a33 a34 (aij = aji) */ 00223 /* a44 */ 00224 00225 /* Packed storage of the upper triangle of A: */ 00226 00227 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */ 00228 00229 /* ===================================================================== */ 00230 00231 /* .. Parameters .. */ 00232 /* .. */ 00233 /* .. Local Scalars .. */ 00234 /* .. */ 00235 /* .. External Functions .. */ 00236 /* .. */ 00237 /* .. External Subroutines .. */ 00238 /* .. */ 00239 /* .. Intrinsic Functions .. */ 00240 /* .. */ 00241 /* .. Executable Statements .. */ 00242 00243 /* Test the input parameters. */ 00244 00245 /* Parameter adjustments */ 00246 --ap; 00247 --afp; 00248 --ipiv; 00249 b_dim1 = *ldb; 00250 b_offset = 1 + b_dim1; 00251 b -= b_offset; 00252 x_dim1 = *ldx; 00253 x_offset = 1 + x_dim1; 00254 x -= x_offset; 00255 --ferr; 00256 --berr; 00257 --work; 00258 --iwork; 00259 00260 /* Function Body */ 00261 *info = 0; 00262 nofact = lsame_(fact, "N"); 00263 if (! nofact && ! lsame_(fact, "F")) { 00264 *info = -1; 00265 } else if (! lsame_(uplo, "U") && ! lsame_(uplo, 00266 "L")) { 00267 *info = -2; 00268 } else if (*n < 0) { 00269 *info = -3; 00270 } else if (*nrhs < 0) { 00271 *info = -4; 00272 } else if (*ldb < max(1,*n)) { 00273 *info = -9; 00274 } else if (*ldx < max(1,*n)) { 00275 *info = -11; 00276 } 00277 if (*info != 0) { 00278 i__1 = -(*info); 00279 xerbla_("DSPSVX", &i__1); 00280 return 0; 00281 } 00282 00283 if (nofact) { 00284 00285 /* Compute the factorization A = U*D*U' or A = L*D*L'. */ 00286 00287 i__1 = *n * (*n + 1) / 2; 00288 dcopy_(&i__1, &ap[1], &c__1, &afp[1], &c__1); 00289 dsptrf_(uplo, n, &afp[1], &ipiv[1], info); 00290 00291 /* Return if INFO is non-zero. */ 00292 00293 if (*info > 0) { 00294 *rcond = 0.; 00295 return 0; 00296 } 00297 } 00298 00299 /* Compute the norm of the matrix A. */ 00300 00301 anorm = dlansp_("I", uplo, n, &ap[1], &work[1]); 00302 00303 /* Compute the reciprocal of the condition number of A. */ 00304 00305 dspcon_(uplo, n, &afp[1], &ipiv[1], &anorm, rcond, &work[1], &iwork[1], 00306 info); 00307 00308 /* Compute the solution vectors X. */ 00309 00310 dlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00311 dsptrs_(uplo, n, nrhs, &afp[1], &ipiv[1], &x[x_offset], ldx, info); 00312 00313 /* Use iterative refinement to improve the computed solutions and */ 00314 /* compute error bounds and backward error estimates for them. */ 00315 00316 dsprfs_(uplo, n, nrhs, &ap[1], &afp[1], &ipiv[1], &b[b_offset], ldb, &x[ 00317 x_offset], ldx, &ferr[1], &berr[1], &work[1], &iwork[1], info); 00318 00319 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00320 00321 if (*rcond < dlamch_("Epsilon")) { 00322 *info = *n + 1; 00323 } 00324 00325 return 0; 00326 00327 /* End of DSPSVX */ 00328 00329 } /* dspsvx_ */