swiftnav: cget52.c Source File

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