00001 /* dsysvx.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 static integer c_n1 = -1; 00020 00021 /* Subroutine */ int dsysvx_(char *fact, char *uplo, integer *n, integer * 00022 nrhs, doublereal *a, integer *lda, doublereal *af, integer *ldaf, 00023 integer *ipiv, doublereal *b, integer *ldb, doublereal *x, integer * 00024 ldx, doublereal *rcond, doublereal *ferr, doublereal *berr, 00025 doublereal *work, integer *lwork, integer *iwork, integer *info) 00026 { 00027 /* System generated locals */ 00028 integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1, 00029 x_offset, i__1, i__2; 00030 00031 /* Local variables */ 00032 integer nb; 00033 extern logical lsame_(char *, char *); 00034 doublereal anorm; 00035 extern doublereal dlamch_(char *); 00036 logical nofact; 00037 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 00038 doublereal *, integer *, doublereal *, integer *), 00039 xerbla_(char *, integer *); 00040 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 00041 integer *, integer *); 00042 extern doublereal dlansy_(char *, char *, integer *, doublereal *, 00043 integer *, doublereal *); 00044 extern /* Subroutine */ int dsycon_(char *, integer *, doublereal *, 00045 integer *, integer *, doublereal *, doublereal *, doublereal *, 00046 integer *, integer *), dsyrfs_(char *, integer *, integer 00047 *, doublereal *, integer *, doublereal *, integer *, integer *, 00048 doublereal *, integer *, doublereal *, integer *, doublereal *, 00049 doublereal *, doublereal *, integer *, integer *), 00050 dsytrf_(char *, integer *, doublereal *, integer *, integer *, 00051 doublereal *, integer *, integer *); 00052 integer lwkopt; 00053 logical lquery; 00054 extern /* Subroutine */ int dsytrs_(char *, integer *, integer *, 00055 doublereal *, integer *, integer *, doublereal *, integer *, 00056 integer *); 00057 00058 00059 /* -- LAPACK driver routine (version 3.2) -- */ 00060 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00061 /* November 2006 */ 00062 00063 /* .. Scalar Arguments .. */ 00064 /* .. */ 00065 /* .. Array Arguments .. */ 00066 /* .. */ 00067 00068 /* Purpose */ 00069 /* ======= */ 00070 00071 /* DSYSVX uses the diagonal pivoting factorization to compute the */ 00072 /* solution to a real system of linear equations A * X = B, */ 00073 /* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS */ 00074 /* matrices. */ 00075 00076 /* Error bounds on the solution and a condition estimate are also */ 00077 /* provided. */ 00078 00079 /* Description */ 00080 /* =========== */ 00081 00082 /* The following steps are performed: */ 00083 00084 /* 1. If FACT = 'N', the diagonal pivoting method is used to factor A. */ 00085 /* The form of the factorization is */ 00086 /* A = U * D * U**T, if UPLO = 'U', or */ 00087 /* A = L * D * L**T, if UPLO = 'L', */ 00088 /* where U (or L) is a product of permutation and unit upper (lower) */ 00089 /* triangular matrices, and D is symmetric and block diagonal with */ 00090 /* 1-by-1 and 2-by-2 diagonal blocks. */ 00091 00092 /* 2. If some D(i,i)=0, so that D is exactly singular, then the routine */ 00093 /* returns with INFO = i. Otherwise, the factored form of A is used */ 00094 /* to estimate the condition number of the matrix A. If the */ 00095 /* reciprocal of the condition number is less than machine precision, */ 00096 /* INFO = N+1 is returned as a warning, but the routine still goes on */ 00097 /* to solve for X and compute error bounds as described below. */ 00098 00099 /* 3. The system of equations is solved for X using the factored form */ 00100 /* of A. */ 00101 00102 /* 4. Iterative refinement is applied to improve the computed solution */ 00103 /* matrix and calculate error bounds and backward error estimates */ 00104 /* for it. */ 00105 00106 /* Arguments */ 00107 /* ========= */ 00108 00109 /* FACT (input) CHARACTER*1 */ 00110 /* Specifies whether or not the factored form of A has been */ 00111 /* supplied on entry. */ 00112 /* = 'F': On entry, AF and IPIV contain the factored form of */ 00113 /* A. AF and IPIV will not be modified. */ 00114 /* = 'N': The matrix A will be copied to AF and factored. */ 00115 00116 /* UPLO (input) CHARACTER*1 */ 00117 /* = 'U': Upper triangle of A is stored; */ 00118 /* = 'L': Lower triangle of A is stored. */ 00119 00120 /* N (input) INTEGER */ 00121 /* The number of linear equations, i.e., the order of the */ 00122 /* matrix A. N >= 0. */ 00123 00124 /* NRHS (input) INTEGER */ 00125 /* The number of right hand sides, i.e., the number of columns */ 00126 /* of the matrices B and X. NRHS >= 0. */ 00127 00128 /* A (input) DOUBLE PRECISION array, dimension (LDA,N) */ 00129 /* The symmetric matrix A. If UPLO = 'U', the leading N-by-N */ 00130 /* upper triangular part of A contains the upper triangular part */ 00131 /* of the matrix A, and the strictly lower triangular part of A */ 00132 /* is not referenced. If UPLO = 'L', the leading N-by-N lower */ 00133 /* triangular part of A contains the lower triangular part of */ 00134 /* the matrix A, and the strictly upper triangular part of A is */ 00135 /* not referenced. */ 00136 00137 /* LDA (input) INTEGER */ 00138 /* The leading dimension of the array A. LDA >= max(1,N). */ 00139 00140 /* AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) */ 00141 /* If FACT = 'F', then AF is an input argument and on entry */ 00142 /* contains the block diagonal matrix D and the multipliers used */ 00143 /* to obtain the factor U or L from the factorization */ 00144 /* A = U*D*U**T or A = L*D*L**T as computed by DSYTRF. */ 00145 00146 /* If FACT = 'N', then AF is an output argument and on exit */ 00147 /* returns the block diagonal matrix D and the multipliers used */ 00148 /* to obtain the factor U or L from the factorization */ 00149 /* A = U*D*U**T or A = L*D*L**T. */ 00150 00151 /* LDAF (input) INTEGER */ 00152 /* The leading dimension of the array AF. LDAF >= max(1,N). */ 00153 00154 /* IPIV (input or output) INTEGER array, dimension (N) */ 00155 /* If FACT = 'F', then IPIV is an input argument and on entry */ 00156 /* contains details of the interchanges and the block structure */ 00157 /* of D, as determined by DSYTRF. */ 00158 /* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */ 00159 /* interchanged and D(k,k) is a 1-by-1 diagonal block. */ 00160 /* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */ 00161 /* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */ 00162 /* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */ 00163 /* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */ 00164 /* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */ 00165 00166 /* If FACT = 'N', then IPIV is an output argument and on exit */ 00167 /* contains details of the interchanges and the block structure */ 00168 /* of D, as determined by DSYTRF. */ 00169 00170 /* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00171 /* The N-by-NRHS right hand side matrix B. */ 00172 00173 /* LDB (input) INTEGER */ 00174 /* The leading dimension of the array B. LDB >= max(1,N). */ 00175 00176 /* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) */ 00177 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. */ 00178 00179 /* LDX (input) INTEGER */ 00180 /* The leading dimension of the array X. LDX >= max(1,N). */ 00181 00182 /* RCOND (output) DOUBLE PRECISION */ 00183 /* The estimate of the reciprocal condition number of the matrix */ 00184 /* A. If RCOND is less than the machine precision (in */ 00185 /* particular, if RCOND = 0), the matrix is singular to working */ 00186 /* precision. This condition is indicated by a return code of */ 00187 /* INFO > 0. */ 00188 00189 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00190 /* The estimated forward error bound for each solution vector */ 00191 /* X(j) (the j-th column of the solution matrix X). */ 00192 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00193 /* is an estimated upper bound for the magnitude of the largest */ 00194 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00195 /* largest element in X(j). The estimate is as reliable as */ 00196 /* the estimate for RCOND, and is almost always a slight */ 00197 /* overestimate of the true error. */ 00198 00199 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00200 /* The componentwise relative backward error of each solution */ 00201 /* vector X(j) (i.e., the smallest relative change in */ 00202 /* any element of A or B that makes X(j) an exact solution). */ 00203 00204 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ 00205 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ 00206 00207 /* LWORK (input) INTEGER */ 00208 /* The length of WORK. LWORK >= max(1,3*N), and for best */ 00209 /* performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where */ 00210 /* NB is the optimal blocksize for DSYTRF. */ 00211 00212 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00213 /* only calculates the optimal size of the WORK array, returns */ 00214 /* this value as the first entry of the WORK array, and no error */ 00215 /* message related to LWORK is issued by XERBLA. */ 00216 00217 /* IWORK (workspace) INTEGER array, dimension (N) */ 00218 00219 /* INFO (output) INTEGER */ 00220 /* = 0: successful exit */ 00221 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00222 /* > 0: if INFO = i, and i is */ 00223 /* <= N: D(i,i) is exactly zero. The factorization */ 00224 /* has been completed but the factor D is exactly */ 00225 /* singular, so the solution and error bounds could */ 00226 /* not be computed. RCOND = 0 is returned. */ 00227 /* = N+1: D is nonsingular, but RCOND is less than machine */ 00228 /* precision, meaning that the matrix is singular */ 00229 /* to working precision. Nevertheless, the */ 00230 /* solution and error bounds are computed because */ 00231 /* there are a number of situations where the */ 00232 /* computed solution can be more accurate than the */ 00233 /* value of RCOND would suggest. */ 00234 00235 /* ===================================================================== */ 00236 00237 /* .. Parameters .. */ 00238 /* .. */ 00239 /* .. Local Scalars .. */ 00240 /* .. */ 00241 /* .. External Functions .. */ 00242 /* .. */ 00243 /* .. External Subroutines .. */ 00244 /* .. */ 00245 /* .. Intrinsic Functions .. */ 00246 /* .. */ 00247 /* .. Executable Statements .. */ 00248 00249 /* Test the input parameters. */ 00250 00251 /* Parameter adjustments */ 00252 a_dim1 = *lda; 00253 a_offset = 1 + a_dim1; 00254 a -= a_offset; 00255 af_dim1 = *ldaf; 00256 af_offset = 1 + af_dim1; 00257 af -= af_offset; 00258 --ipiv; 00259 b_dim1 = *ldb; 00260 b_offset = 1 + b_dim1; 00261 b -= b_offset; 00262 x_dim1 = *ldx; 00263 x_offset = 1 + x_dim1; 00264 x -= x_offset; 00265 --ferr; 00266 --berr; 00267 --work; 00268 --iwork; 00269 00270 /* Function Body */ 00271 *info = 0; 00272 nofact = lsame_(fact, "N"); 00273 lquery = *lwork == -1; 00274 if (! nofact && ! lsame_(fact, "F")) { 00275 *info = -1; 00276 } else if (! lsame_(uplo, "U") && ! lsame_(uplo, 00277 "L")) { 00278 *info = -2; 00279 } else if (*n < 0) { 00280 *info = -3; 00281 } else if (*nrhs < 0) { 00282 *info = -4; 00283 } else if (*lda < max(1,*n)) { 00284 *info = -6; 00285 } else if (*ldaf < max(1,*n)) { 00286 *info = -8; 00287 } else if (*ldb < max(1,*n)) { 00288 *info = -11; 00289 } else if (*ldx < max(1,*n)) { 00290 *info = -13; 00291 } else /* if(complicated condition) */ { 00292 /* Computing MAX */ 00293 i__1 = 1, i__2 = *n * 3; 00294 if (*lwork < max(i__1,i__2) && ! lquery) { 00295 *info = -18; 00296 } 00297 } 00298 00299 if (*info == 0) { 00300 /* Computing MAX */ 00301 i__1 = 1, i__2 = *n * 3; 00302 lwkopt = max(i__1,i__2); 00303 if (nofact) { 00304 nb = ilaenv_(&c__1, "DSYTRF", uplo, n, &c_n1, &c_n1, &c_n1); 00305 /* Computing MAX */ 00306 i__1 = lwkopt, i__2 = *n * nb; 00307 lwkopt = max(i__1,i__2); 00308 } 00309 work[1] = (doublereal) lwkopt; 00310 } 00311 00312 if (*info != 0) { 00313 i__1 = -(*info); 00314 xerbla_("DSYSVX", &i__1); 00315 return 0; 00316 } else if (lquery) { 00317 return 0; 00318 } 00319 00320 if (nofact) { 00321 00322 /* Compute the factorization A = U*D*U' or A = L*D*L'. */ 00323 00324 dlacpy_(uplo, n, n, &a[a_offset], lda, &af[af_offset], ldaf); 00325 dsytrf_(uplo, n, &af[af_offset], ldaf, &ipiv[1], &work[1], lwork, 00326 info); 00327 00328 /* Return if INFO is non-zero. */ 00329 00330 if (*info > 0) { 00331 *rcond = 0.; 00332 return 0; 00333 } 00334 } 00335 00336 /* Compute the norm of the matrix A. */ 00337 00338 anorm = dlansy_("I", uplo, n, &a[a_offset], lda, &work[1]); 00339 00340 /* Compute the reciprocal of the condition number of A. */ 00341 00342 dsycon_(uplo, n, &af[af_offset], ldaf, &ipiv[1], &anorm, rcond, &work[1], 00343 &iwork[1], info); 00344 00345 /* Compute the solution vectors X. */ 00346 00347 dlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00348 dsytrs_(uplo, n, nrhs, &af[af_offset], ldaf, &ipiv[1], &x[x_offset], ldx, 00349 info); 00350 00351 /* Use iterative refinement to improve the computed solutions and */ 00352 /* compute error bounds and backward error estimates for them. */ 00353 00354 dsyrfs_(uplo, n, nrhs, &a[a_offset], lda, &af[af_offset], ldaf, &ipiv[1], 00355 &b[b_offset], ldb, &x[x_offset], ldx, &ferr[1], &berr[1], &work[1] 00356 , &iwork[1], info); 00357 00358 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00359 00360 if (*rcond < dlamch_("Epsilon")) { 00361 *info = *n + 1; 00362 } 00363 00364 work[1] = (doublereal) lwkopt; 00365 00366 return 0; 00367 00368 /* End of DSYSVX */ 00369 00370 } /* dsysvx_ */