cgtrfs.c
Go to the documentation of this file.
00001 /* cgtrfs.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 real c_b18 = -1.f;
00020 static real c_b19 = 1.f;
00021 static complex c_b26 = {1.f,0.f};
00022 
00023 /* Subroutine */ int cgtrfs_(char *trans, integer *n, integer *nrhs, complex *
00024         dl, complex *d__, complex *du, complex *dlf, complex *df, complex *
00025         duf, complex *du2, integer *ipiv, complex *b, integer *ldb, complex *
00026         x, integer *ldx, real *ferr, real *berr, complex *work, real *rwork, 
00027         integer *info)
00028 {
00029     /* System generated locals */
00030     integer b_dim1, b_offset, x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5, 
00031             i__6, i__7, i__8, i__9;
00032     real r__1, r__2, r__3, r__4, r__5, r__6, r__7, r__8, r__9, r__10, r__11, 
00033             r__12, r__13, r__14;
00034     complex q__1;
00035 
00036     /* Builtin functions */
00037     double r_imag(complex *);
00038 
00039     /* Local variables */
00040     integer i__, j;
00041     real s;
00042     integer nz;
00043     real eps;
00044     integer kase;
00045     real safe1, safe2;
00046     extern logical lsame_(char *, char *);
00047     integer isave[3];
00048     extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
00049             complex *, integer *), caxpy_(integer *, complex *, complex *, 
00050             integer *, complex *, integer *);
00051     integer count;
00052     extern /* Subroutine */ int clacn2_(integer *, complex *, complex *, real 
00053             *, integer *, integer *), clagtm_(char *, integer *, integer *, 
00054             real *, complex *, complex *, complex *, complex *, integer *, 
00055             real *, complex *, integer *);
00056     extern doublereal slamch_(char *);
00057     real safmin;
00058     extern /* Subroutine */ int xerbla_(char *, integer *);
00059     logical notran;
00060     char transn[1];
00061     extern /* Subroutine */ int cgttrs_(char *, integer *, integer *, complex 
00062             *, complex *, complex *, complex *, integer *, complex *, integer 
00063             *, integer *);
00064     char transt[1];
00065     real lstres;
00066 
00067 
00068 /*  -- LAPACK routine (version 3.2) -- */
00069 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00070 /*     November 2006 */
00071 
00072 /*     Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. */
00073 
00074 /*     .. Scalar Arguments .. */
00075 /*     .. */
00076 /*     .. Array Arguments .. */
00077 /*     .. */
00078 
00079 /*  Purpose */
00080 /*  ======= */
00081 
00082 /*  CGTRFS improves the computed solution to a system of linear */
00083 /*  equations when the coefficient matrix is tridiagonal, and provides */
00084 /*  error bounds and backward error estimates for the solution. */
00085 
00086 /*  Arguments */
00087 /*  ========= */
00088 
00089 /*  TRANS   (input) CHARACTER*1 */
00090 /*          Specifies the form of the system of equations: */
00091 /*          = 'N':  A * X = B     (No transpose) */
00092 /*          = 'T':  A**T * X = B  (Transpose) */
00093 /*          = 'C':  A**H * X = B  (Conjugate transpose) */
00094 
00095 /*  N       (input) INTEGER */
00096 /*          The order of the matrix A.  N >= 0. */
00097 
00098 /*  NRHS    (input) INTEGER */
00099 /*          The number of right hand sides, i.e., the number of columns */
00100 /*          of the matrix B.  NRHS >= 0. */
00101 
00102 /*  DL      (input) COMPLEX array, dimension (N-1) */
00103 /*          The (n-1) subdiagonal elements of A. */
00104 
00105 /*  D       (input) COMPLEX array, dimension (N) */
00106 /*          The diagonal elements of A. */
00107 
00108 /*  DU      (input) COMPLEX array, dimension (N-1) */
00109 /*          The (n-1) superdiagonal elements of A. */
00110 
00111 /*  DLF     (input) COMPLEX array, dimension (N-1) */
00112 /*          The (n-1) multipliers that define the matrix L from the */
00113 /*          LU factorization of A as computed by CGTTRF. */
00114 
00115 /*  DF      (input) COMPLEX array, dimension (N) */
00116 /*          The n diagonal elements of the upper triangular matrix U from */
00117 /*          the LU factorization of A. */
00118 
00119 /*  DUF     (input) COMPLEX array, dimension (N-1) */
00120 /*          The (n-1) elements of the first superdiagonal of U. */
00121 
00122 /*  DU2     (input) COMPLEX array, dimension (N-2) */
00123 /*          The (n-2) elements of the second superdiagonal of U. */
00124 
00125 /*  IPIV    (input) INTEGER array, dimension (N) */
00126 /*          The pivot indices; for 1 <= i <= n, row i of the matrix was */
00127 /*          interchanged with row IPIV(i).  IPIV(i) will always be either */
00128 /*          i or i+1; IPIV(i) = i indicates a row interchange was not */
00129 /*          required. */
00130 
00131 /*  B       (input) COMPLEX array, dimension (LDB,NRHS) */
00132 /*          The right hand side matrix B. */
00133 
00134 /*  LDB     (input) INTEGER */
00135 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00136 
00137 /*  X       (input/output) COMPLEX array, dimension (LDX,NRHS) */
00138 /*          On entry, the solution matrix X, as computed by CGTTRS. */
00139 /*          On exit, the improved solution matrix X. */
00140 
00141 /*  LDX     (input) INTEGER */
00142 /*          The leading dimension of the array X.  LDX >= max(1,N). */
00143 
00144 /*  FERR    (output) REAL array, dimension (NRHS) */
00145 /*          The estimated forward error bound for each solution vector */
00146 /*          X(j) (the j-th column of the solution matrix X). */
00147 /*          If XTRUE is the true solution corresponding to X(j), FERR(j) */
00148 /*          is an estimated upper bound for the magnitude of the largest */
00149 /*          element in (X(j) - XTRUE) divided by the magnitude of the */
00150 /*          largest element in X(j).  The estimate is as reliable as */
00151 /*          the estimate for RCOND, and is almost always a slight */
00152 /*          overestimate of the true error. */
00153 
00154 /*  BERR    (output) REAL array, dimension (NRHS) */
00155 /*          The componentwise relative backward error of each solution */
00156 /*          vector X(j) (i.e., the smallest relative change in */
00157 /*          any element of A or B that makes X(j) an exact solution). */
00158 
00159 /*  WORK    (workspace) COMPLEX array, dimension (2*N) */
00160 
00161 /*  RWORK   (workspace) REAL array, dimension (N) */
00162 
00163 /*  INFO    (output) INTEGER */
00164 /*          = 0:  successful exit */
00165 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00166 
00167 /*  Internal Parameters */
00168 /*  =================== */
00169 
00170 /*  ITMAX is the maximum number of steps of iterative refinement. */
00171 
00172 /*  ===================================================================== */
00173 
00174 /*     .. Parameters .. */
00175 /*     .. */
00176 /*     .. Local Scalars .. */
00177 /*     .. */
00178 /*     .. Local Arrays .. */
00179 /*     .. */
00180 /*     .. External Subroutines .. */
00181 /*     .. */
00182 /*     .. Intrinsic Functions .. */
00183 /*     .. */
00184 /*     .. External Functions .. */
00185 /*     .. */
00186 /*     .. Statement Functions .. */
00187 /*     .. */
00188 /*     .. Statement Function definitions .. */
00189 /*     .. */
00190 /*     .. Executable Statements .. */
00191 
00192 /*     Test the input parameters. */
00193 
00194     /* Parameter adjustments */
00195     --dl;
00196     --d__;
00197     --du;
00198     --dlf;
00199     --df;
00200     --duf;
00201     --du2;
00202     --ipiv;
00203     b_dim1 = *ldb;
00204     b_offset = 1 + b_dim1;
00205     b -= b_offset;
00206     x_dim1 = *ldx;
00207     x_offset = 1 + x_dim1;
00208     x -= x_offset;
00209     --ferr;
00210     --berr;
00211     --work;
00212     --rwork;
00213 
00214     /* Function Body */
00215     *info = 0;
00216     notran = lsame_(trans, "N");
00217     if (! notran && ! lsame_(trans, "T") && ! lsame_(
00218             trans, "C")) {
00219         *info = -1;
00220     } else if (*n < 0) {
00221         *info = -2;
00222     } else if (*nrhs < 0) {
00223         *info = -3;
00224     } else if (*ldb < max(1,*n)) {
00225         *info = -13;
00226     } else if (*ldx < max(1,*n)) {
00227         *info = -15;
00228     }
00229     if (*info != 0) {
00230         i__1 = -(*info);
00231         xerbla_("CGTRFS", &i__1);
00232         return 0;
00233     }
00234 
00235 /*     Quick return if possible */
00236 
00237     if (*n == 0 || *nrhs == 0) {
00238         i__1 = *nrhs;
00239         for (j = 1; j <= i__1; ++j) {
00240             ferr[j] = 0.f;
00241             berr[j] = 0.f;
00242 /* L10: */
00243         }
00244         return 0;
00245     }
00246 
00247     if (notran) {
00248         *(unsigned char *)transn = 'N';
00249         *(unsigned char *)transt = 'C';
00250     } else {
00251         *(unsigned char *)transn = 'C';
00252         *(unsigned char *)transt = 'N';
00253     }
00254 
00255 /*     NZ = maximum number of nonzero elements in each row of A, plus 1 */
00256 
00257     nz = 4;
00258     eps = slamch_("Epsilon");
00259     safmin = slamch_("Safe minimum");
00260     safe1 = nz * safmin;
00261     safe2 = safe1 / eps;
00262 
00263 /*     Do for each right hand side */
00264 
00265     i__1 = *nrhs;
00266     for (j = 1; j <= i__1; ++j) {
00267 
00268         count = 1;
00269         lstres = 3.f;
00270 L20:
00271 
00272 /*        Loop until stopping criterion is satisfied. */
00273 
00274 /*        Compute residual R = B - op(A) * X, */
00275 /*        where op(A) = A, A**T, or A**H, depending on TRANS. */
00276 
00277         ccopy_(n, &b[j * b_dim1 + 1], &c__1, &work[1], &c__1);
00278         clagtm_(trans, n, &c__1, &c_b18, &dl[1], &d__[1], &du[1], &x[j * 
00279                 x_dim1 + 1], ldx, &c_b19, &work[1], n);
00280 
00281 /*        Compute abs(op(A))*abs(x) + abs(b) for use in the backward */
00282 /*        error bound. */
00283 
00284         if (notran) {
00285             if (*n == 1) {
00286                 i__2 = j * b_dim1 + 1;
00287                 i__3 = j * x_dim1 + 1;
00288                 rwork[1] = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[
00289                         j * b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, 
00290                         dabs(r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * 
00291                         ((r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[j 
00292                         * x_dim1 + 1]), dabs(r__6)));
00293             } else {
00294                 i__2 = j * b_dim1 + 1;
00295                 i__3 = j * x_dim1 + 1;
00296                 i__4 = j * x_dim1 + 2;
00297                 rwork[1] = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[
00298                         j * b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, 
00299                         dabs(r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * 
00300                         ((r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[j 
00301                         * x_dim1 + 1]), dabs(r__6))) + ((r__7 = du[1].r, dabs(
00302                         r__7)) + (r__8 = r_imag(&du[1]), dabs(r__8))) * ((
00303                         r__9 = x[i__4].r, dabs(r__9)) + (r__10 = r_imag(&x[j *
00304                          x_dim1 + 2]), dabs(r__10)));
00305                 i__2 = *n - 1;
00306                 for (i__ = 2; i__ <= i__2; ++i__) {
00307                     i__3 = i__ + j * b_dim1;
00308                     i__4 = i__ - 1;
00309                     i__5 = i__ - 1 + j * x_dim1;
00310                     i__6 = i__;
00311                     i__7 = i__ + j * x_dim1;
00312                     i__8 = i__;
00313                     i__9 = i__ + 1 + j * x_dim1;
00314                     rwork[i__] = (r__1 = b[i__3].r, dabs(r__1)) + (r__2 = 
00315                             r_imag(&b[i__ + j * b_dim1]), dabs(r__2)) + ((
00316                             r__3 = dl[i__4].r, dabs(r__3)) + (r__4 = r_imag(&
00317                             dl[i__ - 1]), dabs(r__4))) * ((r__5 = x[i__5].r, 
00318                             dabs(r__5)) + (r__6 = r_imag(&x[i__ - 1 + j * 
00319                             x_dim1]), dabs(r__6))) + ((r__7 = d__[i__6].r, 
00320                             dabs(r__7)) + (r__8 = r_imag(&d__[i__]), dabs(
00321                             r__8))) * ((r__9 = x[i__7].r, dabs(r__9)) + (
00322                             r__10 = r_imag(&x[i__ + j * x_dim1]), dabs(r__10))
00323                             ) + ((r__11 = du[i__8].r, dabs(r__11)) + (r__12 = 
00324                             r_imag(&du[i__]), dabs(r__12))) * ((r__13 = x[
00325                             i__9].r, dabs(r__13)) + (r__14 = r_imag(&x[i__ + 
00326                             1 + j * x_dim1]), dabs(r__14)));
00327 /* L30: */
00328                 }
00329                 i__2 = *n + j * b_dim1;
00330                 i__3 = *n - 1;
00331                 i__4 = *n - 1 + j * x_dim1;
00332                 i__5 = *n;
00333                 i__6 = *n + j * x_dim1;
00334                 rwork[*n] = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&
00335                         b[*n + j * b_dim1]), dabs(r__2)) + ((r__3 = dl[i__3]
00336                         .r, dabs(r__3)) + (r__4 = r_imag(&dl[*n - 1]), dabs(
00337                         r__4))) * ((r__5 = x[i__4].r, dabs(r__5)) + (r__6 = 
00338                         r_imag(&x[*n - 1 + j * x_dim1]), dabs(r__6))) + ((
00339                         r__7 = d__[i__5].r, dabs(r__7)) + (r__8 = r_imag(&d__[
00340                         *n]), dabs(r__8))) * ((r__9 = x[i__6].r, dabs(r__9)) 
00341                         + (r__10 = r_imag(&x[*n + j * x_dim1]), dabs(r__10)));
00342             }
00343         } else {
00344             if (*n == 1) {
00345                 i__2 = j * b_dim1 + 1;
00346                 i__3 = j * x_dim1 + 1;
00347                 rwork[1] = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[
00348                         j * b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, 
00349                         dabs(r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * 
00350                         ((r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[j 
00351                         * x_dim1 + 1]), dabs(r__6)));
00352             } else {
00353                 i__2 = j * b_dim1 + 1;
00354                 i__3 = j * x_dim1 + 1;
00355                 i__4 = j * x_dim1 + 2;
00356                 rwork[1] = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[
00357                         j * b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, 
00358                         dabs(r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * 
00359                         ((r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[j 
00360                         * x_dim1 + 1]), dabs(r__6))) + ((r__7 = dl[1].r, dabs(
00361                         r__7)) + (r__8 = r_imag(&dl[1]), dabs(r__8))) * ((
00362                         r__9 = x[i__4].r, dabs(r__9)) + (r__10 = r_imag(&x[j *
00363                          x_dim1 + 2]), dabs(r__10)));
00364                 i__2 = *n - 1;
00365                 for (i__ = 2; i__ <= i__2; ++i__) {
00366                     i__3 = i__ + j * b_dim1;
00367                     i__4 = i__ - 1;
00368                     i__5 = i__ - 1 + j * x_dim1;
00369                     i__6 = i__;
00370                     i__7 = i__ + j * x_dim1;
00371                     i__8 = i__;
00372                     i__9 = i__ + 1 + j * x_dim1;
00373                     rwork[i__] = (r__1 = b[i__3].r, dabs(r__1)) + (r__2 = 
00374                             r_imag(&b[i__ + j * b_dim1]), dabs(r__2)) + ((
00375                             r__3 = du[i__4].r, dabs(r__3)) + (r__4 = r_imag(&
00376                             du[i__ - 1]), dabs(r__4))) * ((r__5 = x[i__5].r, 
00377                             dabs(r__5)) + (r__6 = r_imag(&x[i__ - 1 + j * 
00378                             x_dim1]), dabs(r__6))) + ((r__7 = d__[i__6].r, 
00379                             dabs(r__7)) + (r__8 = r_imag(&d__[i__]), dabs(
00380                             r__8))) * ((r__9 = x[i__7].r, dabs(r__9)) + (
00381                             r__10 = r_imag(&x[i__ + j * x_dim1]), dabs(r__10))
00382                             ) + ((r__11 = dl[i__8].r, dabs(r__11)) + (r__12 = 
00383                             r_imag(&dl[i__]), dabs(r__12))) * ((r__13 = x[
00384                             i__9].r, dabs(r__13)) + (r__14 = r_imag(&x[i__ + 
00385                             1 + j * x_dim1]), dabs(r__14)));
00386 /* L40: */
00387                 }
00388                 i__2 = *n + j * b_dim1;
00389                 i__3 = *n - 1;
00390                 i__4 = *n - 1 + j * x_dim1;
00391                 i__5 = *n;
00392                 i__6 = *n + j * x_dim1;
00393                 rwork[*n] = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&
00394                         b[*n + j * b_dim1]), dabs(r__2)) + ((r__3 = du[i__3]
00395                         .r, dabs(r__3)) + (r__4 = r_imag(&du[*n - 1]), dabs(
00396                         r__4))) * ((r__5 = x[i__4].r, dabs(r__5)) + (r__6 = 
00397                         r_imag(&x[*n - 1 + j * x_dim1]), dabs(r__6))) + ((
00398                         r__7 = d__[i__5].r, dabs(r__7)) + (r__8 = r_imag(&d__[
00399                         *n]), dabs(r__8))) * ((r__9 = x[i__6].r, dabs(r__9)) 
00400                         + (r__10 = r_imag(&x[*n + j * x_dim1]), dabs(r__10)));
00401             }
00402         }
00403 
00404 /*        Compute componentwise relative backward error from formula */
00405 
00406 /*        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) */
00407 
00408 /*        where abs(Z) is the componentwise absolute value of the matrix */
00409 /*        or vector Z.  If the i-th component of the denominator is less */
00410 /*        than SAFE2, then SAFE1 is added to the i-th components of the */
00411 /*        numerator and denominator before dividing. */
00412 
00413         s = 0.f;
00414         i__2 = *n;
00415         for (i__ = 1; i__ <= i__2; ++i__) {
00416             if (rwork[i__] > safe2) {
00417 /* Computing MAX */
00418                 i__3 = i__;
00419                 r__3 = s, r__4 = ((r__1 = work[i__3].r, dabs(r__1)) + (r__2 = 
00420                         r_imag(&work[i__]), dabs(r__2))) / rwork[i__];
00421                 s = dmax(r__3,r__4);
00422             } else {
00423 /* Computing MAX */
00424                 i__3 = i__;
00425                 r__3 = s, r__4 = ((r__1 = work[i__3].r, dabs(r__1)) + (r__2 = 
00426                         r_imag(&work[i__]), dabs(r__2)) + safe1) / (rwork[i__]
00427                          + safe1);
00428                 s = dmax(r__3,r__4);
00429             }
00430 /* L50: */
00431         }
00432         berr[j] = s;
00433 
00434 /*        Test stopping criterion. Continue iterating if */
00435 /*           1) The residual BERR(J) is larger than machine epsilon, and */
00436 /*           2) BERR(J) decreased by at least a factor of 2 during the */
00437 /*              last iteration, and */
00438 /*           3) At most ITMAX iterations tried. */
00439 
00440         if (berr[j] > eps && berr[j] * 2.f <= lstres && count <= 5) {
00441 
00442 /*           Update solution and try again. */
00443 
00444             cgttrs_(trans, n, &c__1, &dlf[1], &df[1], &duf[1], &du2[1], &ipiv[
00445                     1], &work[1], n, info);
00446             caxpy_(n, &c_b26, &work[1], &c__1, &x[j * x_dim1 + 1], &c__1);
00447             lstres = berr[j];
00448             ++count;
00449             goto L20;
00450         }
00451 
00452 /*        Bound error from formula */
00453 
00454 /*        norm(X - XTRUE) / norm(X) .le. FERR = */
00455 /*        norm( abs(inv(op(A)))* */
00456 /*           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) */
00457 
00458 /*        where */
00459 /*          norm(Z) is the magnitude of the largest component of Z */
00460 /*          inv(op(A)) is the inverse of op(A) */
00461 /*          abs(Z) is the componentwise absolute value of the matrix or */
00462 /*             vector Z */
00463 /*          NZ is the maximum number of nonzeros in any row of A, plus 1 */
00464 /*          EPS is machine epsilon */
00465 
00466 /*        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) */
00467 /*        is incremented by SAFE1 if the i-th component of */
00468 /*        abs(op(A))*abs(X) + abs(B) is less than SAFE2. */
00469 
00470 /*        Use CLACN2 to estimate the infinity-norm of the matrix */
00471 /*           inv(op(A)) * diag(W), */
00472 /*        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) */
00473 
00474         i__2 = *n;
00475         for (i__ = 1; i__ <= i__2; ++i__) {
00476             if (rwork[i__] > safe2) {
00477                 i__3 = i__;
00478                 rwork[i__] = (r__1 = work[i__3].r, dabs(r__1)) + (r__2 = 
00479                         r_imag(&work[i__]), dabs(r__2)) + nz * eps * rwork[
00480                         i__];
00481             } else {
00482                 i__3 = i__;
00483                 rwork[i__] = (r__1 = work[i__3].r, dabs(r__1)) + (r__2 = 
00484                         r_imag(&work[i__]), dabs(r__2)) + nz * eps * rwork[
00485                         i__] + safe1;
00486             }
00487 /* L60: */
00488         }
00489 
00490         kase = 0;
00491 L70:
00492         clacn2_(n, &work[*n + 1], &work[1], &ferr[j], &kase, isave);
00493         if (kase != 0) {
00494             if (kase == 1) {
00495 
00496 /*              Multiply by diag(W)*inv(op(A)**H). */
00497 
00498                 cgttrs_(transt, n, &c__1, &dlf[1], &df[1], &duf[1], &du2[1], &
00499                         ipiv[1], &work[1], n, info);
00500                 i__2 = *n;
00501                 for (i__ = 1; i__ <= i__2; ++i__) {
00502                     i__3 = i__;
00503                     i__4 = i__;
00504                     i__5 = i__;
00505                     q__1.r = rwork[i__4] * work[i__5].r, q__1.i = rwork[i__4] 
00506                             * work[i__5].i;
00507                     work[i__3].r = q__1.r, work[i__3].i = q__1.i;
00508 /* L80: */
00509                 }
00510             } else {
00511 
00512 /*              Multiply by inv(op(A))*diag(W). */
00513 
00514                 i__2 = *n;
00515                 for (i__ = 1; i__ <= i__2; ++i__) {
00516                     i__3 = i__;
00517                     i__4 = i__;
00518                     i__5 = i__;
00519                     q__1.r = rwork[i__4] * work[i__5].r, q__1.i = rwork[i__4] 
00520                             * work[i__5].i;
00521                     work[i__3].r = q__1.r, work[i__3].i = q__1.i;
00522 /* L90: */
00523                 }
00524                 cgttrs_(transn, n, &c__1, &dlf[1], &df[1], &duf[1], &du2[1], &
00525                         ipiv[1], &work[1], n, info);
00526             }
00527             goto L70;
00528         }
00529 
00530 /*        Normalize error. */
00531 
00532         lstres = 0.f;
00533         i__2 = *n;
00534         for (i__ = 1; i__ <= i__2; ++i__) {
00535 /* Computing MAX */
00536             i__3 = i__ + j * x_dim1;
00537             r__3 = lstres, r__4 = (r__1 = x[i__3].r, dabs(r__1)) + (r__2 = 
00538                     r_imag(&x[i__ + j * x_dim1]), dabs(r__2));
00539             lstres = dmax(r__3,r__4);
00540 /* L100: */
00541         }
00542         if (lstres != 0.f) {
00543             ferr[j] /= lstres;
00544         }
00545 
00546 /* L110: */
00547     }
00548 
00549     return 0;
00550 
00551 /*     End of CGTRFS */
00552 
00553 } /* cgtrfs_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:27