dget03.c
Go to the documentation of this file.
00001 /* dget03.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 doublereal c_b7 = -1.;
00019 static doublereal c_b8 = 0.;
00020 
00021 /* Subroutine */ int dget03_(integer *n, doublereal *a, integer *lda, 
00022         doublereal *ainv, integer *ldainv, doublereal *work, integer *ldwork, 
00023         doublereal *rwork, doublereal *rcond, doublereal *resid)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, ainv_dim1, ainv_offset, work_dim1, work_offset, 
00027             i__1;
00028 
00029     /* Local variables */
00030     integer i__;
00031     doublereal eps;
00032     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
00033             integer *, doublereal *, doublereal *, integer *, doublereal *, 
00034             integer *, doublereal *, doublereal *, integer *);
00035     doublereal anorm;
00036     extern doublereal dlamch_(char *), dlange_(char *, integer *, 
00037             integer *, doublereal *, integer *, doublereal *);
00038     doublereal ainvnm;
00039 
00040 
00041 /*  -- LAPACK test routine (version 3.1) -- */
00042 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00043 /*     November 2006 */
00044 
00045 /*     .. Scalar Arguments .. */
00046 /*     .. */
00047 /*     .. Array Arguments .. */
00048 /*     .. */
00049 
00050 /*  Purpose */
00051 /*  ======= */
00052 
00053 /*  DGET03 computes the residual for a general matrix times its inverse: */
00054 /*     norm( I - AINV*A ) / ( N * norm(A) * norm(AINV) * EPS ), */
00055 /*  where EPS is the machine epsilon. */
00056 
00057 /*  Arguments */
00058 /*  ========== */
00059 
00060 /*  N       (input) INTEGER */
00061 /*          The number of rows and columns of the matrix A.  N >= 0. */
00062 
00063 /*  A       (input) DOUBLE PRECISION array, dimension (LDA,N) */
00064 /*          The original N x N matrix A. */
00065 
00066 /*  LDA     (input) INTEGER */
00067 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00068 
00069 /*  AINV    (input) DOUBLE PRECISION array, dimension (LDAINV,N) */
00070 /*          The inverse of the matrix A. */
00071 
00072 /*  LDAINV  (input) INTEGER */
00073 /*          The leading dimension of the array AINV.  LDAINV >= max(1,N). */
00074 
00075 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,N) */
00076 
00077 /*  LDWORK  (input) INTEGER */
00078 /*          The leading dimension of the array WORK.  LDWORK >= max(1,N). */
00079 
00080 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N) */
00081 
00082 /*  RCOND   (output) DOUBLE PRECISION */
00083 /*          The reciprocal of the condition number of A, computed as */
00084 /*          ( 1/norm(A) ) / norm(AINV). */
00085 
00086 /*  RESID   (output) DOUBLE PRECISION */
00087 /*          norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS ) */
00088 
00089 /*  ===================================================================== */
00090 
00091 /*     .. Parameters .. */
00092 /*     .. */
00093 /*     .. Local Scalars .. */
00094 /*     .. */
00095 /*     .. External Functions .. */
00096 /*     .. */
00097 /*     .. External Subroutines .. */
00098 /*     .. */
00099 /*     .. Intrinsic Functions .. */
00100 /*     .. */
00101 /*     .. Executable Statements .. */
00102 
00103 /*     Quick exit if N = 0. */
00104 
00105     /* Parameter adjustments */
00106     a_dim1 = *lda;
00107     a_offset = 1 + a_dim1;
00108     a -= a_offset;
00109     ainv_dim1 = *ldainv;
00110     ainv_offset = 1 + ainv_dim1;
00111     ainv -= ainv_offset;
00112     work_dim1 = *ldwork;
00113     work_offset = 1 + work_dim1;
00114     work -= work_offset;
00115     --rwork;
00116 
00117     /* Function Body */
00118     if (*n <= 0) {
00119         *rcond = 1.;
00120         *resid = 0.;
00121         return 0;
00122     }
00123 
00124 /*     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. */
00125 
00126     eps = dlamch_("Epsilon");
00127     anorm = dlange_("1", n, n, &a[a_offset], lda, &rwork[1]);
00128     ainvnm = dlange_("1", n, n, &ainv[ainv_offset], ldainv, &rwork[1]);
00129     if (anorm <= 0. || ainvnm <= 0.) {
00130         *rcond = 0.;
00131         *resid = 1. / eps;
00132         return 0;
00133     }
00134     *rcond = 1. / anorm / ainvnm;
00135 
00136 /*     Compute I - A * AINV */
00137 
00138     dgemm_("No transpose", "No transpose", n, n, n, &c_b7, &ainv[ainv_offset], 
00139              ldainv, &a[a_offset], lda, &c_b8, &work[work_offset], ldwork);
00140     i__1 = *n;
00141     for (i__ = 1; i__ <= i__1; ++i__) {
00142         work[i__ + i__ * work_dim1] += 1.;
00143 /* L10: */
00144     }
00145 
00146 /*     Compute norm(I - AINV*A) / (N * norm(A) * norm(AINV) * EPS) */
00147 
00148     *resid = dlange_("1", n, n, &work[work_offset], ldwork, &rwork[1]);
00149 
00150     *resid = *resid * *rcond / eps / (doublereal) (*n);
00151 
00152     return 0;
00153 
00154 /*     End of DGET03 */
00155 
00156 } /* dget03_ */


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