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


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