00001 /* cpbsvx.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 cpbsvx_(char *fact, char *uplo, integer *n, integer *kd, 00021 integer *nrhs, complex *ab, integer *ldab, complex *afb, integer * 00022 ldafb, char *equed, real *s, complex *b, integer *ldb, complex *x, 00023 integer *ldx, real *rcond, real *ferr, real *berr, complex *work, 00024 real *rwork, 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, i__3, i__4, i__5; 00029 real r__1, r__2; 00030 complex q__1; 00031 00032 /* Local variables */ 00033 integer i__, j, j1, j2; 00034 real amax, smin, smax; 00035 extern logical lsame_(char *, char *); 00036 real scond, anorm; 00037 extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 00038 complex *, integer *); 00039 logical equil, rcequ, upper; 00040 extern doublereal clanhb_(char *, char *, integer *, integer *, complex *, 00041 integer *, real *); 00042 extern /* Subroutine */ int claqhb_(char *, integer *, integer *, complex 00043 *, integer *, real *, real *, real *, char *), 00044 cpbcon_(char *, integer *, integer *, complex *, integer *, real * 00045 , real *, complex *, real *, integer *); 00046 extern doublereal slamch_(char *); 00047 logical nofact; 00048 extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 00049 *, integer *, complex *, integer *), xerbla_(char *, 00050 integer *), cpbequ_(char *, integer *, integer *, complex 00051 *, integer *, real *, real *, real *, integer *), cpbrfs_( 00052 char *, integer *, integer *, integer *, complex *, integer *, 00053 complex *, integer *, complex *, integer *, complex *, integer *, 00054 real *, real *, complex *, real *, integer *); 00055 real bignum; 00056 extern /* Subroutine */ int cpbtrf_(char *, integer *, integer *, complex 00057 *, integer *, integer *); 00058 integer infequ; 00059 extern /* Subroutine */ int cpbtrs_(char *, integer *, integer *, integer 00060 *, complex *, integer *, complex *, integer *, integer *); 00061 real 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 /* CPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to */ 00077 /* compute the solution to a complex system of linear equations */ 00078 /* A * X = B, */ 00079 /* where A is an N-by-N Hermitian 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**H * U, if UPLO = 'U', or */ 00100 /* A = L * L**H, 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) COMPLEX array, dimension (LDAB,N) */ 00155 /* On entry, the upper or lower triangle of the Hermitian 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) COMPLEX 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**H*U or A = L*L**H 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**H*U or A = L*L**H. */ 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**H*U or A = L*L**H 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) REAL 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) COMPLEX 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) COMPLEX 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) REAL */ 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) REAL 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) REAL 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) COMPLEX array, dimension (2*N) */ 00244 00245 /* RWORK (workspace) REAL 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 Hermitian 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 --rwork; 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 = slamch_("Safe minimum"); 00336 bignum = 1.f / 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.f; 00362 i__1 = *n; 00363 for (j = 1; j <= i__1; ++j) { 00364 /* Computing MIN */ 00365 r__1 = smin, r__2 = s[j]; 00366 smin = dmin(r__1,r__2); 00367 /* Computing MAX */ 00368 r__1 = smax, r__2 = s[j]; 00369 smax = dmax(r__1,r__2); 00370 /* L10: */ 00371 } 00372 if (smin <= 0.f) { 00373 *info = -11; 00374 } else if (*n > 0) { 00375 scond = dmax(smin,smlnum) / dmin(smax,bignum); 00376 } else { 00377 scond = 1.f; 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_("CPBSVX", &i__1); 00392 return 0; 00393 } 00394 00395 if (equil) { 00396 00397 /* Compute row and column scalings to equilibrate the matrix A. */ 00398 00399 cpbequ_(uplo, n, kd, &ab[ab_offset], ldab, &s[1], &scond, &amax, & 00400 infequ); 00401 if (infequ == 0) { 00402 00403 /* Equilibrate the matrix. */ 00404 00405 claqhb_(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 i__3 = i__ + j * b_dim1; 00419 i__4 = i__; 00420 i__5 = i__ + j * b_dim1; 00421 q__1.r = s[i__4] * b[i__5].r, q__1.i = s[i__4] * b[i__5].i; 00422 b[i__3].r = q__1.r, b[i__3].i = q__1.i; 00423 /* L20: */ 00424 } 00425 /* L30: */ 00426 } 00427 } 00428 00429 if (nofact || equil) { 00430 00431 /* Compute the Cholesky factorization A = U'*U or A = L*L'. */ 00432 00433 if (upper) { 00434 i__1 = *n; 00435 for (j = 1; j <= i__1; ++j) { 00436 /* Computing MAX */ 00437 i__2 = j - *kd; 00438 j1 = max(i__2,1); 00439 i__2 = j - j1 + 1; 00440 ccopy_(&i__2, &ab[*kd + 1 - j + j1 + j * ab_dim1], &c__1, & 00441 afb[*kd + 1 - j + j1 + j * afb_dim1], &c__1); 00442 /* L40: */ 00443 } 00444 } else { 00445 i__1 = *n; 00446 for (j = 1; j <= i__1; ++j) { 00447 /* Computing MIN */ 00448 i__2 = j + *kd; 00449 j2 = min(i__2,*n); 00450 i__2 = j2 - j + 1; 00451 ccopy_(&i__2, &ab[j * ab_dim1 + 1], &c__1, &afb[j * afb_dim1 00452 + 1], &c__1); 00453 /* L50: */ 00454 } 00455 } 00456 00457 cpbtrf_(uplo, n, kd, &afb[afb_offset], ldafb, info); 00458 00459 /* Return if INFO is non-zero. */ 00460 00461 if (*info > 0) { 00462 *rcond = 0.f; 00463 return 0; 00464 } 00465 } 00466 00467 /* Compute the norm of the matrix A. */ 00468 00469 anorm = clanhb_("1", uplo, n, kd, &ab[ab_offset], ldab, &rwork[1]); 00470 00471 /* Compute the reciprocal of the condition number of A. */ 00472 00473 cpbcon_(uplo, n, kd, &afb[afb_offset], ldafb, &anorm, rcond, &work[1], & 00474 rwork[1], info); 00475 00476 /* Compute the solution matrix X. */ 00477 00478 clacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00479 cpbtrs_(uplo, n, kd, nrhs, &afb[afb_offset], ldafb, &x[x_offset], ldx, 00480 info); 00481 00482 /* Use iterative refinement to improve the computed solution and */ 00483 /* compute error bounds and backward error estimates for it. */ 00484 00485 cpbrfs_(uplo, n, kd, nrhs, &ab[ab_offset], ldab, &afb[afb_offset], ldafb, 00486 &b[b_offset], ldb, &x[x_offset], ldx, &ferr[1], &berr[1], &work[1] 00487 , &rwork[1], info); 00488 00489 /* Transform the solution matrix X to a solution of the original */ 00490 /* system. */ 00491 00492 if (rcequ) { 00493 i__1 = *nrhs; 00494 for (j = 1; j <= i__1; ++j) { 00495 i__2 = *n; 00496 for (i__ = 1; i__ <= i__2; ++i__) { 00497 i__3 = i__ + j * x_dim1; 00498 i__4 = i__; 00499 i__5 = i__ + j * x_dim1; 00500 q__1.r = s[i__4] * x[i__5].r, q__1.i = s[i__4] * x[i__5].i; 00501 x[i__3].r = q__1.r, x[i__3].i = q__1.i; 00502 /* L60: */ 00503 } 00504 /* L70: */ 00505 } 00506 i__1 = *nrhs; 00507 for (j = 1; j <= i__1; ++j) { 00508 ferr[j] /= scond; 00509 /* L80: */ 00510 } 00511 } 00512 00513 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00514 00515 if (*rcond < slamch_("Epsilon")) { 00516 *info = *n + 1; 00517 } 00518 00519 return 0; 00520 00521 /* End of CPBSVX */ 00522 00523 } /* cpbsvx_ */