zqrt16.c
Go to the documentation of this file.
00001 /* zqrt16.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 zqrt16_(char *trans, integer *m, integer *n, integer *
00022         nrhs, doublecomplex *a, integer *lda, doublecomplex *x, integer *ldx, 
00023         doublecomplex *b, integer *ldb, doublereal *rwork, doublereal *resid)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, i__1;
00027     doublereal d__1, d__2;
00028     doublecomplex z__1;
00029 
00030     /* Local variables */
00031     integer j, n1, n2;
00032     doublereal eps;
00033     extern logical lsame_(char *, char *);
00034     doublereal anorm, bnorm;
00035     extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 
00036             integer *, doublecomplex *, doublecomplex *, integer *, 
00037             doublecomplex *, integer *, doublecomplex *, doublecomplex *, 
00038             integer *);
00039     doublereal xnorm;
00040     extern doublereal dlamch_(char *), zlange_(char *, integer *, 
00041             integer *, doublecomplex *, integer *, doublereal *), 
00042             dzasum_(integer *, doublecomplex *, integer *);
00043 
00044 
00045 /*  -- LAPACK test routine (version 3.1) -- */
00046 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00047 /*     November 2006 */
00048 
00049 /*     .. Scalar Arguments .. */
00050 /*     .. */
00051 /*     .. Array Arguments .. */
00052 /*     .. */
00053 
00054 /*  Purpose */
00055 /*  ======= */
00056 
00057 /*  ZQRT16 computes the residual for a solution of a system of linear */
00058 /*  equations  A*x = b  or  A'*x = b: */
00059 /*     RESID = norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ), */
00060 /*  where EPS is the machine epsilon. */
00061 
00062 /*  Arguments */
00063 /*  ========= */
00064 
00065 /*  TRANS   (input) CHARACTER*1 */
00066 /*          Specifies the form of the system of equations: */
00067 /*          = 'N':  A *x = b */
00068 /*          = 'T':  A^T*x = b, where A^T is the transpose of A */
00069 /*          = 'C':  A^H*x = b, where A^H is the conjugate transpose of A */
00070 
00071 /*  M       (input) INTEGER */
00072 /*          The number of rows of the matrix A.  M >= 0. */
00073 
00074 /*  N       (input) INTEGER */
00075 /*          The number of columns of the matrix A.  N >= 0. */
00076 
00077 /*  NRHS    (input) INTEGER */
00078 /*          The number of columns of B, the matrix of right hand sides. */
00079 /*          NRHS >= 0. */
00080 
00081 /*  A       (input) COMPLEX*16 array, dimension (LDA,N) */
00082 /*          The original M x N matrix A. */
00083 
00084 /*  LDA     (input) INTEGER */
00085 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00086 
00087 /*  X       (input) COMPLEX*16 array, dimension (LDX,NRHS) */
00088 /*          The computed solution vectors for the system of linear */
00089 /*          equations. */
00090 
00091 /*  LDX     (input) INTEGER */
00092 /*          The leading dimension of the array X.  If TRANS = 'N', */
00093 /*          LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M). */
00094 
00095 /*  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS) */
00096 /*          On entry, the right hand side vectors for the system of */
00097 /*          linear equations. */
00098 /*          On exit, B is overwritten with the difference B - A*X. */
00099 
00100 /*  LDB     (input) INTEGER */
00101 /*          The leading dimension of the array B.  IF TRANS = 'N', */
00102 /*          LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). */
00103 
00104 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (M) */
00105 
00106 /*  RESID   (output) DOUBLE PRECISION */
00107 /*          The maximum over the number of right hand sides of */
00108 /*          norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ). */
00109 
00110 /*  ===================================================================== */
00111 
00112 /*     .. Parameters .. */
00113 /*     .. */
00114 /*     .. Local Scalars .. */
00115 /*     .. */
00116 /*     .. External Functions .. */
00117 /*     .. */
00118 /*     .. External Subroutines .. */
00119 /*     .. */
00120 /*     .. Intrinsic Functions .. */
00121 /*     .. */
00122 /*     .. Executable Statements .. */
00123 
00124 /*     Quick exit if M = 0 or N = 0 or NRHS = 0 */
00125 
00126     /* Parameter adjustments */
00127     a_dim1 = *lda;
00128     a_offset = 1 + a_dim1;
00129     a -= a_offset;
00130     x_dim1 = *ldx;
00131     x_offset = 1 + x_dim1;
00132     x -= x_offset;
00133     b_dim1 = *ldb;
00134     b_offset = 1 + b_dim1;
00135     b -= b_offset;
00136     --rwork;
00137 
00138     /* Function Body */
00139     if (*m <= 0 || *n <= 0 || *nrhs == 0) {
00140         *resid = 0.;
00141         return 0;
00142     }
00143 
00144     if (lsame_(trans, "T") || lsame_(trans, "C")) {
00145         anorm = zlange_("I", m, n, &a[a_offset], lda, &rwork[1]);
00146         n1 = *n;
00147         n2 = *m;
00148     } else {
00149         anorm = zlange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00150         n1 = *m;
00151         n2 = *n;
00152     }
00153 
00154     eps = dlamch_("Epsilon");
00155 
00156 /*     Compute  B - A*X  (or  B - A'*X ) and store in B. */
00157 
00158     z__1.r = -1., z__1.i = -0.;
00159     zgemm_(trans, "No transpose", &n1, nrhs, &n2, &z__1, &a[a_offset], lda, &
00160             x[x_offset], ldx, &c_b1, &b[b_offset], ldb)
00161             ;
00162 
00163 /*     Compute the maximum over the number of right hand sides of */
00164 /*        norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ) . */
00165 
00166     *resid = 0.;
00167     i__1 = *nrhs;
00168     for (j = 1; j <= i__1; ++j) {
00169         bnorm = dzasum_(&n1, &b[j * b_dim1 + 1], &c__1);
00170         xnorm = dzasum_(&n2, &x[j * x_dim1 + 1], &c__1);
00171         if (anorm == 0. && bnorm == 0.) {
00172             *resid = 0.;
00173         } else if (anorm <= 0. || xnorm <= 0.) {
00174             *resid = 1. / eps;
00175         } else {
00176 /* Computing MAX */
00177             d__1 = *resid, d__2 = bnorm / anorm / xnorm / (max(*m,*n) * eps);
00178             *resid = max(d__1,d__2);
00179         }
00180 /* L10: */
00181     }
00182 
00183     return 0;
00184 
00185 /*     End of ZQRT16 */
00186 
00187 } /* zqrt16_ */


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