zget52.c
Go to the documentation of this file.
00001 /* zget52.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 = {0.,0.};
00019 static doublecomplex c_b2 = {1.,0.};
00020 static integer c__1 = 1;
00021 
00022 /* Subroutine */ int zget52_(logical *left, integer *n, doublecomplex *a, 
00023         integer *lda, doublecomplex *b, integer *ldb, doublecomplex *e, 
00024         integer *lde, doublecomplex *alpha, doublecomplex *beta, 
00025         doublecomplex *work, doublereal *rwork, doublereal *result)
00026 {
00027     /* System generated locals */
00028     integer a_dim1, a_offset, b_dim1, b_offset, e_dim1, e_offset, i__1, i__2, 
00029             i__3;
00030     doublereal d__1, d__2, d__3, d__4, d__5, d__6;
00031     doublecomplex z__1;
00032 
00033     /* Builtin functions */
00034     double d_imag(doublecomplex *);
00035     void d_cnjg(doublecomplex *, doublecomplex *);
00036 
00037     /* Local variables */
00038     integer j;
00039     doublereal ulp;
00040     integer jvec;
00041     doublereal temp1;
00042     doublecomplex betai;
00043     doublereal scale, abmax, anorm, bnorm, enorm;
00044     char trans[1];
00045     extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 
00046             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00047             integer *, doublecomplex *, doublecomplex *, integer *);
00048     doublecomplex acoeff, bcoeff;
00049     extern doublereal dlamch_(char *);
00050     doublecomplex alphai;
00051     doublereal alfmax, safmin;
00052     char normab[1];
00053     doublereal safmax, betmax;
00054     extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 
00055             integer *, doublereal *);
00056     doublereal enrmer, errnrm;
00057 
00058 
00059 /*  -- LAPACK test routine (version 3.1) -- */
00060 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00061 /*     November 2006 */
00062 
00063 /*     .. Scalar Arguments .. */
00064 /*     .. */
00065 /*     .. Array Arguments .. */
00066 /*     .. */
00067 
00068 /*  Purpose */
00069 /*  ======= */
00070 
00071 /*  ZGET52  does an eigenvector check for the generalized eigenvalue */
00072 /*  problem. */
00073 
00074 /*  The basic test for right eigenvectors is: */
00075 
00076 /*                            | b(i) A E(i) -  a(i) B E(i) | */
00077 /*          RESULT(1) = max   ------------------------------- */
00078 /*                       i    n ulp max( |b(i) A|, |a(i) B| ) */
00079 
00080 /*  using the 1-norm.  Here, a(i)/b(i) = w is the i-th generalized */
00081 /*  eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th */
00082 /*  generalized eigenvalue of m A - B. */
00083 
00084 /*                          H   H  _      _ */
00085 /*  For left eigenvectors, A , B , a, and b  are used. */
00086 
00087 /*  ZGET52 also tests the normalization of E.  Each eigenvector is */
00088 /*  supposed to be normalized so that the maximum "absolute value" */
00089 /*  of its elements is 1, where in this case, "absolute value" */
00090 /*  of a complex value x is  |Re(x)| + |Im(x)| ; let us call this */
00091 /*  maximum "absolute value" norm of a vector v  M(v). */
00092 /*  If a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate */
00093 /*  vector. The normalization test is: */
00094 
00095 /*          RESULT(2) =      max       | M(v(i)) - 1 | / ( n ulp ) */
00096 /*                     eigenvectors v(i) */
00097 
00098 
00099 /*  Arguments */
00100 /*  ========= */
00101 
00102 /*  LEFT    (input) LOGICAL */
00103 /*          =.TRUE.:  The eigenvectors in the columns of E are assumed */
00104 /*                    to be *left* eigenvectors. */
00105 /*          =.FALSE.: The eigenvectors in the columns of E are assumed */
00106 /*                    to be *right* eigenvectors. */
00107 
00108 /*  N       (input) INTEGER */
00109 /*          The size of the matrices.  If it is zero, ZGET52 does */
00110 /*          nothing.  It must be at least zero. */
00111 
00112 /*  A       (input) COMPLEX*16 array, dimension (LDA, N) */
00113 /*          The matrix A. */
00114 
00115 /*  LDA     (input) INTEGER */
00116 /*          The leading dimension of A.  It must be at least 1 */
00117 /*          and at least N. */
00118 
00119 /*  B       (input) COMPLEX*16 array, dimension (LDB, N) */
00120 /*          The matrix B. */
00121 
00122 /*  LDB     (input) INTEGER */
00123 /*          The leading dimension of B.  It must be at least 1 */
00124 /*          and at least N. */
00125 
00126 /*  E       (input) COMPLEX*16 array, dimension (LDE, N) */
00127 /*          The matrix of eigenvectors.  It must be O( 1 ). */
00128 
00129 /*  LDE     (input) INTEGER */
00130 /*          The leading dimension of E.  It must be at least 1 and at */
00131 /*          least N. */
00132 
00133 /*  ALPHA   (input) COMPLEX*16 array, dimension (N) */
00134 /*          The values a(i) as described above, which, along with b(i), */
00135 /*          define the generalized eigenvalues. */
00136 
00137 /*  BETA    (input) COMPLEX*16 array, dimension (N) */
00138 /*          The values b(i) as described above, which, along with a(i), */
00139 /*          define the generalized eigenvalues. */
00140 
00141 /*  WORK    (workspace) COMPLEX*16 array, dimension (N**2) */
00142 
00143 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N) */
00144 
00145 /*  RESULT  (output) DOUBLE PRECISION array, dimension (2) */
00146 /*          The values computed by the test described above.  If A E or */
00147 /*          B E is likely to overflow, then RESULT(1:2) is set to */
00148 /*          10 / ulp. */
00149 
00150 /*  ===================================================================== */
00151 
00152 /*     .. Parameters .. */
00153 /*     .. */
00154 /*     .. Local Scalars .. */
00155 /*     .. */
00156 /*     .. External Functions .. */
00157 /*     .. */
00158 /*     .. External Subroutines .. */
00159 /*     .. */
00160 /*     .. Intrinsic Functions .. */
00161 /*     .. */
00162 /*     .. Statement Functions .. */
00163 /*     .. */
00164 /*     .. Statement Function definitions .. */
00165 /*     .. */
00166 /*     .. Executable Statements .. */
00167 
00168     /* Parameter adjustments */
00169     a_dim1 = *lda;
00170     a_offset = 1 + a_dim1;
00171     a -= a_offset;
00172     b_dim1 = *ldb;
00173     b_offset = 1 + b_dim1;
00174     b -= b_offset;
00175     e_dim1 = *lde;
00176     e_offset = 1 + e_dim1;
00177     e -= e_offset;
00178     --alpha;
00179     --beta;
00180     --work;
00181     --rwork;
00182     --result;
00183 
00184     /* Function Body */
00185     result[1] = 0.;
00186     result[2] = 0.;
00187     if (*n <= 0) {
00188         return 0;
00189     }
00190 
00191     safmin = dlamch_("Safe minimum");
00192     safmax = 1. / safmin;
00193     ulp = dlamch_("Epsilon") * dlamch_("Base");
00194 
00195     if (*left) {
00196         *(unsigned char *)trans = 'C';
00197         *(unsigned char *)normab = 'I';
00198     } else {
00199         *(unsigned char *)trans = 'N';
00200         *(unsigned char *)normab = 'O';
00201     }
00202 
00203 /*     Norm of A, B, and E: */
00204 
00205 /* Computing MAX */
00206     d__1 = zlange_(normab, n, n, &a[a_offset], lda, &rwork[1]);
00207     anorm = max(d__1,safmin);
00208 /* Computing MAX */
00209     d__1 = zlange_(normab, n, n, &b[b_offset], ldb, &rwork[1]);
00210     bnorm = max(d__1,safmin);
00211 /* Computing MAX */
00212     d__1 = zlange_("O", n, n, &e[e_offset], lde, &rwork[1]);
00213     enorm = max(d__1,ulp);
00214     alfmax = safmax / max(1.,bnorm);
00215     betmax = safmax / max(1.,anorm);
00216 
00217 /*     Compute error matrix. */
00218 /*     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| ) */
00219 
00220     i__1 = *n;
00221     for (jvec = 1; jvec <= i__1; ++jvec) {
00222         i__2 = jvec;
00223         alphai.r = alpha[i__2].r, alphai.i = alpha[i__2].i;
00224         i__2 = jvec;
00225         betai.r = beta[i__2].r, betai.i = beta[i__2].i;
00226 /* Computing MAX */
00227         d__5 = (d__1 = alphai.r, abs(d__1)) + (d__2 = d_imag(&alphai), abs(
00228                 d__2)), d__6 = (d__3 = betai.r, abs(d__3)) + (d__4 = d_imag(&
00229                 betai), abs(d__4));
00230         abmax = max(d__5,d__6);
00231         if ((d__1 = alphai.r, abs(d__1)) + (d__2 = d_imag(&alphai), abs(d__2))
00232                  > alfmax || (d__3 = betai.r, abs(d__3)) + (d__4 = d_imag(&
00233                 betai), abs(d__4)) > betmax || abmax < 1.) {
00234             scale = 1. / max(abmax,safmin);
00235             z__1.r = scale * alphai.r, z__1.i = scale * alphai.i;
00236             alphai.r = z__1.r, alphai.i = z__1.i;
00237             z__1.r = scale * betai.r, z__1.i = scale * betai.i;
00238             betai.r = z__1.r, betai.i = z__1.i;
00239         }
00240 /* Computing MAX */
00241         d__5 = ((d__1 = alphai.r, abs(d__1)) + (d__2 = d_imag(&alphai), abs(
00242                 d__2))) * bnorm, d__6 = ((d__3 = betai.r, abs(d__3)) + (d__4 =
00243                  d_imag(&betai), abs(d__4))) * anorm, d__5 = max(d__5,d__6);
00244         scale = 1. / max(d__5,safmin);
00245         z__1.r = scale * betai.r, z__1.i = scale * betai.i;
00246         acoeff.r = z__1.r, acoeff.i = z__1.i;
00247         z__1.r = scale * alphai.r, z__1.i = scale * alphai.i;
00248         bcoeff.r = z__1.r, bcoeff.i = z__1.i;
00249         if (*left) {
00250             d_cnjg(&z__1, &acoeff);
00251             acoeff.r = z__1.r, acoeff.i = z__1.i;
00252             d_cnjg(&z__1, &bcoeff);
00253             bcoeff.r = z__1.r, bcoeff.i = z__1.i;
00254         }
00255         zgemv_(trans, n, n, &acoeff, &a[a_offset], lda, &e[jvec * e_dim1 + 1], 
00256                  &c__1, &c_b1, &work[*n * (jvec - 1) + 1], &c__1);
00257         z__1.r = -bcoeff.r, z__1.i = -bcoeff.i;
00258         zgemv_(trans, n, n, &z__1, &b[b_offset], lda, &e[jvec * e_dim1 + 1], &
00259                 c__1, &c_b2, &work[*n * (jvec - 1) + 1], &c__1);
00260 /* L10: */
00261     }
00262 
00263     errnrm = zlange_("One", n, n, &work[1], n, &rwork[1]) / enorm;
00264 
00265 /*     Compute RESULT(1) */
00266 
00267     result[1] = errnrm / ulp;
00268 
00269 /*     Normalization of E: */
00270 
00271     enrmer = 0.;
00272     i__1 = *n;
00273     for (jvec = 1; jvec <= i__1; ++jvec) {
00274         temp1 = 0.;
00275         i__2 = *n;
00276         for (j = 1; j <= i__2; ++j) {
00277 /* Computing MAX */
00278             i__3 = j + jvec * e_dim1;
00279             d__3 = temp1, d__4 = (d__1 = e[i__3].r, abs(d__1)) + (d__2 = 
00280                     d_imag(&e[j + jvec * e_dim1]), abs(d__2));
00281             temp1 = max(d__3,d__4);
00282 /* L20: */
00283         }
00284 /* Computing MAX */
00285         d__1 = enrmer, d__2 = temp1 - 1.;
00286         enrmer = max(d__1,d__2);
00287 /* L30: */
00288     }
00289 
00290 /*     Compute RESULT(2) : the normalization error in E. */
00291 
00292     result[2] = enrmer / ((doublereal) (*n) * ulp);
00293 
00294     return 0;
00295 
00296 /*     End of ZGET52 */
00297 
00298 } /* zget52_ */

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