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


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