sget23.c
Go to the documentation of this file.
00001 /* sget23.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 integer c__4 = 4;
00020 
00021 /* Subroutine */ int sget23_(logical *comp, char *balanc, integer *jtype, 
00022         real *thresh, integer *iseed, integer *nounit, integer *n, real *a, 
00023         integer *lda, real *h__, real *wr, real *wi, real *wr1, real *wi1, 
00024         real *vl, integer *ldvl, real *vr, integer *ldvr, real *lre, integer *
00025         ldlre, real *rcondv, real *rcndv1, real *rcdvin, real *rconde, real *
00026         rcnde1, real *rcdein, real *scale, real *scale1, real *result, real *
00027         work, integer *lwork, integer *iwork, integer *info)
00028 {
00029     /* Initialized data */
00030 
00031     static char sens[1*2] = "N" "V";
00032 
00033     /* Format strings */
00034     static char fmt_9998[] = "(\002 SGET23: \002,a,\002 returned INFO=\002,i"
00035             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, BALANC = "
00036             "\002,a,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
00037     static char fmt_9999[] = "(\002 SGET23: \002,a,\002 returned INFO=\002,i"
00038             "6,\002.\002,/9x,\002N=\002,i6,\002, INPUT EXAMPLE NUMBER = \002,"
00039             "i4)";
00040 
00041     /* System generated locals */
00042     integer a_dim1, a_offset, h_dim1, h_offset, lre_dim1, lre_offset, vl_dim1,
00043              vl_offset, vr_dim1, vr_offset, i__1, i__2, i__3;
00044     real r__1, r__2, r__3, r__4, r__5;
00045 
00046     /* Builtin functions */
00047     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00048 
00049     /* Local variables */
00050     integer i__, j;
00051     real v;
00052     integer jj, ihi, ilo;
00053     real dum[1], eps, res[2], tol, ulp, vmx;
00054     integer ihi1, ilo1, kmin;
00055     real vmax, tnrm, vrmx, vtst;
00056     extern doublereal snrm2_(integer *, real *, integer *);
00057     logical balok, nobal;
00058     real abnrm;
00059     extern logical lsame_(char *, char *);
00060     integer iinfo;
00061     extern /* Subroutine */ int sget22_(char *, char *, char *, integer *, 
00062             real *, integer *, real *, integer *, real *, real *, real *, 
00063             real *);
00064     char sense[1];
00065     integer isens;
00066     real vimin, tolin, vrmin, abnrm1;
00067     extern doublereal slapy2_(real *, real *), slamch_(char *);
00068     extern /* Subroutine */ int xerbla_(char *, integer *), slacpy_(
00069             char *, integer *, integer *, real *, integer *, real *, integer *
00070 );
00071     integer isensm;
00072     extern /* Subroutine */ int sgeevx_(char *, char *, char *, char *, 
00073             integer *, real *, integer *, real *, real *, real *, integer *, 
00074             real *, integer *, integer *, integer *, real *, real *, real *, 
00075             real *, real *, integer *, integer *, integer *);
00076     real smlnum, ulpinv;
00077 
00078     /* Fortran I/O blocks */
00079     static cilist io___14 = { 0, 0, 0, fmt_9998, 0 };
00080     static cilist io___15 = { 0, 0, 0, fmt_9999, 0 };
00081     static cilist io___28 = { 0, 0, 0, fmt_9998, 0 };
00082     static cilist io___29 = { 0, 0, 0, fmt_9999, 0 };
00083     static cilist io___30 = { 0, 0, 0, fmt_9998, 0 };
00084     static cilist io___31 = { 0, 0, 0, fmt_9999, 0 };
00085     static cilist io___32 = { 0, 0, 0, fmt_9998, 0 };
00086     static cilist io___33 = { 0, 0, 0, fmt_9999, 0 };
00087     static cilist io___34 = { 0, 0, 0, fmt_9999, 0 };
00088 
00089 
00090 
00091 /*  -- LAPACK test routine (version 3.1) -- */
00092 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00093 /*     November 2006 */
00094 
00095 /*     .. Scalar Arguments .. */
00096 /*     .. */
00097 /*     .. Array Arguments .. */
00098 /*     .. */
00099 
00100 /*  Purpose */
00101 /*  ======= */
00102 
00103 /*     SGET23  checks the nonsymmetric eigenvalue problem driver SGEEVX. */
00104 /*     If COMP = .FALSE., the first 8 of the following tests will be */
00105 /*     performed on the input matrix A, and also test 9 if LWORK is */
00106 /*     sufficiently large. */
00107 /*     if COMP is .TRUE. all 11 tests will be performed. */
00108 
00109 /*     (1)     | A * VR - VR * W | / ( n |A| ulp ) */
00110 
00111 /*       Here VR is the matrix of unit right eigenvectors. */
00112 /*       W is a block diagonal matrix, with a 1x1 block for each */
00113 /*       real eigenvalue and a 2x2 block for each complex conjugate */
00114 /*       pair.  If eigenvalues j and j+1 are a complex conjugate pair, */
00115 /*       so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the */
00116 /*       2 x 2 block corresponding to the pair will be: */
00117 
00118 /*               (  wr  wi  ) */
00119 /*               ( -wi  wr  ) */
00120 
00121 /*       Such a block multiplying an n x 2 matrix  ( ur ui ) on the */
00122 /*       right will be the same as multiplying  ur + i*ui  by  wr + i*wi. */
00123 
00124 /*     (2)     | A**H * VL - VL * W**H | / ( n |A| ulp ) */
00125 
00126 /*       Here VL is the matrix of unit left eigenvectors, A**H is the */
00127 /*       conjugate transpose of A, and W is as above. */
00128 
00129 /*     (3)     | |VR(i)| - 1 | / ulp and largest component real */
00130 
00131 /*       VR(i) denotes the i-th column of VR. */
00132 
00133 /*     (4)     | |VL(i)| - 1 | / ulp and largest component real */
00134 
00135 /*       VL(i) denotes the i-th column of VL. */
00136 
00137 /*     (5)     0 if W(full) = W(partial), 1/ulp otherwise */
00138 
00139 /*       W(full) denotes the eigenvalues computed when VR, VL, RCONDV */
00140 /*       and RCONDE are also computed, and W(partial) denotes the */
00141 /*       eigenvalues computed when only some of VR, VL, RCONDV, and */
00142 /*       RCONDE are computed. */
00143 
00144 /*     (6)     0 if VR(full) = VR(partial), 1/ulp otherwise */
00145 
00146 /*       VR(full) denotes the right eigenvectors computed when VL, RCONDV */
00147 /*       and RCONDE are computed, and VR(partial) denotes the result */
00148 /*       when only some of VL and RCONDV are computed. */
00149 
00150 /*     (7)     0 if VL(full) = VL(partial), 1/ulp otherwise */
00151 
00152 /*       VL(full) denotes the left eigenvectors computed when VR, RCONDV */
00153 /*       and RCONDE are computed, and VL(partial) denotes the result */
00154 /*       when only some of VR and RCONDV are computed. */
00155 
00156 /*     (8)     0 if SCALE, ILO, IHI, ABNRM (full) = */
00157 /*                  SCALE, ILO, IHI, ABNRM (partial) */
00158 /*             1/ulp otherwise */
00159 
00160 /*       SCALE, ILO, IHI and ABNRM describe how the matrix is balanced. */
00161 /*       (full) is when VR, VL, RCONDE and RCONDV are also computed, and */
00162 /*       (partial) is when some are not computed. */
00163 
00164 /*     (9)     0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise */
00165 
00166 /*       RCONDV(full) denotes the reciprocal condition numbers of the */
00167 /*       right eigenvectors computed when VR, VL and RCONDE are also */
00168 /*       computed. RCONDV(partial) denotes the reciprocal condition */
00169 /*       numbers when only some of VR, VL and RCONDE are computed. */
00170 
00171 /*    (10)     |RCONDV - RCDVIN| / cond(RCONDV) */
00172 
00173 /*       RCONDV is the reciprocal right eigenvector condition number */
00174 /*       computed by SGEEVX and RCDVIN (the precomputed true value) */
00175 /*       is supplied as input. cond(RCONDV) is the condition number of */
00176 /*       RCONDV, and takes errors in computing RCONDV into account, so */
00177 /*       that the resulting quantity should be O(ULP). cond(RCONDV) is */
00178 /*       essentially given by norm(A)/RCONDE. */
00179 
00180 /*    (11)     |RCONDE - RCDEIN| / cond(RCONDE) */
00181 
00182 /*       RCONDE is the reciprocal eigenvalue condition number */
00183 /*       computed by SGEEVX and RCDEIN (the precomputed true value) */
00184 /*       is supplied as input.  cond(RCONDE) is the condition number */
00185 /*       of RCONDE, and takes errors in computing RCONDE into account, */
00186 /*       so that the resulting quantity should be O(ULP). cond(RCONDE) */
00187 /*       is essentially given by norm(A)/RCONDV. */
00188 
00189 /*  Arguments */
00190 /*  ========= */
00191 
00192 /*  COMP    (input) LOGICAL */
00193 /*          COMP describes which input tests to perform: */
00194 /*            = .FALSE. if the computed condition numbers are not to */
00195 /*                      be tested against RCDVIN and RCDEIN */
00196 /*            = .TRUE.  if they are to be compared */
00197 
00198 /*  BALANC  (input) CHARACTER */
00199 /*          Describes the balancing option to be tested. */
00200 /*            = 'N' for no permuting or diagonal scaling */
00201 /*            = 'P' for permuting but no diagonal scaling */
00202 /*            = 'S' for no permuting but diagonal scaling */
00203 /*            = 'B' for permuting and diagonal scaling */
00204 
00205 /*  JTYPE   (input) INTEGER */
00206 /*          Type of input matrix. Used to label output if error occurs. */
00207 
00208 /*  THRESH  (input) REAL */
00209 /*          A test will count as "failed" if the "error", computed as */
00210 /*          described above, exceeds THRESH.  Note that the error */
00211 /*          is scaled to be O(1), so THRESH should be a reasonably */
00212 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00213 /*          it should not depend on the precision (single vs. double) */
00214 /*          or the size of the matrix.  It must be at least zero. */
00215 
00216 /*  ISEED   (input) INTEGER array, dimension (4) */
00217 /*          If COMP = .FALSE., the random number generator seed */
00218 /*          used to produce matrix. */
00219 /*          If COMP = .TRUE., ISEED(1) = the number of the example. */
00220 /*          Used to label output if error occurs. */
00221 
00222 /*  NOUNIT  (input) INTEGER */
00223 /*          The FORTRAN unit number for printing out error messages */
00224 /*          (e.g., if a routine returns INFO not equal to 0.) */
00225 
00226 /*  N       (input) INTEGER */
00227 /*          The dimension of A. N must be at least 0. */
00228 
00229 /*  A       (input/output) REAL array, dimension (LDA,N) */
00230 /*          Used to hold the matrix whose eigenvalues are to be */
00231 /*          computed. */
00232 
00233 /*  LDA     (input) INTEGER */
00234 /*          The leading dimension of A, and H. LDA must be at */
00235 /*          least 1 and at least N. */
00236 
00237 /*  H       (workspace) REAL array, dimension (LDA,N) */
00238 /*          Another copy of the test matrix A, modified by SGEEVX. */
00239 
00240 /*  WR      (workspace) REAL array, dimension (N) */
00241 /*  WI      (workspace) REAL array, dimension (N) */
00242 /*          The real and imaginary parts of the eigenvalues of A. */
00243 /*          On exit, WR + WI*i are the eigenvalues of the matrix in A. */
00244 
00245 /*  WR1     (workspace) REAL array, dimension (N) */
00246 /*  WI1     (workspace) REAL array, dimension (N) */
00247 /*          Like WR, WI, these arrays contain the eigenvalues of A, */
00248 /*          but those computed when SGEEVX only computes a partial */
00249 /*          eigendecomposition, i.e. not the eigenvalues and left */
00250 /*          and right eigenvectors. */
00251 
00252 /*  VL      (workspace) REAL array, dimension (LDVL,N) */
00253 /*          VL holds the computed left eigenvectors. */
00254 
00255 /*  LDVL    (input) INTEGER */
00256 /*          Leading dimension of VL. Must be at least max(1,N). */
00257 
00258 /*  VR      (workspace) REAL array, dimension (LDVR,N) */
00259 /*          VR holds the computed right eigenvectors. */
00260 
00261 /*  LDVR    (input) INTEGER */
00262 /*          Leading dimension of VR. Must be at least max(1,N). */
00263 
00264 /*  LRE     (workspace) REAL array, dimension (LDLRE,N) */
00265 /*          LRE holds the computed right or left eigenvectors. */
00266 
00267 /*  LDLRE   (input) INTEGER */
00268 /*          Leading dimension of LRE. Must be at least max(1,N). */
00269 
00270 /*  RCONDV  (workspace) REAL array, dimension (N) */
00271 /*          RCONDV holds the computed reciprocal condition numbers */
00272 /*          for eigenvectors. */
00273 
00274 /*  RCNDV1  (workspace) REAL array, dimension (N) */
00275 /*          RCNDV1 holds more computed reciprocal condition numbers */
00276 /*          for eigenvectors. */
00277 
00278 /*  RCDVIN  (input) REAL array, dimension (N) */
00279 /*          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal */
00280 /*          condition numbers for eigenvectors to be compared with */
00281 /*          RCONDV. */
00282 
00283 /*  RCONDE  (workspace) REAL array, dimension (N) */
00284 /*          RCONDE holds the computed reciprocal condition numbers */
00285 /*          for eigenvalues. */
00286 
00287 /*  RCNDE1  (workspace) REAL array, dimension (N) */
00288 /*          RCNDE1 holds more computed reciprocal condition numbers */
00289 /*          for eigenvalues. */
00290 
00291 /*  RCDEIN  (input) REAL array, dimension (N) */
00292 /*          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal */
00293 /*          condition numbers for eigenvalues to be compared with */
00294 /*          RCONDE. */
00295 
00296 /*  SCALE   (workspace) REAL array, dimension (N) */
00297 /*          Holds information describing balancing of matrix. */
00298 
00299 /*  SCALE1  (workspace) REAL array, dimension (N) */
00300 /*          Holds information describing balancing of matrix. */
00301 
00302 /*  RESULT  (output) REAL array, dimension (11) */
00303 /*          The values computed by the 11 tests described above. */
00304 /*          The values are currently limited to 1/ulp, to avoid */
00305 /*          overflow. */
00306 
00307 /*  WORK    (workspace) REAL array, dimension (LWORK) */
00308 
00309 /*  LWORK   (input) INTEGER */
00310 /*          The number of entries in WORK.  This must be at least */
00311 /*          3*N, and 6*N+N**2 if tests 9, 10 or 11 are to be performed. */
00312 
00313 /*  IWORK   (workspace) INTEGER array, dimension (2*N) */
00314 
00315 /*  INFO    (output) INTEGER */
00316 /*          If 0,  successful exit. */
00317 /*          If <0, input parameter -INFO had an incorrect value. */
00318 /*          If >0, SGEEVX returned an error code, the absolute */
00319 /*                 value of which is returned. */
00320 
00321 /*  ===================================================================== */
00322 
00323 
00324 /*     .. Parameters .. */
00325 /*     .. */
00326 /*     .. Local Scalars .. */
00327 /*     .. */
00328 /*     .. Local Arrays .. */
00329 /*     .. */
00330 /*     .. External Functions .. */
00331 /*     .. */
00332 /*     .. External Subroutines .. */
00333 /*     .. */
00334 /*     .. Intrinsic Functions .. */
00335 /*     .. */
00336 /*     .. Data statements .. */
00337     /* Parameter adjustments */
00338     --iseed;
00339     h_dim1 = *lda;
00340     h_offset = 1 + h_dim1;
00341     h__ -= h_offset;
00342     a_dim1 = *lda;
00343     a_offset = 1 + a_dim1;
00344     a -= a_offset;
00345     --wr;
00346     --wi;
00347     --wr1;
00348     --wi1;
00349     vl_dim1 = *ldvl;
00350     vl_offset = 1 + vl_dim1;
00351     vl -= vl_offset;
00352     vr_dim1 = *ldvr;
00353     vr_offset = 1 + vr_dim1;
00354     vr -= vr_offset;
00355     lre_dim1 = *ldlre;
00356     lre_offset = 1 + lre_dim1;
00357     lre -= lre_offset;
00358     --rcondv;
00359     --rcndv1;
00360     --rcdvin;
00361     --rconde;
00362     --rcnde1;
00363     --rcdein;
00364     --scale;
00365     --scale1;
00366     --result;
00367     --work;
00368     --iwork;
00369 
00370     /* Function Body */
00371 /*     .. */
00372 /*     .. Executable Statements .. */
00373 
00374 /*     Check for errors */
00375 
00376     nobal = lsame_(balanc, "N");
00377     balok = nobal || lsame_(balanc, "P") || lsame_(
00378             balanc, "S") || lsame_(balanc, "B");
00379     *info = 0;
00380     if (! balok) {
00381         *info = -2;
00382     } else if (*thresh < 0.f) {
00383         *info = -4;
00384     } else if (*nounit <= 0) {
00385         *info = -6;
00386     } else if (*n < 0) {
00387         *info = -7;
00388     } else if (*lda < 1 || *lda < *n) {
00389         *info = -9;
00390     } else if (*ldvl < 1 || *ldvl < *n) {
00391         *info = -16;
00392     } else if (*ldvr < 1 || *ldvr < *n) {
00393         *info = -18;
00394     } else if (*ldlre < 1 || *ldlre < *n) {
00395         *info = -20;
00396     } else if (*lwork < *n * 3 || *comp && *lwork < *n * 6 + *n * *n) {
00397         *info = -31;
00398     }
00399 
00400     if (*info != 0) {
00401         i__1 = -(*info);
00402         xerbla_("SGET23", &i__1);
00403         return 0;
00404     }
00405 
00406 /*     Quick return if nothing to do */
00407 
00408     for (i__ = 1; i__ <= 11; ++i__) {
00409         result[i__] = -1.f;
00410 /* L10: */
00411     }
00412 
00413     if (*n == 0) {
00414         return 0;
00415     }
00416 
00417 /*     More Important constants */
00418 
00419     ulp = slamch_("Precision");
00420     smlnum = slamch_("S");
00421     ulpinv = 1.f / ulp;
00422 
00423 /*     Compute eigenvalues and eigenvectors, and test them */
00424 
00425     if (*lwork >= *n * 6 + *n * *n) {
00426         *(unsigned char *)sense = 'B';
00427         isensm = 2;
00428     } else {
00429         *(unsigned char *)sense = 'E';
00430         isensm = 1;
00431     }
00432     slacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00433     sgeevx_(balanc, "V", "V", sense, n, &h__[h_offset], lda, &wr[1], &wi[1], &
00434             vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &ilo, &ihi, &scale[1], 
00435             &abnrm, &rconde[1], &rcondv[1], &work[1], lwork, &iwork[1], &
00436             iinfo);
00437     if (iinfo != 0) {
00438         result[1] = ulpinv;
00439         if (*jtype != 22) {
00440             io___14.ciunit = *nounit;
00441             s_wsfe(&io___14);
00442             do_fio(&c__1, "SGEEVX1", (ftnlen)7);
00443             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00444             do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00445             do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00446             do_fio(&c__1, balanc, (ftnlen)1);
00447             do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00448             e_wsfe();
00449         } else {
00450             io___15.ciunit = *nounit;
00451             s_wsfe(&io___15);
00452             do_fio(&c__1, "SGEEVX1", (ftnlen)7);
00453             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00454             do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00455             do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00456             e_wsfe();
00457         }
00458         *info = abs(iinfo);
00459         return 0;
00460     }
00461 
00462 /*     Do Test (1) */
00463 
00464     sget22_("N", "N", "N", n, &a[a_offset], lda, &vr[vr_offset], ldvr, &wr[1], 
00465              &wi[1], &work[1], res);
00466     result[1] = res[0];
00467 
00468 /*     Do Test (2) */
00469 
00470     sget22_("T", "N", "T", n, &a[a_offset], lda, &vl[vl_offset], ldvl, &wr[1], 
00471              &wi[1], &work[1], res);
00472     result[2] = res[0];
00473 
00474 /*     Do Test (3) */
00475 
00476     i__1 = *n;
00477     for (j = 1; j <= i__1; ++j) {
00478         tnrm = 1.f;
00479         if (wi[j] == 0.f) {
00480             tnrm = snrm2_(n, &vr[j * vr_dim1 + 1], &c__1);
00481         } else if (wi[j] > 0.f) {
00482             r__1 = snrm2_(n, &vr[j * vr_dim1 + 1], &c__1);
00483             r__2 = snrm2_(n, &vr[(j + 1) * vr_dim1 + 1], &c__1);
00484             tnrm = slapy2_(&r__1, &r__2);
00485         }
00486 /* Computing MAX */
00487 /* Computing MIN */
00488         r__4 = ulpinv, r__5 = (r__1 = tnrm - 1.f, dabs(r__1)) / ulp;
00489         r__2 = result[3], r__3 = dmin(r__4,r__5);
00490         result[3] = dmax(r__2,r__3);
00491         if (wi[j] > 0.f) {
00492             vmx = 0.f;
00493             vrmx = 0.f;
00494             i__2 = *n;
00495             for (jj = 1; jj <= i__2; ++jj) {
00496                 vtst = slapy2_(&vr[jj + j * vr_dim1], &vr[jj + (j + 1) * 
00497                         vr_dim1]);
00498                 if (vtst > vmx) {
00499                     vmx = vtst;
00500                 }
00501                 if (vr[jj + (j + 1) * vr_dim1] == 0.f && (r__1 = vr[jj + j * 
00502                         vr_dim1], dabs(r__1)) > vrmx) {
00503                     vrmx = (r__2 = vr[jj + j * vr_dim1], dabs(r__2));
00504                 }
00505 /* L20: */
00506             }
00507             if (vrmx / vmx < 1.f - ulp * 2.f) {
00508                 result[3] = ulpinv;
00509             }
00510         }
00511 /* L30: */
00512     }
00513 
00514 /*     Do Test (4) */
00515 
00516     i__1 = *n;
00517     for (j = 1; j <= i__1; ++j) {
00518         tnrm = 1.f;
00519         if (wi[j] == 0.f) {
00520             tnrm = snrm2_(n, &vl[j * vl_dim1 + 1], &c__1);
00521         } else if (wi[j] > 0.f) {
00522             r__1 = snrm2_(n, &vl[j * vl_dim1 + 1], &c__1);
00523             r__2 = snrm2_(n, &vl[(j + 1) * vl_dim1 + 1], &c__1);
00524             tnrm = slapy2_(&r__1, &r__2);
00525         }
00526 /* Computing MAX */
00527 /* Computing MIN */
00528         r__4 = ulpinv, r__5 = (r__1 = tnrm - 1.f, dabs(r__1)) / ulp;
00529         r__2 = result[4], r__3 = dmin(r__4,r__5);
00530         result[4] = dmax(r__2,r__3);
00531         if (wi[j] > 0.f) {
00532             vmx = 0.f;
00533             vrmx = 0.f;
00534             i__2 = *n;
00535             for (jj = 1; jj <= i__2; ++jj) {
00536                 vtst = slapy2_(&vl[jj + j * vl_dim1], &vl[jj + (j + 1) * 
00537                         vl_dim1]);
00538                 if (vtst > vmx) {
00539                     vmx = vtst;
00540                 }
00541                 if (vl[jj + (j + 1) * vl_dim1] == 0.f && (r__1 = vl[jj + j * 
00542                         vl_dim1], dabs(r__1)) > vrmx) {
00543                     vrmx = (r__2 = vl[jj + j * vl_dim1], dabs(r__2));
00544                 }
00545 /* L40: */
00546             }
00547             if (vrmx / vmx < 1.f - ulp * 2.f) {
00548                 result[4] = ulpinv;
00549             }
00550         }
00551 /* L50: */
00552     }
00553 
00554 /*     Test for all options of computing condition numbers */
00555 
00556     i__1 = isensm;
00557     for (isens = 1; isens <= i__1; ++isens) {
00558 
00559         *(unsigned char *)sense = *(unsigned char *)&sens[isens - 1];
00560 
00561 /*        Compute eigenvalues only, and test them */
00562 
00563         slacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00564         sgeevx_(balanc, "N", "N", sense, n, &h__[h_offset], lda, &wr1[1], &
00565                 wi1[1], dum, &c__1, dum, &c__1, &ilo1, &ihi1, &scale1[1], &
00566                 abnrm1, &rcnde1[1], &rcndv1[1], &work[1], lwork, &iwork[1], &
00567                 iinfo);
00568         if (iinfo != 0) {
00569             result[1] = ulpinv;
00570             if (*jtype != 22) {
00571                 io___28.ciunit = *nounit;
00572                 s_wsfe(&io___28);
00573                 do_fio(&c__1, "SGEEVX2", (ftnlen)7);
00574                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00575                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00576                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00577                 do_fio(&c__1, balanc, (ftnlen)1);
00578                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00579                 e_wsfe();
00580             } else {
00581                 io___29.ciunit = *nounit;
00582                 s_wsfe(&io___29);
00583                 do_fio(&c__1, "SGEEVX2", (ftnlen)7);
00584                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00585                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00586                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00587                 e_wsfe();
00588             }
00589             *info = abs(iinfo);
00590             goto L190;
00591         }
00592 
00593 /*        Do Test (5) */
00594 
00595         i__2 = *n;
00596         for (j = 1; j <= i__2; ++j) {
00597             if (wr[j] != wr1[j] || wi[j] != wi1[j]) {
00598                 result[5] = ulpinv;
00599             }
00600 /* L60: */
00601         }
00602 
00603 /*        Do Test (8) */
00604 
00605         if (! nobal) {
00606             i__2 = *n;
00607             for (j = 1; j <= i__2; ++j) {
00608                 if (scale[j] != scale1[j]) {
00609                     result[8] = ulpinv;
00610                 }
00611 /* L70: */
00612             }
00613             if (ilo != ilo1) {
00614                 result[8] = ulpinv;
00615             }
00616             if (ihi != ihi1) {
00617                 result[8] = ulpinv;
00618             }
00619             if (abnrm != abnrm1) {
00620                 result[8] = ulpinv;
00621             }
00622         }
00623 
00624 /*        Do Test (9) */
00625 
00626         if (isens == 2 && *n > 1) {
00627             i__2 = *n;
00628             for (j = 1; j <= i__2; ++j) {
00629                 if (rcondv[j] != rcndv1[j]) {
00630                     result[9] = ulpinv;
00631                 }
00632 /* L80: */
00633             }
00634         }
00635 
00636 /*        Compute eigenvalues and right eigenvectors, and test them */
00637 
00638         slacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00639         sgeevx_(balanc, "N", "V", sense, n, &h__[h_offset], lda, &wr1[1], &
00640                 wi1[1], dum, &c__1, &lre[lre_offset], ldlre, &ilo1, &ihi1, &
00641                 scale1[1], &abnrm1, &rcnde1[1], &rcndv1[1], &work[1], lwork, &
00642                 iwork[1], &iinfo);
00643         if (iinfo != 0) {
00644             result[1] = ulpinv;
00645             if (*jtype != 22) {
00646                 io___30.ciunit = *nounit;
00647                 s_wsfe(&io___30);
00648                 do_fio(&c__1, "SGEEVX3", (ftnlen)7);
00649                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00650                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00651                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00652                 do_fio(&c__1, balanc, (ftnlen)1);
00653                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00654                 e_wsfe();
00655             } else {
00656                 io___31.ciunit = *nounit;
00657                 s_wsfe(&io___31);
00658                 do_fio(&c__1, "SGEEVX3", (ftnlen)7);
00659                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00660                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00661                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00662                 e_wsfe();
00663             }
00664             *info = abs(iinfo);
00665             goto L190;
00666         }
00667 
00668 /*        Do Test (5) again */
00669 
00670         i__2 = *n;
00671         for (j = 1; j <= i__2; ++j) {
00672             if (wr[j] != wr1[j] || wi[j] != wi1[j]) {
00673                 result[5] = ulpinv;
00674             }
00675 /* L90: */
00676         }
00677 
00678 /*        Do Test (6) */
00679 
00680         i__2 = *n;
00681         for (j = 1; j <= i__2; ++j) {
00682             i__3 = *n;
00683             for (jj = 1; jj <= i__3; ++jj) {
00684                 if (vr[j + jj * vr_dim1] != lre[j + jj * lre_dim1]) {
00685                     result[6] = ulpinv;
00686                 }
00687 /* L100: */
00688             }
00689 /* L110: */
00690         }
00691 
00692 /*        Do Test (8) again */
00693 
00694         if (! nobal) {
00695             i__2 = *n;
00696             for (j = 1; j <= i__2; ++j) {
00697                 if (scale[j] != scale1[j]) {
00698                     result[8] = ulpinv;
00699                 }
00700 /* L120: */
00701             }
00702             if (ilo != ilo1) {
00703                 result[8] = ulpinv;
00704             }
00705             if (ihi != ihi1) {
00706                 result[8] = ulpinv;
00707             }
00708             if (abnrm != abnrm1) {
00709                 result[8] = ulpinv;
00710             }
00711         }
00712 
00713 /*        Do Test (9) again */
00714 
00715         if (isens == 2 && *n > 1) {
00716             i__2 = *n;
00717             for (j = 1; j <= i__2; ++j) {
00718                 if (rcondv[j] != rcndv1[j]) {
00719                     result[9] = ulpinv;
00720                 }
00721 /* L130: */
00722             }
00723         }
00724 
00725 /*        Compute eigenvalues and left eigenvectors, and test them */
00726 
00727         slacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00728         sgeevx_(balanc, "V", "N", sense, n, &h__[h_offset], lda, &wr1[1], &
00729                 wi1[1], &lre[lre_offset], ldlre, dum, &c__1, &ilo1, &ihi1, &
00730                 scale1[1], &abnrm1, &rcnde1[1], &rcndv1[1], &work[1], lwork, &
00731                 iwork[1], &iinfo);
00732         if (iinfo != 0) {
00733             result[1] = ulpinv;
00734             if (*jtype != 22) {
00735                 io___32.ciunit = *nounit;
00736                 s_wsfe(&io___32);
00737                 do_fio(&c__1, "SGEEVX4", (ftnlen)7);
00738                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00739                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00740                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00741                 do_fio(&c__1, balanc, (ftnlen)1);
00742                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00743                 e_wsfe();
00744             } else {
00745                 io___33.ciunit = *nounit;
00746                 s_wsfe(&io___33);
00747                 do_fio(&c__1, "SGEEVX4", (ftnlen)7);
00748                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00749                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00750                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00751                 e_wsfe();
00752             }
00753             *info = abs(iinfo);
00754             goto L190;
00755         }
00756 
00757 /*        Do Test (5) again */
00758 
00759         i__2 = *n;
00760         for (j = 1; j <= i__2; ++j) {
00761             if (wr[j] != wr1[j] || wi[j] != wi1[j]) {
00762                 result[5] = ulpinv;
00763             }
00764 /* L140: */
00765         }
00766 
00767 /*        Do Test (7) */
00768 
00769         i__2 = *n;
00770         for (j = 1; j <= i__2; ++j) {
00771             i__3 = *n;
00772             for (jj = 1; jj <= i__3; ++jj) {
00773                 if (vl[j + jj * vl_dim1] != lre[j + jj * lre_dim1]) {
00774                     result[7] = ulpinv;
00775                 }
00776 /* L150: */
00777             }
00778 /* L160: */
00779         }
00780 
00781 /*        Do Test (8) again */
00782 
00783         if (! nobal) {
00784             i__2 = *n;
00785             for (j = 1; j <= i__2; ++j) {
00786                 if (scale[j] != scale1[j]) {
00787                     result[8] = ulpinv;
00788                 }
00789 /* L170: */
00790             }
00791             if (ilo != ilo1) {
00792                 result[8] = ulpinv;
00793             }
00794             if (ihi != ihi1) {
00795                 result[8] = ulpinv;
00796             }
00797             if (abnrm != abnrm1) {
00798                 result[8] = ulpinv;
00799             }
00800         }
00801 
00802 /*        Do Test (9) again */
00803 
00804         if (isens == 2 && *n > 1) {
00805             i__2 = *n;
00806             for (j = 1; j <= i__2; ++j) {
00807                 if (rcondv[j] != rcndv1[j]) {
00808                     result[9] = ulpinv;
00809                 }
00810 /* L180: */
00811             }
00812         }
00813 
00814 L190:
00815 
00816 /* L200: */
00817         ;
00818     }
00819 
00820 /*     If COMP, compare condition numbers to precomputed ones */
00821 
00822     if (*comp) {
00823         slacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00824         sgeevx_("N", "V", "V", "B", n, &h__[h_offset], lda, &wr[1], &wi[1], &
00825                 vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &ilo, &ihi, &scale[
00826                 1], &abnrm, &rconde[1], &rcondv[1], &work[1], lwork, &iwork[1]
00827 , &iinfo);
00828         if (iinfo != 0) {
00829             result[1] = ulpinv;
00830             io___34.ciunit = *nounit;
00831             s_wsfe(&io___34);
00832             do_fio(&c__1, "SGEEVX5", (ftnlen)7);
00833             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00834             do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00835             do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00836             e_wsfe();
00837             *info = abs(iinfo);
00838             goto L250;
00839         }
00840 
00841 /*        Sort eigenvalues and condition numbers lexicographically */
00842 /*        to compare with inputs */
00843 
00844         i__1 = *n - 1;
00845         for (i__ = 1; i__ <= i__1; ++i__) {
00846             kmin = i__;
00847             vrmin = wr[i__];
00848             vimin = wi[i__];
00849             i__2 = *n;
00850             for (j = i__ + 1; j <= i__2; ++j) {
00851                 if (wr[j] < vrmin) {
00852                     kmin = j;
00853                     vrmin = wr[j];
00854                     vimin = wi[j];
00855                 }
00856 /* L210: */
00857             }
00858             wr[kmin] = wr[i__];
00859             wi[kmin] = wi[i__];
00860             wr[i__] = vrmin;
00861             wi[i__] = vimin;
00862             vrmin = rconde[kmin];
00863             rconde[kmin] = rconde[i__];
00864             rconde[i__] = vrmin;
00865             vrmin = rcondv[kmin];
00866             rcondv[kmin] = rcondv[i__];
00867             rcondv[i__] = vrmin;
00868 /* L220: */
00869         }
00870 
00871 /*        Compare condition numbers for eigenvectors */
00872 /*        taking their condition numbers into account */
00873 
00874         result[10] = 0.f;
00875         eps = dmax(5.9605e-8f,ulp);
00876 /* Computing MAX */
00877         r__1 = (real) (*n) * eps * abnrm;
00878         v = dmax(r__1,smlnum);
00879         if (abnrm == 0.f) {
00880             v = 1.f;
00881         }
00882         i__1 = *n;
00883         for (i__ = 1; i__ <= i__1; ++i__) {
00884             if (v > rcondv[i__] * rconde[i__]) {
00885                 tol = rcondv[i__];
00886             } else {
00887                 tol = v / rconde[i__];
00888             }
00889             if (v > rcdvin[i__] * rcdein[i__]) {
00890                 tolin = rcdvin[i__];
00891             } else {
00892                 tolin = v / rcdein[i__];
00893             }
00894 /* Computing MAX */
00895             r__1 = tol, r__2 = smlnum / eps;
00896             tol = dmax(r__1,r__2);
00897 /* Computing MAX */
00898             r__1 = tolin, r__2 = smlnum / eps;
00899             tolin = dmax(r__1,r__2);
00900             if (eps * (rcdvin[i__] - tolin) > rcondv[i__] + tol) {
00901                 vmax = 1.f / eps;
00902             } else if (rcdvin[i__] - tolin > rcondv[i__] + tol) {
00903                 vmax = (rcdvin[i__] - tolin) / (rcondv[i__] + tol);
00904             } else if (rcdvin[i__] + tolin < eps * (rcondv[i__] - tol)) {
00905                 vmax = 1.f / eps;
00906             } else if (rcdvin[i__] + tolin < rcondv[i__] - tol) {
00907                 vmax = (rcondv[i__] - tol) / (rcdvin[i__] + tolin);
00908             } else {
00909                 vmax = 1.f;
00910             }
00911             result[10] = dmax(result[10],vmax);
00912 /* L230: */
00913         }
00914 
00915 /*        Compare condition numbers for eigenvalues */
00916 /*        taking their condition numbers into account */
00917 
00918         result[11] = 0.f;
00919         i__1 = *n;
00920         for (i__ = 1; i__ <= i__1; ++i__) {
00921             if (v > rcondv[i__]) {
00922                 tol = 1.f;
00923             } else {
00924                 tol = v / rcondv[i__];
00925             }
00926             if (v > rcdvin[i__]) {
00927                 tolin = 1.f;
00928             } else {
00929                 tolin = v / rcdvin[i__];
00930             }
00931 /* Computing MAX */
00932             r__1 = tol, r__2 = smlnum / eps;
00933             tol = dmax(r__1,r__2);
00934 /* Computing MAX */
00935             r__1 = tolin, r__2 = smlnum / eps;
00936             tolin = dmax(r__1,r__2);
00937             if (eps * (rcdein[i__] - tolin) > rconde[i__] + tol) {
00938                 vmax = 1.f / eps;
00939             } else if (rcdein[i__] - tolin > rconde[i__] + tol) {
00940                 vmax = (rcdein[i__] - tolin) / (rconde[i__] + tol);
00941             } else if (rcdein[i__] + tolin < eps * (rconde[i__] - tol)) {
00942                 vmax = 1.f / eps;
00943             } else if (rcdein[i__] + tolin < rconde[i__] - tol) {
00944                 vmax = (rconde[i__] - tol) / (rcdein[i__] + tolin);
00945             } else {
00946                 vmax = 1.f;
00947             }
00948             result[11] = dmax(result[11],vmax);
00949 /* L240: */
00950         }
00951 L250:
00952 
00953         ;
00954     }
00955 
00956 
00957     return 0;
00958 
00959 /*     End of SGET23 */
00960 
00961 } /* sget23_ */


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