sget52.c
Go to the documentation of this file.
00001 /* sget52.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 integer c__1 = 1;
00019 static real c_b12 = 0.f;
00020 static real c_b15 = 1.f;
00021 
00022 /* Subroutine */ int sget52_(logical *left, integer *n, real *a, integer *lda, 
00023          real *b, integer *ldb, real *e, integer *lde, real *alphar, real *
00024         alphai, real *beta, real *work, 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     real r__1, r__2, r__3, r__4;
00029 
00030     /* Local variables */
00031     integer j;
00032     real ulp;
00033     integer jvec;
00034     real temp1, acoef, scale, abmax, salfi, sbeta, salfr, anorm, bnorm, enorm;
00035     extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 
00036             real *, integer *, real *, integer *, real *, real *, integer *);
00037     char trans[1];
00038     real bcoefi, bcoefr, alfmax;
00039     extern doublereal slamch_(char *), slange_(char *, integer *, 
00040             integer *, real *, integer *, real *);
00041     real safmin;
00042     char normab[1];
00043     real safmax, betmax, enrmer;
00044     logical ilcplx;
00045     real errnrm;
00046 
00047 
00048 /*  -- LAPACK test routine (version 3.1) -- */
00049 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00050 /*     November 2006 */
00051 
00052 /*     .. Scalar Arguments .. */
00053 /*     .. */
00054 /*     .. Array Arguments .. */
00055 /*     .. */
00056 
00057 /*  Purpose */
00058 /*  ======= */
00059 
00060 /*  SGET52  does an eigenvector check for the generalized eigenvalue */
00061 /*  problem. */
00062 
00063 /*  The basic test for right eigenvectors is: */
00064 
00065 /*                            | b(j) A E(j) -  a(j) B E(j) | */
00066 /*          RESULT(1) = max   ------------------------------- */
00067 /*                       j    n ulp max( |b(j) A|, |a(j) B| ) */
00068 
00069 /*  using the 1-norm.  Here, a(j)/b(j) = w is the j-th generalized */
00070 /*  eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th */
00071 /*  generalized eigenvalue of m A - B. */
00072 
00073 /*  For real eigenvalues, the test is straightforward.  For complex */
00074 /*  eigenvalues, E(j) and a(j) are complex, represented by */
00075 /*  Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that */
00076 /*  eigenvector becomes */
00077 
00078 /*                  max( |Wr|, |Wi| ) */
00079 /*      -------------------------------------------- */
00080 /*      n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| ) */
00081 
00082 /*  where */
00083 
00084 /*      Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j) */
00085 
00086 /*      Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j) */
00087 
00088 /*                          T   T  _ */
00089 /*  For left eigenvectors, A , B , a, and b  are used. */
00090 
00091 /*  SGET52 also tests the normalization of E.  Each eigenvector is */
00092 /*  supposed to be normalized so that the maximum "absolute value" */
00093 /*  of its elements is 1, where in this case, "absolute value" */
00094 /*  of a complex value x is  |Re(x)| + |Im(x)| ; let us call this */
00095 /*  maximum "absolute value" norm of a vector v  M(v). */
00096 /*  if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate */
00097 /*  vector.  The normalization test is: */
00098 
00099 /*          RESULT(2) =      max       | M(v(j)) - 1 | / ( n ulp ) */
00100 /*                     eigenvectors v(j) */
00101 
00102 /*  Arguments */
00103 /*  ========= */
00104 
00105 /*  LEFT    (input) LOGICAL */
00106 /*          =.TRUE.:  The eigenvectors in the columns of E are assumed */
00107 /*                    to be *left* eigenvectors. */
00108 /*          =.FALSE.: The eigenvectors in the columns of E are assumed */
00109 /*                    to be *right* eigenvectors. */
00110 
00111 /*  N       (input) INTEGER */
00112 /*          The size of the matrices.  If it is zero, SGET52 does */
00113 /*          nothing.  It must be at least zero. */
00114 
00115 /*  A       (input) REAL array, dimension (LDA, N) */
00116 /*          The matrix A. */
00117 
00118 /*  LDA     (input) INTEGER */
00119 /*          The leading dimension of A.  It must be at least 1 */
00120 /*          and at least N. */
00121 
00122 /*  B       (input) REAL array, dimension (LDB, N) */
00123 /*          The matrix B. */
00124 
00125 /*  LDB     (input) INTEGER */
00126 /*          The leading dimension of B.  It must be at least 1 */
00127 /*          and at least N. */
00128 
00129 /*  E       (input) REAL array, dimension (LDE, N) */
00130 /*          The matrix of eigenvectors.  It must be O( 1 ).  Complex */
00131 /*          eigenvalues and eigenvectors always come in pairs, the */
00132 /*          eigenvalue and its conjugate being stored in adjacent */
00133 /*          elements of ALPHAR, ALPHAI, and BETA.  Thus, if a(j)/b(j) */
00134 /*          and a(j+1)/b(j+1) are a complex conjugate pair of */
00135 /*          generalized eigenvalues, then E(,j) contains the real part */
00136 /*          of the eigenvector and E(,j+1) contains the imaginary part. */
00137 /*          Note that whether E(,j) is a real eigenvector or part of a */
00138 /*          complex one is specified by whether ALPHAI(j) is zero or not. */
00139 
00140 /*  LDE     (input) INTEGER */
00141 /*          The leading dimension of E.  It must be at least 1 and at */
00142 /*          least N. */
00143 
00144 /*  ALPHAR  (input) REAL array, dimension (N) */
00145 /*          The real parts of the values a(j) as described above, which, */
00146 /*          along with b(j), define the generalized eigenvalues. */
00147 /*          Complex eigenvalues always come in complex conjugate pairs */
00148 /*          a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent */
00149 /*          elements in ALPHAR, ALPHAI, and BETA.  Thus, if the j-th */
00150 /*          and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1) */
00151 /*          is assumed to be equal to ALPHAR(j)/BETA(j). */
00152 
00153 /*  ALPHAI  (input) REAL array, dimension (N) */
00154 /*          The imaginary parts of the values a(j) as described above, */
00155 /*          which, along with b(j), define the generalized eigenvalues. */
00156 /*          If ALPHAI(j)=0, then the eigenvalue is real, otherwise it */
00157 /*          is part of a complex conjugate pair.  Complex eigenvalues */
00158 /*          always come in complex conjugate pairs a(j)/b(j) and */
00159 /*          a(j+1)/b(j+1), which are stored in adjacent elements in */
00160 /*          ALPHAR, ALPHAI, and BETA.  Thus, if the j-th and (j+1)-st */
00161 /*          eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to */
00162 /*          be equal to  -ALPHAI(j)/BETA(j).  Also, nonzero values in */
00163 /*          ALPHAI are assumed to always come in adjacent pairs. */
00164 
00165 /*  BETA    (input) REAL array, dimension (N) */
00166 /*          The values b(j) as described above, which, along with a(j), */
00167 /*          define the generalized eigenvalues. */
00168 
00169 /*  WORK    (workspace) REAL array, dimension (N**2+N) */
00170 
00171 /*  RESULT  (output) REAL array, dimension (2) */
00172 /*          The values computed by the test described above.  If A E or */
00173 /*          B E is likely to overflow, then RESULT(1:2) is set to */
00174 /*          10 / ulp. */
00175 
00176 /*  ===================================================================== */
00177 
00178 /*     .. Parameters .. */
00179 /*     .. */
00180 /*     .. Local Scalars .. */
00181 /*     .. */
00182 /*     .. External Functions .. */
00183 /*     .. */
00184 /*     .. External Subroutines .. */
00185 /*     .. */
00186 /*     .. Intrinsic Functions .. */
00187 /*     .. */
00188 /*     .. Executable Statements .. */
00189 
00190     /* Parameter adjustments */
00191     a_dim1 = *lda;
00192     a_offset = 1 + a_dim1;
00193     a -= a_offset;
00194     b_dim1 = *ldb;
00195     b_offset = 1 + b_dim1;
00196     b -= b_offset;
00197     e_dim1 = *lde;
00198     e_offset = 1 + e_dim1;
00199     e -= e_offset;
00200     --alphar;
00201     --alphai;
00202     --beta;
00203     --work;
00204     --result;
00205 
00206     /* Function Body */
00207     result[1] = 0.f;
00208     result[2] = 0.f;
00209     if (*n <= 0) {
00210         return 0;
00211     }
00212 
00213     safmin = slamch_("Safe minimum");
00214     safmax = 1.f / safmin;
00215     ulp = slamch_("Epsilon") * slamch_("Base");
00216 
00217     if (*left) {
00218         *(unsigned char *)trans = 'T';
00219         *(unsigned char *)normab = 'I';
00220     } else {
00221         *(unsigned char *)trans = 'N';
00222         *(unsigned char *)normab = 'O';
00223     }
00224 
00225 /*     Norm of A, B, and E: */
00226 
00227 /* Computing MAX */
00228     r__1 = slange_(normab, n, n, &a[a_offset], lda, &work[1]);
00229     anorm = dmax(r__1,safmin);
00230 /* Computing MAX */
00231     r__1 = slange_(normab, n, n, &b[b_offset], ldb, &work[1]);
00232     bnorm = dmax(r__1,safmin);
00233 /* Computing MAX */
00234     r__1 = slange_("O", n, n, &e[e_offset], lde, &work[1]);
00235     enorm = dmax(r__1,ulp);
00236     alfmax = safmax / dmax(1.f,bnorm);
00237     betmax = safmax / dmax(1.f,anorm);
00238 
00239 /*     Compute error matrix. */
00240 /*     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| ) */
00241 
00242     ilcplx = FALSE_;
00243     i__1 = *n;
00244     for (jvec = 1; jvec <= i__1; ++jvec) {
00245         if (ilcplx) {
00246 
00247 /*           2nd Eigenvalue/-vector of pair -- do nothing */
00248 
00249             ilcplx = FALSE_;
00250         } else {
00251             salfr = alphar[jvec];
00252             salfi = alphai[jvec];
00253             sbeta = beta[jvec];
00254             if (salfi == 0.f) {
00255 
00256 /*              Real eigenvalue and -vector */
00257 
00258 /* Computing MAX */
00259                 r__1 = dabs(salfr), r__2 = dabs(sbeta);
00260                 abmax = dmax(r__1,r__2);
00261                 if (dabs(salfr) > alfmax || dabs(sbeta) > betmax || abmax < 
00262                         1.f) {
00263                     scale = 1.f / dmax(abmax,safmin);
00264                     salfr = scale * salfr;
00265                     sbeta = scale * sbeta;
00266                 }
00267 /* Computing MAX */
00268                 r__1 = dabs(salfr) * bnorm, r__2 = dabs(sbeta) * anorm, r__1 =
00269                          max(r__1,r__2);
00270                 scale = 1.f / dmax(r__1,safmin);
00271                 acoef = scale * sbeta;
00272                 bcoefr = scale * salfr;
00273                 sgemv_(trans, n, n, &acoef, &a[a_offset], lda, &e[jvec * 
00274                         e_dim1 + 1], &c__1, &c_b12, &work[*n * (jvec - 1) + 1]
00275 , &c__1);
00276                 r__1 = -bcoefr;
00277                 sgemv_(trans, n, n, &r__1, &b[b_offset], lda, &e[jvec * 
00278                         e_dim1 + 1], &c__1, &c_b15, &work[*n * (jvec - 1) + 1]
00279 , &c__1);
00280             } else {
00281 
00282 /*              Complex conjugate pair */
00283 
00284                 ilcplx = TRUE_;
00285                 if (jvec == *n) {
00286                     result[1] = 10.f / ulp;
00287                     return 0;
00288                 }
00289 /* Computing MAX */
00290                 r__1 = dabs(salfr) + dabs(salfi), r__2 = dabs(sbeta);
00291                 abmax = dmax(r__1,r__2);
00292                 if (dabs(salfr) + dabs(salfi) > alfmax || dabs(sbeta) > 
00293                         betmax || abmax < 1.f) {
00294                     scale = 1.f / dmax(abmax,safmin);
00295                     salfr = scale * salfr;
00296                     salfi = scale * salfi;
00297                     sbeta = scale * sbeta;
00298                 }
00299 /* Computing MAX */
00300                 r__1 = (dabs(salfr) + dabs(salfi)) * bnorm, r__2 = dabs(sbeta)
00301                          * anorm, r__1 = max(r__1,r__2);
00302                 scale = 1.f / dmax(r__1,safmin);
00303                 acoef = scale * sbeta;
00304                 bcoefr = scale * salfr;
00305                 bcoefi = scale * salfi;
00306                 if (*left) {
00307                     bcoefi = -bcoefi;
00308                 }
00309 
00310                 sgemv_(trans, n, n, &acoef, &a[a_offset], lda, &e[jvec * 
00311                         e_dim1 + 1], &c__1, &c_b12, &work[*n * (jvec - 1) + 1]
00312 , &c__1);
00313                 r__1 = -bcoefr;
00314                 sgemv_(trans, n, n, &r__1, &b[b_offset], lda, &e[jvec * 
00315                         e_dim1 + 1], &c__1, &c_b15, &work[*n * (jvec - 1) + 1]
00316 , &c__1);
00317                 sgemv_(trans, n, n, &bcoefi, &b[b_offset], lda, &e[(jvec + 1) 
00318                         * e_dim1 + 1], &c__1, &c_b15, &work[*n * (jvec - 1) + 
00319                         1], &c__1);
00320 
00321                 sgemv_(trans, n, n, &acoef, &a[a_offset], lda, &e[(jvec + 1) *
00322                          e_dim1 + 1], &c__1, &c_b12, &work[*n * jvec + 1], &
00323                         c__1);
00324                 r__1 = -bcoefi;
00325                 sgemv_(trans, n, n, &r__1, &b[b_offset], lda, &e[jvec * 
00326                         e_dim1 + 1], &c__1, &c_b15, &work[*n * jvec + 1], &
00327                         c__1);
00328                 r__1 = -bcoefr;
00329                 sgemv_(trans, n, n, &r__1, &b[b_offset], lda, &e[(jvec + 1) * 
00330                         e_dim1 + 1], &c__1, &c_b15, &work[*n * jvec + 1], &
00331                         c__1);
00332             }
00333         }
00334 /* L10: */
00335     }
00336 
00337 /* Computing 2nd power */
00338     i__1 = *n;
00339     errnrm = slange_("One", n, n, &work[1], n, &work[i__1 * i__1 + 1]) / enorm;
00340 
00341 /*     Compute RESULT(1) */
00342 
00343     result[1] = errnrm / ulp;
00344 
00345 /*     Normalization of E: */
00346 
00347     enrmer = 0.f;
00348     ilcplx = FALSE_;
00349     i__1 = *n;
00350     for (jvec = 1; jvec <= i__1; ++jvec) {
00351         if (ilcplx) {
00352             ilcplx = FALSE_;
00353         } else {
00354             temp1 = 0.f;
00355             if (alphai[jvec] == 0.f) {
00356                 i__2 = *n;
00357                 for (j = 1; j <= i__2; ++j) {
00358 /* Computing MAX */
00359                     r__2 = temp1, r__3 = (r__1 = e[j + jvec * e_dim1], dabs(
00360                             r__1));
00361                     temp1 = dmax(r__2,r__3);
00362 /* L20: */
00363                 }
00364 /* Computing MAX */
00365                 r__1 = enrmer, r__2 = temp1 - 1.f;
00366                 enrmer = dmax(r__1,r__2);
00367             } else {
00368                 ilcplx = TRUE_;
00369                 i__2 = *n;
00370                 for (j = 1; j <= i__2; ++j) {
00371 /* Computing MAX */
00372                     r__3 = temp1, r__4 = (r__1 = e[j + jvec * e_dim1], dabs(
00373                             r__1)) + (r__2 = e[j + (jvec + 1) * e_dim1], dabs(
00374                             r__2));
00375                     temp1 = dmax(r__3,r__4);
00376 /* L30: */
00377                 }
00378 /* Computing MAX */
00379                 r__1 = enrmer, r__2 = temp1 - 1.f;
00380                 enrmer = dmax(r__1,r__2);
00381             }
00382         }
00383 /* L40: */
00384     }
00385 
00386 /*     Compute RESULT(2) : the normalization error in E. */
00387 
00388     result[2] = enrmer / ((real) (*n) * ulp);
00389 
00390     return 0;
00391 
00392 /*     End of SGET52 */
00393 
00394 } /* sget52_ */

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