00001 /* cgbsvx.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 cgbsvx_(char *fact, char *trans, integer *n, integer *kl, 00021 integer *ku, integer *nrhs, complex *ab, integer *ldab, complex *afb, 00022 integer *ldafb, integer *ipiv, char *equed, real *r__, real *c__, 00023 complex *b, integer *ldb, complex *x, integer *ldx, real *rcond, real 00024 *ferr, real *berr, complex *work, 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 /* Builtin functions */ 00033 double c_abs(complex *); 00034 00035 /* Local variables */ 00036 integer i__, j, j1, j2; 00037 real amax; 00038 char norm[1]; 00039 extern logical lsame_(char *, char *); 00040 real rcmin, rcmax, anorm; 00041 extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 00042 complex *, integer *); 00043 logical equil; 00044 extern doublereal clangb_(char *, integer *, integer *, integer *, 00045 complex *, integer *, real *); 00046 extern /* Subroutine */ int claqgb_(integer *, integer *, integer *, 00047 integer *, complex *, integer *, real *, real *, real *, real *, 00048 real *, char *), cgbcon_(char *, integer *, integer *, 00049 integer *, complex *, integer *, integer *, real *, real *, 00050 complex *, real *, integer *); 00051 real colcnd; 00052 extern doublereal clantb_(char *, char *, char *, integer *, integer *, 00053 complex *, integer *, real *); 00054 extern /* Subroutine */ int cgbequ_(integer *, integer *, integer *, 00055 integer *, complex *, integer *, real *, real *, real *, real *, 00056 real *, integer *); 00057 extern doublereal slamch_(char *); 00058 extern /* Subroutine */ int cgbrfs_(char *, integer *, integer *, integer 00059 *, integer *, complex *, integer *, complex *, integer *, integer 00060 *, complex *, integer *, complex *, integer *, real *, real *, 00061 complex *, real *, integer *), cgbtrf_(integer *, integer 00062 *, integer *, integer *, complex *, integer *, integer *, integer 00063 *); 00064 logical nofact; 00065 extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 00066 *, integer *, complex *, integer *), xerbla_(char *, 00067 integer *); 00068 real bignum; 00069 extern /* Subroutine */ int cgbtrs_(char *, integer *, integer *, integer 00070 *, integer *, complex *, integer *, integer *, complex *, integer 00071 *, integer *); 00072 integer infequ; 00073 logical colequ; 00074 real rowcnd; 00075 logical notran; 00076 real smlnum; 00077 logical rowequ; 00078 real rpvgrw; 00079 00080 00081 /* -- LAPACK driver routine (version 3.2) -- */ 00082 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00083 /* November 2006 */ 00084 00085 /* .. Scalar Arguments .. */ 00086 /* .. */ 00087 /* .. Array Arguments .. */ 00088 /* .. */ 00089 00090 /* Purpose */ 00091 /* ======= */ 00092 00093 /* CGBSVX uses the LU factorization to compute the solution to a complex */ 00094 /* system of linear equations A * X = B, A**T * X = B, or A**H * X = B, */ 00095 /* where A is a band matrix of order N with KL subdiagonals and KU */ 00096 /* superdiagonals, and X and B are N-by-NRHS matrices. */ 00097 00098 /* Error bounds on the solution and a condition estimate are also */ 00099 /* provided. */ 00100 00101 /* Description */ 00102 /* =========== */ 00103 00104 /* The following steps are performed by this subroutine: */ 00105 00106 /* 1. If FACT = 'E', real scaling factors are computed to equilibrate */ 00107 /* the system: */ 00108 /* TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B */ 00109 /* TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B */ 00110 /* TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B */ 00111 /* Whether or not the system will be equilibrated depends on the */ 00112 /* scaling of the matrix A, but if equilibration is used, A is */ 00113 /* overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N') */ 00114 /* or diag(C)*B (if TRANS = 'T' or 'C'). */ 00115 00116 /* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the */ 00117 /* matrix A (after equilibration if FACT = 'E') as */ 00118 /* A = L * U, */ 00119 /* where L is a product of permutation and unit lower triangular */ 00120 /* matrices with KL subdiagonals, and U is upper triangular with */ 00121 /* KL+KU superdiagonals. */ 00122 00123 /* 3. If some U(i,i)=0, so that U is exactly singular, then the routine */ 00124 /* returns with INFO = i. Otherwise, the factored form of A is used */ 00125 /* to estimate the condition number of the matrix A. If the */ 00126 /* reciprocal of the condition number is less than machine precision, */ 00127 /* INFO = N+1 is returned as a warning, but the routine still goes on */ 00128 /* to solve for X and compute error bounds as described below. */ 00129 00130 /* 4. The system of equations is solved for X using the factored form */ 00131 /* of A. */ 00132 00133 /* 5. Iterative refinement is applied to improve the computed solution */ 00134 /* matrix and calculate error bounds and backward error estimates */ 00135 /* for it. */ 00136 00137 /* 6. If equilibration was used, the matrix X is premultiplied by */ 00138 /* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so */ 00139 /* that it solves the original system before equilibration. */ 00140 00141 /* Arguments */ 00142 /* ========= */ 00143 00144 /* FACT (input) CHARACTER*1 */ 00145 /* Specifies whether or not the factored form of the matrix A is */ 00146 /* supplied on entry, and if not, whether the matrix A should be */ 00147 /* equilibrated before it is factored. */ 00148 /* = 'F': On entry, AFB and IPIV contain the factored form of */ 00149 /* A. If EQUED is not 'N', the matrix A has been */ 00150 /* equilibrated with scaling factors given by R and C. */ 00151 /* AB, AFB, and IPIV are not modified. */ 00152 /* = 'N': The matrix A will be copied to AFB and factored. */ 00153 /* = 'E': The matrix A will be equilibrated if necessary, then */ 00154 /* copied to AFB and factored. */ 00155 00156 /* TRANS (input) CHARACTER*1 */ 00157 /* Specifies the form of the system of equations. */ 00158 /* = 'N': A * X = B (No transpose) */ 00159 /* = 'T': A**T * X = B (Transpose) */ 00160 /* = 'C': A**H * X = B (Conjugate transpose) */ 00161 00162 /* N (input) INTEGER */ 00163 /* The number of linear equations, i.e., the order of the */ 00164 /* matrix A. N >= 0. */ 00165 00166 /* KL (input) INTEGER */ 00167 /* The number of subdiagonals within the band of A. KL >= 0. */ 00168 00169 /* KU (input) INTEGER */ 00170 /* The number of superdiagonals within the band of A. KU >= 0. */ 00171 00172 /* NRHS (input) INTEGER */ 00173 /* The number of right hand sides, i.e., the number of columns */ 00174 /* of the matrices B and X. NRHS >= 0. */ 00175 00176 /* AB (input/output) COMPLEX array, dimension (LDAB,N) */ 00177 /* On entry, the matrix A in band storage, in rows 1 to KL+KU+1. */ 00178 /* The j-th column of A is stored in the j-th column of the */ 00179 /* array AB as follows: */ 00180 /* AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) */ 00181 00182 /* If FACT = 'F' and EQUED is not 'N', then A must have been */ 00183 /* equilibrated by the scaling factors in R and/or C. AB is not */ 00184 /* modified if FACT = 'F' or 'N', or if FACT = 'E' and */ 00185 /* EQUED = 'N' on exit. */ 00186 00187 /* On exit, if EQUED .ne. 'N', A is scaled as follows: */ 00188 /* EQUED = 'R': A := diag(R) * A */ 00189 /* EQUED = 'C': A := A * diag(C) */ 00190 /* EQUED = 'B': A := diag(R) * A * diag(C). */ 00191 00192 /* LDAB (input) INTEGER */ 00193 /* The leading dimension of the array AB. LDAB >= KL+KU+1. */ 00194 00195 /* AFB (input or output) COMPLEX array, dimension (LDAFB,N) */ 00196 /* If FACT = 'F', then AFB is an input argument and on entry */ 00197 /* contains details of the LU factorization of the band matrix */ 00198 /* A, as computed by CGBTRF. U is stored as an upper triangular */ 00199 /* band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, */ 00200 /* and the multipliers used during the factorization are stored */ 00201 /* in rows KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is */ 00202 /* the factored form of the equilibrated matrix A. */ 00203 00204 /* If FACT = 'N', then AFB is an output argument and on exit */ 00205 /* returns details of the LU factorization of A. */ 00206 00207 /* If FACT = 'E', then AFB is an output argument and on exit */ 00208 /* returns details of the LU factorization of the equilibrated */ 00209 /* matrix A (see the description of AB for the form of the */ 00210 /* equilibrated matrix). */ 00211 00212 /* LDAFB (input) INTEGER */ 00213 /* The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. */ 00214 00215 /* IPIV (input or output) INTEGER array, dimension (N) */ 00216 /* If FACT = 'F', then IPIV is an input argument and on entry */ 00217 /* contains the pivot indices from the factorization A = L*U */ 00218 /* as computed by CGBTRF; row i of the matrix was interchanged */ 00219 /* with row IPIV(i). */ 00220 00221 /* If FACT = 'N', then IPIV is an output argument and on exit */ 00222 /* contains the pivot indices from the factorization A = L*U */ 00223 /* of the original matrix A. */ 00224 00225 /* If FACT = 'E', then IPIV is an output argument and on exit */ 00226 /* contains the pivot indices from the factorization A = L*U */ 00227 /* of the equilibrated matrix A. */ 00228 00229 /* EQUED (input or output) CHARACTER*1 */ 00230 /* Specifies the form of equilibration that was done. */ 00231 /* = 'N': No equilibration (always true if FACT = 'N'). */ 00232 /* = 'R': Row equilibration, i.e., A has been premultiplied by */ 00233 /* diag(R). */ 00234 /* = 'C': Column equilibration, i.e., A has been postmultiplied */ 00235 /* by diag(C). */ 00236 /* = 'B': Both row and column equilibration, i.e., A has been */ 00237 /* replaced by diag(R) * A * diag(C). */ 00238 /* EQUED is an input argument if FACT = 'F'; otherwise, it is an */ 00239 /* output argument. */ 00240 00241 /* R (input or output) REAL array, dimension (N) */ 00242 /* The row scale factors for A. If EQUED = 'R' or 'B', A is */ 00243 /* multiplied on the left by diag(R); if EQUED = 'N' or 'C', R */ 00244 /* is not accessed. R is an input argument if FACT = 'F'; */ 00245 /* otherwise, R is an output argument. If FACT = 'F' and */ 00246 /* EQUED = 'R' or 'B', each element of R must be positive. */ 00247 00248 /* C (input or output) REAL array, dimension (N) */ 00249 /* The column scale factors for A. If EQUED = 'C' or 'B', A is */ 00250 /* multiplied on the right by diag(C); if EQUED = 'N' or 'R', C */ 00251 /* is not accessed. C is an input argument if FACT = 'F'; */ 00252 /* otherwise, C is an output argument. If FACT = 'F' and */ 00253 /* EQUED = 'C' or 'B', each element of C must be positive. */ 00254 00255 /* B (input/output) COMPLEX array, dimension (LDB,NRHS) */ 00256 /* On entry, the right hand side matrix B. */ 00257 /* On exit, */ 00258 /* if EQUED = 'N', B is not modified; */ 00259 /* if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by */ 00260 /* diag(R)*B; */ 00261 /* if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is */ 00262 /* overwritten by diag(C)*B. */ 00263 00264 /* LDB (input) INTEGER */ 00265 /* The leading dimension of the array B. LDB >= max(1,N). */ 00266 00267 /* X (output) COMPLEX array, dimension (LDX,NRHS) */ 00268 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X */ 00269 /* to the original system of equations. Note that A and B are */ 00270 /* modified on exit if EQUED .ne. 'N', and the solution to the */ 00271 /* equilibrated system is inv(diag(C))*X if TRANS = 'N' and */ 00272 /* EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' */ 00273 /* and EQUED = 'R' or 'B'. */ 00274 00275 /* LDX (input) INTEGER */ 00276 /* The leading dimension of the array X. LDX >= max(1,N). */ 00277 00278 /* RCOND (output) REAL */ 00279 /* The estimate of the reciprocal condition number of the matrix */ 00280 /* A after equilibration (if done). If RCOND is less than the */ 00281 /* machine precision (in particular, if RCOND = 0), the matrix */ 00282 /* is singular to working precision. This condition is */ 00283 /* indicated by a return code of INFO > 0. */ 00284 00285 /* FERR (output) REAL array, dimension (NRHS) */ 00286 /* The estimated forward error bound for each solution vector */ 00287 /* X(j) (the j-th column of the solution matrix X). */ 00288 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00289 /* is an estimated upper bound for the magnitude of the largest */ 00290 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00291 /* largest element in X(j). The estimate is as reliable as */ 00292 /* the estimate for RCOND, and is almost always a slight */ 00293 /* overestimate of the true error. */ 00294 00295 /* BERR (output) REAL array, dimension (NRHS) */ 00296 /* The componentwise relative backward error of each solution */ 00297 /* vector X(j) (i.e., the smallest relative change in */ 00298 /* any element of A or B that makes X(j) an exact solution). */ 00299 00300 /* WORK (workspace) COMPLEX array, dimension (2*N) */ 00301 00302 /* RWORK (workspace/output) REAL array, dimension (N) */ 00303 /* On exit, RWORK(1) contains the reciprocal pivot growth */ 00304 /* factor norm(A)/norm(U). The "max absolute element" norm is */ 00305 /* used. If RWORK(1) is much less than 1, then the stability */ 00306 /* of the LU factorization of the (equilibrated) matrix A */ 00307 /* could be poor. This also means that the solution X, condition */ 00308 /* estimator RCOND, and forward error bound FERR could be */ 00309 /* unreliable. If factorization fails with 0<INFO<=N, then */ 00310 /* RWORK(1) contains the reciprocal pivot growth factor for the */ 00311 /* leading INFO columns of A. */ 00312 00313 /* INFO (output) INTEGER */ 00314 /* = 0: successful exit */ 00315 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00316 /* > 0: if INFO = i, and i is */ 00317 /* <= N: U(i,i) is exactly zero. The factorization */ 00318 /* has been completed, but the factor U is exactly */ 00319 /* singular, so the solution and error bounds */ 00320 /* could not be computed. RCOND = 0 is returned. */ 00321 /* = N+1: U is nonsingular, but RCOND is less than machine */ 00322 /* precision, meaning that the matrix is singular */ 00323 /* to working precision. Nevertheless, the */ 00324 /* solution and error bounds are computed because */ 00325 /* there are a number of situations where the */ 00326 /* computed solution can be more accurate than the */ 00327 /* value of RCOND would suggest. */ 00328 00329 /* ===================================================================== */ 00330 /* Moved setting of INFO = N+1 so INFO does not subsequently get */ 00331 /* overwritten. Sven, 17 Mar 05. */ 00332 /* ===================================================================== */ 00333 00334 /* .. Parameters .. */ 00335 /* .. */ 00336 /* .. Local Scalars .. */ 00337 /* .. */ 00338 /* .. External Functions .. */ 00339 /* .. */ 00340 /* .. External Subroutines .. */ 00341 /* .. */ 00342 /* .. Intrinsic Functions .. */ 00343 /* .. */ 00344 /* .. Executable Statements .. */ 00345 00346 /* Parameter adjustments */ 00347 ab_dim1 = *ldab; 00348 ab_offset = 1 + ab_dim1; 00349 ab -= ab_offset; 00350 afb_dim1 = *ldafb; 00351 afb_offset = 1 + afb_dim1; 00352 afb -= afb_offset; 00353 --ipiv; 00354 --r__; 00355 --c__; 00356 b_dim1 = *ldb; 00357 b_offset = 1 + b_dim1; 00358 b -= b_offset; 00359 x_dim1 = *ldx; 00360 x_offset = 1 + x_dim1; 00361 x -= x_offset; 00362 --ferr; 00363 --berr; 00364 --work; 00365 --rwork; 00366 00367 /* Function Body */ 00368 *info = 0; 00369 nofact = lsame_(fact, "N"); 00370 equil = lsame_(fact, "E"); 00371 notran = lsame_(trans, "N"); 00372 if (nofact || equil) { 00373 *(unsigned char *)equed = 'N'; 00374 rowequ = FALSE_; 00375 colequ = FALSE_; 00376 } else { 00377 rowequ = lsame_(equed, "R") || lsame_(equed, 00378 "B"); 00379 colequ = lsame_(equed, "C") || lsame_(equed, 00380 "B"); 00381 smlnum = slamch_("Safe minimum"); 00382 bignum = 1.f / smlnum; 00383 } 00384 00385 /* Test the input parameters. */ 00386 00387 if (! nofact && ! equil && ! lsame_(fact, "F")) { 00388 *info = -1; 00389 } else if (! notran && ! lsame_(trans, "T") && ! 00390 lsame_(trans, "C")) { 00391 *info = -2; 00392 } else if (*n < 0) { 00393 *info = -3; 00394 } else if (*kl < 0) { 00395 *info = -4; 00396 } else if (*ku < 0) { 00397 *info = -5; 00398 } else if (*nrhs < 0) { 00399 *info = -6; 00400 } else if (*ldab < *kl + *ku + 1) { 00401 *info = -8; 00402 } else if (*ldafb < (*kl << 1) + *ku + 1) { 00403 *info = -10; 00404 } else if (lsame_(fact, "F") && ! (rowequ || colequ 00405 || lsame_(equed, "N"))) { 00406 *info = -12; 00407 } else { 00408 if (rowequ) { 00409 rcmin = bignum; 00410 rcmax = 0.f; 00411 i__1 = *n; 00412 for (j = 1; j <= i__1; ++j) { 00413 /* Computing MIN */ 00414 r__1 = rcmin, r__2 = r__[j]; 00415 rcmin = dmin(r__1,r__2); 00416 /* Computing MAX */ 00417 r__1 = rcmax, r__2 = r__[j]; 00418 rcmax = dmax(r__1,r__2); 00419 /* L10: */ 00420 } 00421 if (rcmin <= 0.f) { 00422 *info = -13; 00423 } else if (*n > 0) { 00424 rowcnd = dmax(rcmin,smlnum) / dmin(rcmax,bignum); 00425 } else { 00426 rowcnd = 1.f; 00427 } 00428 } 00429 if (colequ && *info == 0) { 00430 rcmin = bignum; 00431 rcmax = 0.f; 00432 i__1 = *n; 00433 for (j = 1; j <= i__1; ++j) { 00434 /* Computing MIN */ 00435 r__1 = rcmin, r__2 = c__[j]; 00436 rcmin = dmin(r__1,r__2); 00437 /* Computing MAX */ 00438 r__1 = rcmax, r__2 = c__[j]; 00439 rcmax = dmax(r__1,r__2); 00440 /* L20: */ 00441 } 00442 if (rcmin <= 0.f) { 00443 *info = -14; 00444 } else if (*n > 0) { 00445 colcnd = dmax(rcmin,smlnum) / dmin(rcmax,bignum); 00446 } else { 00447 colcnd = 1.f; 00448 } 00449 } 00450 if (*info == 0) { 00451 if (*ldb < max(1,*n)) { 00452 *info = -16; 00453 } else if (*ldx < max(1,*n)) { 00454 *info = -18; 00455 } 00456 } 00457 } 00458 00459 if (*info != 0) { 00460 i__1 = -(*info); 00461 xerbla_("CGBSVX", &i__1); 00462 return 0; 00463 } 00464 00465 if (equil) { 00466 00467 /* Compute row and column scalings to equilibrate the matrix A. */ 00468 00469 cgbequ_(n, n, kl, ku, &ab[ab_offset], ldab, &r__[1], &c__[1], &rowcnd, 00470 &colcnd, &amax, &infequ); 00471 if (infequ == 0) { 00472 00473 /* Equilibrate the matrix. */ 00474 00475 claqgb_(n, n, kl, ku, &ab[ab_offset], ldab, &r__[1], &c__[1], & 00476 rowcnd, &colcnd, &amax, equed); 00477 rowequ = lsame_(equed, "R") || lsame_(equed, 00478 "B"); 00479 colequ = lsame_(equed, "C") || lsame_(equed, 00480 "B"); 00481 } 00482 } 00483 00484 /* Scale the right hand side. */ 00485 00486 if (notran) { 00487 if (rowequ) { 00488 i__1 = *nrhs; 00489 for (j = 1; j <= i__1; ++j) { 00490 i__2 = *n; 00491 for (i__ = 1; i__ <= i__2; ++i__) { 00492 i__3 = i__ + j * b_dim1; 00493 i__4 = i__; 00494 i__5 = i__ + j * b_dim1; 00495 q__1.r = r__[i__4] * b[i__5].r, q__1.i = r__[i__4] * b[ 00496 i__5].i; 00497 b[i__3].r = q__1.r, b[i__3].i = q__1.i; 00498 /* L30: */ 00499 } 00500 /* L40: */ 00501 } 00502 } 00503 } else if (colequ) { 00504 i__1 = *nrhs; 00505 for (j = 1; j <= i__1; ++j) { 00506 i__2 = *n; 00507 for (i__ = 1; i__ <= i__2; ++i__) { 00508 i__3 = i__ + j * b_dim1; 00509 i__4 = i__; 00510 i__5 = i__ + j * b_dim1; 00511 q__1.r = c__[i__4] * b[i__5].r, q__1.i = c__[i__4] * b[i__5] 00512 .i; 00513 b[i__3].r = q__1.r, b[i__3].i = q__1.i; 00514 /* L50: */ 00515 } 00516 /* L60: */ 00517 } 00518 } 00519 00520 if (nofact || equil) { 00521 00522 /* Compute the LU factorization of the band matrix A. */ 00523 00524 i__1 = *n; 00525 for (j = 1; j <= i__1; ++j) { 00526 /* Computing MAX */ 00527 i__2 = j - *ku; 00528 j1 = max(i__2,1); 00529 /* Computing MIN */ 00530 i__2 = j + *kl; 00531 j2 = min(i__2,*n); 00532 i__2 = j2 - j1 + 1; 00533 ccopy_(&i__2, &ab[*ku + 1 - j + j1 + j * ab_dim1], &c__1, &afb[* 00534 kl + *ku + 1 - j + j1 + j * afb_dim1], &c__1); 00535 /* L70: */ 00536 } 00537 00538 cgbtrf_(n, n, kl, ku, &afb[afb_offset], ldafb, &ipiv[1], info); 00539 00540 /* Return if INFO is non-zero. */ 00541 00542 if (*info > 0) { 00543 00544 /* Compute the reciprocal pivot growth factor of the */ 00545 /* leading rank-deficient INFO columns of A. */ 00546 00547 anorm = 0.f; 00548 i__1 = *info; 00549 for (j = 1; j <= i__1; ++j) { 00550 /* Computing MAX */ 00551 i__2 = *ku + 2 - j; 00552 /* Computing MIN */ 00553 i__4 = *n + *ku + 1 - j, i__5 = *kl + *ku + 1; 00554 i__3 = min(i__4,i__5); 00555 for (i__ = max(i__2,1); i__ <= i__3; ++i__) { 00556 /* Computing MAX */ 00557 r__1 = anorm, r__2 = c_abs(&ab[i__ + j * ab_dim1]); 00558 anorm = dmax(r__1,r__2); 00559 /* L80: */ 00560 } 00561 /* L90: */ 00562 } 00563 /* Computing MIN */ 00564 i__3 = *info - 1, i__2 = *kl + *ku; 00565 i__1 = min(i__3,i__2); 00566 /* Computing MAX */ 00567 i__4 = 1, i__5 = *kl + *ku + 2 - *info; 00568 rpvgrw = clantb_("M", "U", "N", info, &i__1, &afb[max(i__4, i__5) 00569 + afb_dim1], ldafb, &rwork[1]); 00570 if (rpvgrw == 0.f) { 00571 rpvgrw = 1.f; 00572 } else { 00573 rpvgrw = anorm / rpvgrw; 00574 } 00575 rwork[1] = rpvgrw; 00576 *rcond = 0.f; 00577 return 0; 00578 } 00579 } 00580 00581 /* Compute the norm of the matrix A and the */ 00582 /* reciprocal pivot growth factor RPVGRW. */ 00583 00584 if (notran) { 00585 *(unsigned char *)norm = '1'; 00586 } else { 00587 *(unsigned char *)norm = 'I'; 00588 } 00589 anorm = clangb_(norm, n, kl, ku, &ab[ab_offset], ldab, &rwork[1]); 00590 i__1 = *kl + *ku; 00591 rpvgrw = clantb_("M", "U", "N", n, &i__1, &afb[afb_offset], ldafb, &rwork[ 00592 1]); 00593 if (rpvgrw == 0.f) { 00594 rpvgrw = 1.f; 00595 } else { 00596 rpvgrw = clangb_("M", n, kl, ku, &ab[ab_offset], ldab, &rwork[1]) / rpvgrw; 00597 } 00598 00599 /* Compute the reciprocal of the condition number of A. */ 00600 00601 cgbcon_(norm, n, kl, ku, &afb[afb_offset], ldafb, &ipiv[1], &anorm, rcond, 00602 &work[1], &rwork[1], info); 00603 00604 /* Compute the solution matrix X. */ 00605 00606 clacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00607 cgbtrs_(trans, n, kl, ku, nrhs, &afb[afb_offset], ldafb, &ipiv[1], &x[ 00608 x_offset], ldx, info); 00609 00610 /* Use iterative refinement to improve the computed solution and */ 00611 /* compute error bounds and backward error estimates for it. */ 00612 00613 cgbrfs_(trans, n, kl, ku, nrhs, &ab[ab_offset], ldab, &afb[afb_offset], 00614 ldafb, &ipiv[1], &b[b_offset], ldb, &x[x_offset], ldx, &ferr[1], & 00615 berr[1], &work[1], &rwork[1], info); 00616 00617 /* Transform the solution matrix X to a solution of the original */ 00618 /* system. */ 00619 00620 if (notran) { 00621 if (colequ) { 00622 i__1 = *nrhs; 00623 for (j = 1; j <= i__1; ++j) { 00624 i__3 = *n; 00625 for (i__ = 1; i__ <= i__3; ++i__) { 00626 i__2 = i__ + j * x_dim1; 00627 i__4 = i__; 00628 i__5 = i__ + j * x_dim1; 00629 q__1.r = c__[i__4] * x[i__5].r, q__1.i = c__[i__4] * x[ 00630 i__5].i; 00631 x[i__2].r = q__1.r, x[i__2].i = q__1.i; 00632 /* L100: */ 00633 } 00634 /* L110: */ 00635 } 00636 i__1 = *nrhs; 00637 for (j = 1; j <= i__1; ++j) { 00638 ferr[j] /= colcnd; 00639 /* L120: */ 00640 } 00641 } 00642 } else if (rowequ) { 00643 i__1 = *nrhs; 00644 for (j = 1; j <= i__1; ++j) { 00645 i__3 = *n; 00646 for (i__ = 1; i__ <= i__3; ++i__) { 00647 i__2 = i__ + j * x_dim1; 00648 i__4 = i__; 00649 i__5 = i__ + j * x_dim1; 00650 q__1.r = r__[i__4] * x[i__5].r, q__1.i = r__[i__4] * x[i__5] 00651 .i; 00652 x[i__2].r = q__1.r, x[i__2].i = q__1.i; 00653 /* L130: */ 00654 } 00655 /* L140: */ 00656 } 00657 i__1 = *nrhs; 00658 for (j = 1; j <= i__1; ++j) { 00659 ferr[j] /= rowcnd; 00660 /* L150: */ 00661 } 00662 } 00663 00664 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00665 00666 if (*rcond < slamch_("Epsilon")) { 00667 *info = *n + 1; 00668 } 00669 00670 rwork[1] = rpvgrw; 00671 return 0; 00672 00673 /* End of CGBSVX */ 00674 00675 } /* cgbsvx_ */