00001 /* dpbsvx.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 dpbsvx_(char *fact, char *uplo, integer *n, integer *kd, 00021 integer *nrhs, doublereal *ab, integer *ldab, doublereal *afb, 00022 integer *ldafb, char *equed, doublereal *s, doublereal *b, integer * 00023 ldb, doublereal *x, integer *ldx, doublereal *rcond, doublereal *ferr, 00024 doublereal *berr, doublereal *work, integer *iwork, integer *info) 00025 { 00026 /* System generated locals */ 00027 integer ab_dim1, ab_offset, afb_dim1, afb_offset, b_dim1, b_offset, 00028 x_dim1, x_offset, i__1, i__2; 00029 doublereal d__1, d__2; 00030 00031 /* Local variables */ 00032 integer i__, j, j1, j2; 00033 doublereal amax, smin, smax; 00034 extern logical lsame_(char *, char *); 00035 doublereal scond, anorm; 00036 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 00037 doublereal *, integer *); 00038 logical equil, rcequ, upper; 00039 extern doublereal dlamch_(char *), dlansb_(char *, char *, 00040 integer *, integer *, doublereal *, integer *, doublereal *); 00041 extern /* Subroutine */ int dpbcon_(char *, integer *, integer *, 00042 doublereal *, integer *, doublereal *, doublereal *, doublereal *, 00043 integer *, integer *), dlaqsb_(char *, integer *, 00044 integer *, doublereal *, integer *, doublereal *, doublereal *, 00045 doublereal *, char *); 00046 logical nofact; 00047 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 00048 doublereal *, integer *, doublereal *, integer *), 00049 xerbla_(char *, integer *), dpbequ_(char *, integer *, 00050 integer *, doublereal *, integer *, doublereal *, doublereal *, 00051 doublereal *, integer *); 00052 doublereal bignum; 00053 extern /* Subroutine */ int dpbrfs_(char *, integer *, integer *, integer 00054 *, doublereal *, integer *, doublereal *, integer *, doublereal *, 00055 integer *, doublereal *, integer *, doublereal *, doublereal *, 00056 doublereal *, integer *, integer *), dpbtrf_(char *, 00057 integer *, integer *, doublereal *, integer *, integer *); 00058 integer infequ; 00059 extern /* Subroutine */ int dpbtrs_(char *, integer *, integer *, integer 00060 *, doublereal *, integer *, doublereal *, integer *, integer *); 00061 doublereal smlnum; 00062 00063 00064 /* -- LAPACK driver routine (version 3.2) -- */ 00065 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00066 /* November 2006 */ 00067 00068 /* .. Scalar Arguments .. */ 00069 /* .. */ 00070 /* .. Array Arguments .. */ 00071 /* .. */ 00072 00073 /* Purpose */ 00074 /* ======= */ 00075 00076 /* DPBSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to */ 00077 /* compute the solution to a real system of linear equations */ 00078 /* A * X = B, */ 00079 /* where A is an N-by-N symmetric positive definite band matrix and X */ 00080 /* and B are N-by-NRHS matrices. */ 00081 00082 /* Error bounds on the solution and a condition estimate are also */ 00083 /* provided. */ 00084 00085 /* Description */ 00086 /* =========== */ 00087 00088 /* The following steps are performed: */ 00089 00090 /* 1. If FACT = 'E', real scaling factors are computed to equilibrate */ 00091 /* the system: */ 00092 /* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B */ 00093 /* Whether or not the system will be equilibrated depends on the */ 00094 /* scaling of the matrix A, but if equilibration is used, A is */ 00095 /* overwritten by diag(S)*A*diag(S) and B by diag(S)*B. */ 00096 00097 /* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to */ 00098 /* factor the matrix A (after equilibration if FACT = 'E') as */ 00099 /* A = U**T * U, if UPLO = 'U', or */ 00100 /* A = L * L**T, if UPLO = 'L', */ 00101 /* where U is an upper triangular band matrix, and L is a lower */ 00102 /* triangular band matrix. */ 00103 00104 /* 3. If the leading i-by-i principal minor is not positive definite, */ 00105 /* then the routine returns with INFO = i. Otherwise, the factored */ 00106 /* form of A is used to estimate the condition number of the matrix */ 00107 /* A. If the reciprocal of the condition number is less than machine */ 00108 /* precision, INFO = N+1 is returned as a warning, but the routine */ 00109 /* still goes on to solve for X and compute error bounds as */ 00110 /* described below. */ 00111 00112 /* 4. The system of equations is solved for X using the factored form */ 00113 /* of A. */ 00114 00115 /* 5. Iterative refinement is applied to improve the computed solution */ 00116 /* matrix and calculate error bounds and backward error estimates */ 00117 /* for it. */ 00118 00119 /* 6. If equilibration was used, the matrix X is premultiplied by */ 00120 /* diag(S) so that it solves the original system before */ 00121 /* equilibration. */ 00122 00123 /* Arguments */ 00124 /* ========= */ 00125 00126 /* FACT (input) CHARACTER*1 */ 00127 /* Specifies whether or not the factored form of the matrix A is */ 00128 /* supplied on entry, and if not, whether the matrix A should be */ 00129 /* equilibrated before it is factored. */ 00130 /* = 'F': On entry, AFB contains the factored form of A. */ 00131 /* If EQUED = 'Y', the matrix A has been equilibrated */ 00132 /* with scaling factors given by S. AB and AFB will not */ 00133 /* be modified. */ 00134 /* = 'N': The matrix A will be copied to AFB and factored. */ 00135 /* = 'E': The matrix A will be equilibrated if necessary, then */ 00136 /* copied to AFB and factored. */ 00137 00138 /* UPLO (input) CHARACTER*1 */ 00139 /* = 'U': Upper triangle of A is stored; */ 00140 /* = 'L': Lower triangle of A is stored. */ 00141 00142 /* N (input) INTEGER */ 00143 /* The number of linear equations, i.e., the order of the */ 00144 /* matrix A. N >= 0. */ 00145 00146 /* KD (input) INTEGER */ 00147 /* The number of superdiagonals of the matrix A if UPLO = 'U', */ 00148 /* or the number of subdiagonals if UPLO = 'L'. KD >= 0. */ 00149 00150 /* NRHS (input) INTEGER */ 00151 /* The number of right-hand sides, i.e., the number of columns */ 00152 /* of the matrices B and X. NRHS >= 0. */ 00153 00154 /* AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) */ 00155 /* On entry, the upper or lower triangle of the symmetric band */ 00156 /* matrix A, stored in the first KD+1 rows of the array, except */ 00157 /* if FACT = 'F' and EQUED = 'Y', then A must contain the */ 00158 /* equilibrated matrix diag(S)*A*diag(S). The j-th column of A */ 00159 /* is stored in the j-th column of the array AB as follows: */ 00160 /* if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j; */ 00161 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(N,j+KD). */ 00162 /* See below for further details. */ 00163 00164 /* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by */ 00165 /* diag(S)*A*diag(S). */ 00166 00167 /* LDAB (input) INTEGER */ 00168 /* The leading dimension of the array A. LDAB >= KD+1. */ 00169 00170 /* AFB (input or output) DOUBLE PRECISION array, dimension (LDAFB,N) */ 00171 /* If FACT = 'F', then AFB is an input argument and on entry */ 00172 /* contains the triangular factor U or L from the Cholesky */ 00173 /* factorization A = U**T*U or A = L*L**T of the band matrix */ 00174 /* A, in the same storage format as A (see AB). If EQUED = 'Y', */ 00175 /* then AFB is the factored form of the equilibrated matrix A. */ 00176 00177 /* If FACT = 'N', then AFB is an output argument and on exit */ 00178 /* returns the triangular factor U or L from the Cholesky */ 00179 /* factorization A = U**T*U or A = L*L**T. */ 00180 00181 /* If FACT = 'E', then AFB is an output argument and on exit */ 00182 /* returns the triangular factor U or L from the Cholesky */ 00183 /* factorization A = U**T*U or A = L*L**T of the equilibrated */ 00184 /* matrix A (see the description of A for the form of the */ 00185 /* equilibrated matrix). */ 00186 00187 /* LDAFB (input) INTEGER */ 00188 /* The leading dimension of the array AFB. LDAFB >= KD+1. */ 00189 00190 /* EQUED (input or output) CHARACTER*1 */ 00191 /* Specifies the form of equilibration that was done. */ 00192 /* = 'N': No equilibration (always true if FACT = 'N'). */ 00193 /* = 'Y': Equilibration was done, i.e., A has been replaced by */ 00194 /* diag(S) * A * diag(S). */ 00195 /* EQUED is an input argument if FACT = 'F'; otherwise, it is an */ 00196 /* output argument. */ 00197 00198 /* S (input or output) DOUBLE PRECISION array, dimension (N) */ 00199 /* The scale factors for A; not accessed if EQUED = 'N'. S is */ 00200 /* an input argument if FACT = 'F'; otherwise, S is an output */ 00201 /* argument. If FACT = 'F' and EQUED = 'Y', each element of S */ 00202 /* must be positive. */ 00203 00204 /* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00205 /* On entry, the N-by-NRHS right hand side matrix B. */ 00206 /* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', */ 00207 /* B is overwritten by diag(S) * B. */ 00208 00209 /* LDB (input) INTEGER */ 00210 /* The leading dimension of the array B. LDB >= max(1,N). */ 00211 00212 /* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) */ 00213 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to */ 00214 /* the original system of equations. Note that if EQUED = 'Y', */ 00215 /* A and B are modified on exit, and the solution to the */ 00216 /* equilibrated system is inv(diag(S))*X. */ 00217 00218 /* LDX (input) INTEGER */ 00219 /* The leading dimension of the array X. LDX >= max(1,N). */ 00220 00221 /* RCOND (output) DOUBLE PRECISION */ 00222 /* The estimate of the reciprocal condition number of the matrix */ 00223 /* A after equilibration (if done). If RCOND is less than the */ 00224 /* machine precision (in particular, if RCOND = 0), the matrix */ 00225 /* is singular to working precision. This condition is */ 00226 /* indicated by a return code of INFO > 0. */ 00227 00228 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00229 /* The estimated forward error bound for each solution vector */ 00230 /* X(j) (the j-th column of the solution matrix X). */ 00231 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00232 /* is an estimated upper bound for the magnitude of the largest */ 00233 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00234 /* largest element in X(j). The estimate is as reliable as */ 00235 /* the estimate for RCOND, and is almost always a slight */ 00236 /* overestimate of the true error. */ 00237 00238 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00239 /* The componentwise relative backward error of each solution */ 00240 /* vector X(j) (i.e., the smallest relative change in */ 00241 /* any element of A or B that makes X(j) an exact solution). */ 00242 00243 /* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */ 00244 00245 /* IWORK (workspace) INTEGER array, dimension (N) */ 00246 00247 /* INFO (output) INTEGER */ 00248 /* = 0: successful exit */ 00249 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00250 /* > 0: if INFO = i, and i is */ 00251 /* <= N: the leading minor of order i of A is */ 00252 /* not positive definite, so the factorization */ 00253 /* could not be completed, and the solution has not */ 00254 /* been computed. RCOND = 0 is returned. */ 00255 /* = N+1: U is nonsingular, but RCOND is less than machine */ 00256 /* precision, meaning that the matrix is singular */ 00257 /* to working precision. Nevertheless, the */ 00258 /* solution and error bounds are computed because */ 00259 /* there are a number of situations where the */ 00260 /* computed solution can be more accurate than the */ 00261 /* value of RCOND would suggest. */ 00262 00263 /* Further Details */ 00264 /* =============== */ 00265 00266 /* The band storage scheme is illustrated by the following example, when */ 00267 /* N = 6, KD = 2, and UPLO = 'U': */ 00268 00269 /* Two-dimensional storage of the symmetric matrix A: */ 00270 00271 /* a11 a12 a13 */ 00272 /* a22 a23 a24 */ 00273 /* a33 a34 a35 */ 00274 /* a44 a45 a46 */ 00275 /* a55 a56 */ 00276 /* (aij=conjg(aji)) a66 */ 00277 00278 /* Band storage of the upper triangle of A: */ 00279 00280 /* * * a13 a24 a35 a46 */ 00281 /* * a12 a23 a34 a45 a56 */ 00282 /* a11 a22 a33 a44 a55 a66 */ 00283 00284 /* Similarly, if UPLO = 'L' the format of A is as follows: */ 00285 00286 /* a11 a22 a33 a44 a55 a66 */ 00287 /* a21 a32 a43 a54 a65 * */ 00288 /* a31 a42 a53 a64 * * */ 00289 00290 /* Array elements marked * are not used by the routine. */ 00291 00292 /* ===================================================================== */ 00293 00294 /* .. Parameters .. */ 00295 /* .. */ 00296 /* .. Local Scalars .. */ 00297 /* .. */ 00298 /* .. External Functions .. */ 00299 /* .. */ 00300 /* .. External Subroutines .. */ 00301 /* .. */ 00302 /* .. Intrinsic Functions .. */ 00303 /* .. */ 00304 /* .. Executable Statements .. */ 00305 00306 /* Parameter adjustments */ 00307 ab_dim1 = *ldab; 00308 ab_offset = 1 + ab_dim1; 00309 ab -= ab_offset; 00310 afb_dim1 = *ldafb; 00311 afb_offset = 1 + afb_dim1; 00312 afb -= afb_offset; 00313 --s; 00314 b_dim1 = *ldb; 00315 b_offset = 1 + b_dim1; 00316 b -= b_offset; 00317 x_dim1 = *ldx; 00318 x_offset = 1 + x_dim1; 00319 x -= x_offset; 00320 --ferr; 00321 --berr; 00322 --work; 00323 --iwork; 00324 00325 /* Function Body */ 00326 *info = 0; 00327 nofact = lsame_(fact, "N"); 00328 equil = lsame_(fact, "E"); 00329 upper = lsame_(uplo, "U"); 00330 if (nofact || equil) { 00331 *(unsigned char *)equed = 'N'; 00332 rcequ = FALSE_; 00333 } else { 00334 rcequ = lsame_(equed, "Y"); 00335 smlnum = dlamch_("Safe minimum"); 00336 bignum = 1. / smlnum; 00337 } 00338 00339 /* Test the input parameters. */ 00340 00341 if (! nofact && ! equil && ! lsame_(fact, "F")) { 00342 *info = -1; 00343 } else if (! upper && ! lsame_(uplo, "L")) { 00344 *info = -2; 00345 } else if (*n < 0) { 00346 *info = -3; 00347 } else if (*kd < 0) { 00348 *info = -4; 00349 } else if (*nrhs < 0) { 00350 *info = -5; 00351 } else if (*ldab < *kd + 1) { 00352 *info = -7; 00353 } else if (*ldafb < *kd + 1) { 00354 *info = -9; 00355 } else if (lsame_(fact, "F") && ! (rcequ || lsame_( 00356 equed, "N"))) { 00357 *info = -10; 00358 } else { 00359 if (rcequ) { 00360 smin = bignum; 00361 smax = 0.; 00362 i__1 = *n; 00363 for (j = 1; j <= i__1; ++j) { 00364 /* Computing MIN */ 00365 d__1 = smin, d__2 = s[j]; 00366 smin = min(d__1,d__2); 00367 /* Computing MAX */ 00368 d__1 = smax, d__2 = s[j]; 00369 smax = max(d__1,d__2); 00370 /* L10: */ 00371 } 00372 if (smin <= 0.) { 00373 *info = -11; 00374 } else if (*n > 0) { 00375 scond = max(smin,smlnum) / min(smax,bignum); 00376 } else { 00377 scond = 1.; 00378 } 00379 } 00380 if (*info == 0) { 00381 if (*ldb < max(1,*n)) { 00382 *info = -13; 00383 } else if (*ldx < max(1,*n)) { 00384 *info = -15; 00385 } 00386 } 00387 } 00388 00389 if (*info != 0) { 00390 i__1 = -(*info); 00391 xerbla_("DPBSVX", &i__1); 00392 return 0; 00393 } 00394 00395 if (equil) { 00396 00397 /* Compute row and column scalings to equilibrate the matrix A. */ 00398 00399 dpbequ_(uplo, n, kd, &ab[ab_offset], ldab, &s[1], &scond, &amax, & 00400 infequ); 00401 if (infequ == 0) { 00402 00403 /* Equilibrate the matrix. */ 00404 00405 dlaqsb_(uplo, n, kd, &ab[ab_offset], ldab, &s[1], &scond, &amax, 00406 equed); 00407 rcequ = lsame_(equed, "Y"); 00408 } 00409 } 00410 00411 /* Scale the right-hand side. */ 00412 00413 if (rcequ) { 00414 i__1 = *nrhs; 00415 for (j = 1; j <= i__1; ++j) { 00416 i__2 = *n; 00417 for (i__ = 1; i__ <= i__2; ++i__) { 00418 b[i__ + j * b_dim1] = s[i__] * b[i__ + j * b_dim1]; 00419 /* L20: */ 00420 } 00421 /* L30: */ 00422 } 00423 } 00424 00425 if (nofact || equil) { 00426 00427 /* Compute the Cholesky factorization A = U'*U or A = L*L'. */ 00428 00429 if (upper) { 00430 i__1 = *n; 00431 for (j = 1; j <= i__1; ++j) { 00432 /* Computing MAX */ 00433 i__2 = j - *kd; 00434 j1 = max(i__2,1); 00435 i__2 = j - j1 + 1; 00436 dcopy_(&i__2, &ab[*kd + 1 - j + j1 + j * ab_dim1], &c__1, & 00437 afb[*kd + 1 - j + j1 + j * afb_dim1], &c__1); 00438 /* L40: */ 00439 } 00440 } else { 00441 i__1 = *n; 00442 for (j = 1; j <= i__1; ++j) { 00443 /* Computing MIN */ 00444 i__2 = j + *kd; 00445 j2 = min(i__2,*n); 00446 i__2 = j2 - j + 1; 00447 dcopy_(&i__2, &ab[j * ab_dim1 + 1], &c__1, &afb[j * afb_dim1 00448 + 1], &c__1); 00449 /* L50: */ 00450 } 00451 } 00452 00453 dpbtrf_(uplo, n, kd, &afb[afb_offset], ldafb, info); 00454 00455 /* Return if INFO is non-zero. */ 00456 00457 if (*info > 0) { 00458 *rcond = 0.; 00459 return 0; 00460 } 00461 } 00462 00463 /* Compute the norm of the matrix A. */ 00464 00465 anorm = dlansb_("1", uplo, n, kd, &ab[ab_offset], ldab, &work[1]); 00466 00467 /* Compute the reciprocal of the condition number of A. */ 00468 00469 dpbcon_(uplo, n, kd, &afb[afb_offset], ldafb, &anorm, rcond, &work[1], & 00470 iwork[1], info); 00471 00472 /* Compute the solution matrix X. */ 00473 00474 dlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00475 dpbtrs_(uplo, n, kd, nrhs, &afb[afb_offset], ldafb, &x[x_offset], ldx, 00476 info); 00477 00478 /* Use iterative refinement to improve the computed solution and */ 00479 /* compute error bounds and backward error estimates for it. */ 00480 00481 dpbrfs_(uplo, n, kd, nrhs, &ab[ab_offset], ldab, &afb[afb_offset], ldafb, 00482 &b[b_offset], ldb, &x[x_offset], ldx, &ferr[1], &berr[1], &work[1] 00483 , &iwork[1], info); 00484 00485 /* Transform the solution matrix X to a solution of the original */ 00486 /* system. */ 00487 00488 if (rcequ) { 00489 i__1 = *nrhs; 00490 for (j = 1; j <= i__1; ++j) { 00491 i__2 = *n; 00492 for (i__ = 1; i__ <= i__2; ++i__) { 00493 x[i__ + j * x_dim1] = s[i__] * x[i__ + j * x_dim1]; 00494 /* L60: */ 00495 } 00496 /* L70: */ 00497 } 00498 i__1 = *nrhs; 00499 for (j = 1; j <= i__1; ++j) { 00500 ferr[j] /= scond; 00501 /* L80: */ 00502 } 00503 } 00504 00505 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00506 00507 if (*rcond < dlamch_("Epsilon")) { 00508 *info = *n + 1; 00509 } 00510 00511 return 0; 00512 00513 /* End of DPBSVX */ 00514 00515 } /* dpbsvx_ */