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


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