zporfs.c
Go to the documentation of this file.
00001 /* zporfs.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 doublecomplex c_b1 = {1.,0.};
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int zporfs_(char *uplo, integer *n, integer *nrhs, 
00022         doublecomplex *a, integer *lda, doublecomplex *af, integer *ldaf, 
00023         doublecomplex *b, integer *ldb, doublecomplex *x, integer *ldx, 
00024         doublereal *ferr, doublereal *berr, doublecomplex *work, doublereal *
00025         rwork, integer *info)
00026 {
00027     /* System generated locals */
00028     integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1, 
00029             x_offset, i__1, i__2, i__3, i__4, i__5;
00030     doublereal d__1, d__2, d__3, d__4;
00031     doublecomplex z__1;
00032 
00033     /* Builtin functions */
00034     double d_imag(doublecomplex *);
00035 
00036     /* Local variables */
00037     integer i__, j, k;
00038     doublereal s, xk;
00039     integer nz;
00040     doublereal eps;
00041     integer kase;
00042     doublereal safe1, safe2;
00043     extern logical lsame_(char *, char *);
00044     integer isave[3], count;
00045     extern /* Subroutine */ int zhemv_(char *, integer *, doublecomplex *, 
00046             doublecomplex *, integer *, doublecomplex *, integer *, 
00047             doublecomplex *, doublecomplex *, integer *);
00048     logical upper;
00049     extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
00050             doublecomplex *, integer *), zaxpy_(integer *, doublecomplex *, 
00051             doublecomplex *, integer *, doublecomplex *, integer *), zlacn2_(
00052             integer *, doublecomplex *, doublecomplex *, doublereal *, 
00053             integer *, integer *);
00054     extern doublereal dlamch_(char *);
00055     doublereal safmin;
00056     extern /* Subroutine */ int xerbla_(char *, integer *);
00057     doublereal lstres;
00058     extern /* Subroutine */ int zpotrs_(char *, integer *, integer *, 
00059             doublecomplex *, integer *, doublecomplex *, integer *, integer *);
00060 
00061 
00062 /*  -- LAPACK routine (version 3.2) -- */
00063 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00064 /*     November 2006 */
00065 
00066 /*     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH. */
00067 
00068 /*     .. Scalar Arguments .. */
00069 /*     .. */
00070 /*     .. Array Arguments .. */
00071 /*     .. */
00072 
00073 /*  Purpose */
00074 /*  ======= */
00075 
00076 /*  ZPORFS improves the computed solution to a system of linear */
00077 /*  equations when the coefficient matrix is Hermitian positive definite, */
00078 /*  and provides error bounds and backward error estimates for the */
00079 /*  solution. */
00080 
00081 /*  Arguments */
00082 /*  ========= */
00083 
00084 /*  UPLO    (input) CHARACTER*1 */
00085 /*          = 'U':  Upper triangle of A is stored; */
00086 /*          = 'L':  Lower triangle of A is stored. */
00087 
00088 /*  N       (input) INTEGER */
00089 /*          The order of the matrix A.  N >= 0. */
00090 
00091 /*  NRHS    (input) INTEGER */
00092 /*          The number of right hand sides, i.e., the number of columns */
00093 /*          of the matrices B and X.  NRHS >= 0. */
00094 
00095 /*  A       (input) COMPLEX*16 array, dimension (LDA,N) */
00096 /*          The Hermitian matrix A.  If UPLO = 'U', the leading N-by-N */
00097 /*          upper triangular part of A contains the upper triangular part */
00098 /*          of the matrix A, and the strictly lower triangular part of A */
00099 /*          is not referenced.  If UPLO = 'L', the leading N-by-N lower */
00100 /*          triangular part of A contains the lower triangular part of */
00101 /*          the matrix A, and the strictly upper triangular part of A is */
00102 /*          not referenced. */
00103 
00104 /*  LDA     (input) INTEGER */
00105 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00106 
00107 /*  AF      (input) COMPLEX*16 array, dimension (LDAF,N) */
00108 /*          The triangular factor U or L from the Cholesky factorization */
00109 /*          A = U**H*U or A = L*L**H, as computed by ZPOTRF. */
00110 
00111 /*  LDAF    (input) INTEGER */
00112 /*          The leading dimension of the array AF.  LDAF >= max(1,N). */
00113 
00114 /*  B       (input) COMPLEX*16 array, dimension (LDB,NRHS) */
00115 /*          The right hand side matrix B. */
00116 
00117 /*  LDB     (input) INTEGER */
00118 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00119 
00120 /*  X       (input/output) COMPLEX*16 array, dimension (LDX,NRHS) */
00121 /*          On entry, the solution matrix X, as computed by ZPOTRS. */
00122 /*          On exit, the improved solution matrix X. */
00123 
00124 /*  LDX     (input) INTEGER */
00125 /*          The leading dimension of the array X.  LDX >= max(1,N). */
00126 
00127 /*  FERR    (output) DOUBLE PRECISION array, dimension (NRHS) */
00128 /*          The estimated forward error bound for each solution vector */
00129 /*          X(j) (the j-th column of the solution matrix X). */
00130 /*          If XTRUE is the true solution corresponding to X(j), FERR(j) */
00131 /*          is an estimated upper bound for the magnitude of the largest */
00132 /*          element in (X(j) - XTRUE) divided by the magnitude of the */
00133 /*          largest element in X(j).  The estimate is as reliable as */
00134 /*          the estimate for RCOND, and is almost always a slight */
00135 /*          overestimate of the true error. */
00136 
00137 /*  BERR    (output) DOUBLE PRECISION array, dimension (NRHS) */
00138 /*          The componentwise relative backward error of each solution */
00139 /*          vector X(j) (i.e., the smallest relative change in */
00140 /*          any element of A or B that makes X(j) an exact solution). */
00141 
00142 /*  WORK    (workspace) COMPLEX*16 array, dimension (2*N) */
00143 
00144 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N) */
00145 
00146 /*  INFO    (output) INTEGER */
00147 /*          = 0:  successful exit */
00148 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00149 
00150 /*  Internal Parameters */
00151 /*  =================== */
00152 
00153 /*  ITMAX is the maximum number of steps of iterative refinement. */
00154 
00155 /*  ==================================================================== */
00156 
00157 /*     .. Parameters .. */
00158 /*     .. */
00159 /*     .. Local Scalars .. */
00160 /*     .. */
00161 /*     .. Local Arrays .. */
00162 /*     .. */
00163 /*     .. External Subroutines .. */
00164 /*     .. */
00165 /*     .. Intrinsic Functions .. */
00166 /*     .. */
00167 /*     .. External Functions .. */
00168 /*     .. */
00169 /*     .. Statement Functions .. */
00170 /*     .. */
00171 /*     .. Statement Function definitions .. */
00172 /*     .. */
00173 /*     .. Executable Statements .. */
00174 
00175 /*     Test the input parameters. */
00176 
00177     /* Parameter adjustments */
00178     a_dim1 = *lda;
00179     a_offset = 1 + a_dim1;
00180     a -= a_offset;
00181     af_dim1 = *ldaf;
00182     af_offset = 1 + af_dim1;
00183     af -= af_offset;
00184     b_dim1 = *ldb;
00185     b_offset = 1 + b_dim1;
00186     b -= b_offset;
00187     x_dim1 = *ldx;
00188     x_offset = 1 + x_dim1;
00189     x -= x_offset;
00190     --ferr;
00191     --berr;
00192     --work;
00193     --rwork;
00194 
00195     /* Function Body */
00196     *info = 0;
00197     upper = lsame_(uplo, "U");
00198     if (! upper && ! lsame_(uplo, "L")) {
00199         *info = -1;
00200     } else if (*n < 0) {
00201         *info = -2;
00202     } else if (*nrhs < 0) {
00203         *info = -3;
00204     } else if (*lda < max(1,*n)) {
00205         *info = -5;
00206     } else if (*ldaf < max(1,*n)) {
00207         *info = -7;
00208     } else if (*ldb < max(1,*n)) {
00209         *info = -9;
00210     } else if (*ldx < max(1,*n)) {
00211         *info = -11;
00212     }
00213     if (*info != 0) {
00214         i__1 = -(*info);
00215         xerbla_("ZPORFS", &i__1);
00216         return 0;
00217     }
00218 
00219 /*     Quick return if possible */
00220 
00221     if (*n == 0 || *nrhs == 0) {
00222         i__1 = *nrhs;
00223         for (j = 1; j <= i__1; ++j) {
00224             ferr[j] = 0.;
00225             berr[j] = 0.;
00226 /* L10: */
00227         }
00228         return 0;
00229     }
00230 
00231 /*     NZ = maximum number of nonzero elements in each row of A, plus 1 */
00232 
00233     nz = *n + 1;
00234     eps = dlamch_("Epsilon");
00235     safmin = dlamch_("Safe minimum");
00236     safe1 = nz * safmin;
00237     safe2 = safe1 / eps;
00238 
00239 /*     Do for each right hand side */
00240 
00241     i__1 = *nrhs;
00242     for (j = 1; j <= i__1; ++j) {
00243 
00244         count = 1;
00245         lstres = 3.;
00246 L20:
00247 
00248 /*        Loop until stopping criterion is satisfied. */
00249 
00250 /*        Compute residual R = B - A * X */
00251 
00252         zcopy_(n, &b[j * b_dim1 + 1], &c__1, &work[1], &c__1);
00253         z__1.r = -1., z__1.i = -0.;
00254         zhemv_(uplo, n, &z__1, &a[a_offset], lda, &x[j * x_dim1 + 1], &c__1, &
00255                 c_b1, &work[1], &c__1);
00256 
00257 /*        Compute componentwise relative backward error from formula */
00258 
00259 /*        max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) ) */
00260 
00261 /*        where abs(Z) is the componentwise absolute value of the matrix */
00262 /*        or vector Z.  If the i-th component of the denominator is less */
00263 /*        than SAFE2, then SAFE1 is added to the i-th components of the */
00264 /*        numerator and denominator before dividing. */
00265 
00266         i__2 = *n;
00267         for (i__ = 1; i__ <= i__2; ++i__) {
00268             i__3 = i__ + j * b_dim1;
00269             rwork[i__] = (d__1 = b[i__3].r, abs(d__1)) + (d__2 = d_imag(&b[
00270                     i__ + j * b_dim1]), abs(d__2));
00271 /* L30: */
00272         }
00273 
00274 /*        Compute abs(A)*abs(X) + abs(B). */
00275 
00276         if (upper) {
00277             i__2 = *n;
00278             for (k = 1; k <= i__2; ++k) {
00279                 s = 0.;
00280                 i__3 = k + j * x_dim1;
00281                 xk = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[k + j *
00282                          x_dim1]), abs(d__2));
00283                 i__3 = k - 1;
00284                 for (i__ = 1; i__ <= i__3; ++i__) {
00285                     i__4 = i__ + k * a_dim1;
00286                     rwork[i__] += ((d__1 = a[i__4].r, abs(d__1)) + (d__2 = 
00287                             d_imag(&a[i__ + k * a_dim1]), abs(d__2))) * xk;
00288                     i__4 = i__ + k * a_dim1;
00289                     i__5 = i__ + j * x_dim1;
00290                     s += ((d__1 = a[i__4].r, abs(d__1)) + (d__2 = d_imag(&a[
00291                             i__ + k * a_dim1]), abs(d__2))) * ((d__3 = x[i__5]
00292                             .r, abs(d__3)) + (d__4 = d_imag(&x[i__ + j * 
00293                             x_dim1]), abs(d__4)));
00294 /* L40: */
00295                 }
00296                 i__3 = k + k * a_dim1;
00297                 rwork[k] = rwork[k] + (d__1 = a[i__3].r, abs(d__1)) * xk + s;
00298 /* L50: */
00299             }
00300         } else {
00301             i__2 = *n;
00302             for (k = 1; k <= i__2; ++k) {
00303                 s = 0.;
00304                 i__3 = k + j * x_dim1;
00305                 xk = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[k + j *
00306                          x_dim1]), abs(d__2));
00307                 i__3 = k + k * a_dim1;
00308                 rwork[k] += (d__1 = a[i__3].r, abs(d__1)) * xk;
00309                 i__3 = *n;
00310                 for (i__ = k + 1; i__ <= i__3; ++i__) {
00311                     i__4 = i__ + k * a_dim1;
00312                     rwork[i__] += ((d__1 = a[i__4].r, abs(d__1)) + (d__2 = 
00313                             d_imag(&a[i__ + k * a_dim1]), abs(d__2))) * xk;
00314                     i__4 = i__ + k * a_dim1;
00315                     i__5 = i__ + j * x_dim1;
00316                     s += ((d__1 = a[i__4].r, abs(d__1)) + (d__2 = d_imag(&a[
00317                             i__ + k * a_dim1]), abs(d__2))) * ((d__3 = x[i__5]
00318                             .r, abs(d__3)) + (d__4 = d_imag(&x[i__ + j * 
00319                             x_dim1]), abs(d__4)));
00320 /* L60: */
00321                 }
00322                 rwork[k] += s;
00323 /* L70: */
00324             }
00325         }
00326         s = 0.;
00327         i__2 = *n;
00328         for (i__ = 1; i__ <= i__2; ++i__) {
00329             if (rwork[i__] > safe2) {
00330 /* Computing MAX */
00331                 i__3 = i__;
00332                 d__3 = s, d__4 = ((d__1 = work[i__3].r, abs(d__1)) + (d__2 = 
00333                         d_imag(&work[i__]), abs(d__2))) / rwork[i__];
00334                 s = max(d__3,d__4);
00335             } else {
00336 /* Computing MAX */
00337                 i__3 = i__;
00338                 d__3 = s, d__4 = ((d__1 = work[i__3].r, abs(d__1)) + (d__2 = 
00339                         d_imag(&work[i__]), abs(d__2)) + safe1) / (rwork[i__] 
00340                         + safe1);
00341                 s = max(d__3,d__4);
00342             }
00343 /* L80: */
00344         }
00345         berr[j] = s;
00346 
00347 /*        Test stopping criterion. Continue iterating if */
00348 /*           1) The residual BERR(J) is larger than machine epsilon, and */
00349 /*           2) BERR(J) decreased by at least a factor of 2 during the */
00350 /*              last iteration, and */
00351 /*           3) At most ITMAX iterations tried. */
00352 
00353         if (berr[j] > eps && berr[j] * 2. <= lstres && count <= 5) {
00354 
00355 /*           Update solution and try again. */
00356 
00357             zpotrs_(uplo, n, &c__1, &af[af_offset], ldaf, &work[1], n, info);
00358             zaxpy_(n, &c_b1, &work[1], &c__1, &x[j * x_dim1 + 1], &c__1);
00359             lstres = berr[j];
00360             ++count;
00361             goto L20;
00362         }
00363 
00364 /*        Bound error from formula */
00365 
00366 /*        norm(X - XTRUE) / norm(X) .le. FERR = */
00367 /*        norm( abs(inv(A))* */
00368 /*           ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X) */
00369 
00370 /*        where */
00371 /*          norm(Z) is the magnitude of the largest component of Z */
00372 /*          inv(A) is the inverse of A */
00373 /*          abs(Z) is the componentwise absolute value of the matrix or */
00374 /*             vector Z */
00375 /*          NZ is the maximum number of nonzeros in any row of A, plus 1 */
00376 /*          EPS is machine epsilon */
00377 
00378 /*        The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B)) */
00379 /*        is incremented by SAFE1 if the i-th component of */
00380 /*        abs(A)*abs(X) + abs(B) is less than SAFE2. */
00381 
00382 /*        Use ZLACN2 to estimate the infinity-norm of the matrix */
00383 /*           inv(A) * diag(W), */
00384 /*        where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) */
00385 
00386         i__2 = *n;
00387         for (i__ = 1; i__ <= i__2; ++i__) {
00388             if (rwork[i__] > safe2) {
00389                 i__3 = i__;
00390                 rwork[i__] = (d__1 = work[i__3].r, abs(d__1)) + (d__2 = 
00391                         d_imag(&work[i__]), abs(d__2)) + nz * eps * rwork[i__]
00392                         ;
00393             } else {
00394                 i__3 = i__;
00395                 rwork[i__] = (d__1 = work[i__3].r, abs(d__1)) + (d__2 = 
00396                         d_imag(&work[i__]), abs(d__2)) + nz * eps * rwork[i__]
00397                          + safe1;
00398             }
00399 /* L90: */
00400         }
00401 
00402         kase = 0;
00403 L100:
00404         zlacn2_(n, &work[*n + 1], &work[1], &ferr[j], &kase, isave);
00405         if (kase != 0) {
00406             if (kase == 1) {
00407 
00408 /*              Multiply by diag(W)*inv(A'). */
00409 
00410                 zpotrs_(uplo, n, &c__1, &af[af_offset], ldaf, &work[1], n, 
00411                         info);
00412                 i__2 = *n;
00413                 for (i__ = 1; i__ <= i__2; ++i__) {
00414                     i__3 = i__;
00415                     i__4 = i__;
00416                     i__5 = i__;
00417                     z__1.r = rwork[i__4] * work[i__5].r, z__1.i = rwork[i__4] 
00418                             * work[i__5].i;
00419                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00420 /* L110: */
00421                 }
00422             } else if (kase == 2) {
00423 
00424 /*              Multiply by inv(A)*diag(W). */
00425 
00426                 i__2 = *n;
00427                 for (i__ = 1; i__ <= i__2; ++i__) {
00428                     i__3 = i__;
00429                     i__4 = i__;
00430                     i__5 = i__;
00431                     z__1.r = rwork[i__4] * work[i__5].r, z__1.i = rwork[i__4] 
00432                             * work[i__5].i;
00433                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00434 /* L120: */
00435                 }
00436                 zpotrs_(uplo, n, &c__1, &af[af_offset], ldaf, &work[1], n, 
00437                         info);
00438             }
00439             goto L100;
00440         }
00441 
00442 /*        Normalize error. */
00443 
00444         lstres = 0.;
00445         i__2 = *n;
00446         for (i__ = 1; i__ <= i__2; ++i__) {
00447 /* Computing MAX */
00448             i__3 = i__ + j * x_dim1;
00449             d__3 = lstres, d__4 = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = 
00450                     d_imag(&x[i__ + j * x_dim1]), abs(d__2));
00451             lstres = max(d__3,d__4);
00452 /* L130: */
00453         }
00454         if (lstres != 0.) {
00455             ferr[j] /= lstres;
00456         }
00457 
00458 /* L140: */
00459     }
00460 
00461     return 0;
00462 
00463 /*     End of ZPORFS */
00464 
00465 } /* zporfs_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:43