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