zget24.c
Go to the documentation of this file.
00001 /* zget24.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 /* Common Block Declarations */
00017 
00018 struct {
00019     integer selopt, seldim;
00020     logical selval[20];
00021     doublereal selwr[20], selwi[20];
00022 } sslct_;
00023 
00024 #define sslct_1 sslct_
00025 
00026 /* Table of constant values */
00027 
00028 static doublecomplex c_b1 = {0.,0.};
00029 static doublecomplex c_b2 = {1.,0.};
00030 static integer c__1 = 1;
00031 static integer c__4 = 4;
00032 
00033 /* Subroutine */ int zget24_(logical *comp, integer *jtype, doublereal *
00034         thresh, integer *iseed, integer *nounit, integer *n, doublecomplex *a, 
00035          integer *lda, doublecomplex *h__, doublecomplex *ht, doublecomplex *
00036         w, doublecomplex *wt, doublecomplex *wtmp, doublecomplex *vs, integer 
00037         *ldvs, doublecomplex *vs1, doublereal *rcdein, doublereal *rcdvin, 
00038         integer *nslct, integer *islct, integer *isrt, doublereal *result, 
00039         doublecomplex *work, integer *lwork, doublereal *rwork, logical *
00040         bwork, integer *info)
00041 {
00042     /* Format strings */
00043     static char fmt_9998[] = "(\002 ZGET24: \002,a,\002 returned INFO=\002,i"
00044             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, ISEED="
00045             "(\002,3(i5,\002,\002),i5,\002)\002)";
00046     static char fmt_9999[] = "(\002 ZGET24: \002,a,\002 returned INFO=\002,i"
00047             "6,\002.\002,/9x,\002N=\002,i6,\002, INPUT EXAMPLE NUMBER = \002,"
00048             "i4)";
00049 
00050     /* System generated locals */
00051     integer a_dim1, a_offset, h_dim1, h_offset, ht_dim1, ht_offset, vs_dim1, 
00052             vs_offset, vs1_dim1, vs1_offset, i__1, i__2, i__3, i__4;
00053     doublereal d__1, d__2;
00054     doublecomplex z__1;
00055 
00056     /* Builtin functions */
00057     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00058     double d_imag(doublecomplex *);
00059 
00060     /* Local variables */
00061     integer i__, j;
00062     doublereal v, eps, tol, ulp;
00063     integer sdim, kmin;
00064     doublecomplex ctmp;
00065     integer itmp, ipnt[20], rsub;
00066     char sort[1];
00067     integer sdim1, iinfo;
00068     doublereal anorm;
00069     extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 
00070             integer *, doublecomplex *, doublecomplex *, integer *, 
00071             doublecomplex *, integer *, doublecomplex *, doublecomplex *, 
00072             integer *);
00073     doublereal tolin;
00074     integer isort;
00075     extern /* Subroutine */ int zunt01_(char *, integer *, integer *, 
00076             doublecomplex *, integer *, doublecomplex *, integer *, 
00077             doublereal *, doublereal *);
00078     doublereal wnorm;
00079     extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
00080             doublecomplex *, integer *);
00081     doublereal rcnde1, rcndv1;
00082     extern doublereal dlamch_(char *);
00083     doublereal rconde;
00084     extern /* Subroutine */ int xerbla_(char *, integer *);
00085     extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 
00086             integer *, doublereal *);
00087     integer knteig;
00088     doublereal rcondv, vricmp;
00089     extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 
00090             doublecomplex *, integer *, doublecomplex *, integer *);
00091     doublereal vrimin;
00092     extern logical zslect_(doublecomplex *);
00093     extern /* Subroutine */ int zgeesx_(char *, char *, L_fp, char *, integer 
00094             *, doublecomplex *, integer *, integer *, doublecomplex *, 
00095             doublecomplex *, integer *, doublereal *, doublereal *, 
00096             doublecomplex *, integer *, doublereal *, logical *, integer *);
00097     doublereal smlnum, ulpinv;
00098 
00099     /* Fortran I/O blocks */
00100     static cilist io___12 = { 0, 0, 0, fmt_9998, 0 };
00101     static cilist io___13 = { 0, 0, 0, fmt_9999, 0 };
00102     static cilist io___17 = { 0, 0, 0, fmt_9998, 0 };
00103     static cilist io___18 = { 0, 0, 0, fmt_9999, 0 };
00104     static cilist io___21 = { 0, 0, 0, fmt_9998, 0 };
00105     static cilist io___22 = { 0, 0, 0, fmt_9999, 0 };
00106     static cilist io___25 = { 0, 0, 0, fmt_9998, 0 };
00107     static cilist io___26 = { 0, 0, 0, fmt_9999, 0 };
00108     static cilist io___27 = { 0, 0, 0, fmt_9998, 0 };
00109     static cilist io___28 = { 0, 0, 0, fmt_9999, 0 };
00110     static cilist io___29 = { 0, 0, 0, fmt_9998, 0 };
00111     static cilist io___30 = { 0, 0, 0, fmt_9999, 0 };
00112     static cilist io___31 = { 0, 0, 0, fmt_9998, 0 };
00113     static cilist io___32 = { 0, 0, 0, fmt_9999, 0 };
00114     static cilist io___33 = { 0, 0, 0, fmt_9998, 0 };
00115     static cilist io___34 = { 0, 0, 0, fmt_9999, 0 };
00116     static cilist io___42 = { 0, 0, 0, fmt_9999, 0 };
00117 
00118 
00119 
00120 /*  -- LAPACK test routine (version 3.1) -- */
00121 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00122 /*     November 2006 */
00123 
00124 /*     .. Scalar Arguments .. */
00125 /*     .. */
00126 /*     .. Array Arguments .. */
00127 /*     .. */
00128 
00129 /*  Purpose */
00130 /*  ======= */
00131 
00132 /*     ZGET24 checks the nonsymmetric eigenvalue (Schur form) problem */
00133 /*     expert driver ZGEESX. */
00134 
00135 /*     If COMP = .FALSE., the first 13 of the following tests will be */
00136 /*     be performed on the input matrix A, and also tests 14 and 15 */
00137 /*     if LWORK is sufficiently large. */
00138 /*     If COMP = .TRUE., all 17 test will be performed. */
00139 
00140 /*     (1)     0 if T is in Schur form, 1/ulp otherwise */
00141 /*            (no sorting of eigenvalues) */
00142 
00143 /*     (2)     | A - VS T VS' | / ( n |A| ulp ) */
00144 
00145 /*       Here VS is the matrix of Schur eigenvectors, and T is in Schur */
00146 /*       form  (no sorting of eigenvalues). */
00147 
00148 /*     (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues). */
00149 
00150 /*     (4)     0     if W are eigenvalues of T */
00151 /*             1/ulp otherwise */
00152 /*             (no sorting of eigenvalues) */
00153 
00154 /*     (5)     0     if T(with VS) = T(without VS), */
00155 /*             1/ulp otherwise */
00156 /*             (no sorting of eigenvalues) */
00157 
00158 /*     (6)     0     if eigenvalues(with VS) = eigenvalues(without VS), */
00159 /*             1/ulp otherwise */
00160 /*             (no sorting of eigenvalues) */
00161 
00162 /*     (7)     0 if T is in Schur form, 1/ulp otherwise */
00163 /*             (with sorting of eigenvalues) */
00164 
00165 /*     (8)     | A - VS T VS' | / ( n |A| ulp ) */
00166 
00167 /*       Here VS is the matrix of Schur eigenvectors, and T is in Schur */
00168 /*       form  (with sorting of eigenvalues). */
00169 
00170 /*     (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues). */
00171 
00172 /*     (10)    0     if W are eigenvalues of T */
00173 /*             1/ulp otherwise */
00174 /*             If workspace sufficient, also compare W with and */
00175 /*             without reciprocal condition numbers */
00176 /*             (with sorting of eigenvalues) */
00177 
00178 /*     (11)    0     if T(with VS) = T(without VS), */
00179 /*             1/ulp otherwise */
00180 /*             If workspace sufficient, also compare T with and without */
00181 /*             reciprocal condition numbers */
00182 /*             (with sorting of eigenvalues) */
00183 
00184 /*     (12)    0     if eigenvalues(with VS) = eigenvalues(without VS), */
00185 /*             1/ulp otherwise */
00186 /*             If workspace sufficient, also compare VS with and without */
00187 /*             reciprocal condition numbers */
00188 /*             (with sorting of eigenvalues) */
00189 
00190 /*     (13)    if sorting worked and SDIM is the number of */
00191 /*             eigenvalues which were SELECTed */
00192 /*             If workspace sufficient, also compare SDIM with and */
00193 /*             without reciprocal condition numbers */
00194 
00195 /*     (14)    if RCONDE the same no matter if VS and/or RCONDV computed */
00196 
00197 /*     (15)    if RCONDV the same no matter if VS and/or RCONDE computed */
00198 
00199 /*     (16)  |RCONDE - RCDEIN| / cond(RCONDE) */
00200 
00201 /*        RCONDE is the reciprocal average eigenvalue condition number */
00202 /*        computed by ZGEESX and RCDEIN (the precomputed true value) */
00203 /*        is supplied as input.  cond(RCONDE) is the condition number */
00204 /*        of RCONDE, and takes errors in computing RCONDE into account, */
00205 /*        so that the resulting quantity should be O(ULP). cond(RCONDE) */
00206 /*        is essentially given by norm(A)/RCONDV. */
00207 
00208 /*     (17)  |RCONDV - RCDVIN| / cond(RCONDV) */
00209 
00210 /*        RCONDV is the reciprocal right invariant subspace condition */
00211 /*        number computed by ZGEESX and RCDVIN (the precomputed true */
00212 /*        value) is supplied as input. cond(RCONDV) is the condition */
00213 /*        number of RCONDV, and takes errors in computing RCONDV into */
00214 /*        account, so that the resulting quantity should be O(ULP). */
00215 /*        cond(RCONDV) is essentially given by norm(A)/RCONDE. */
00216 
00217 /*  Arguments */
00218 /*  ========= */
00219 
00220 /*  COMP    (input) LOGICAL */
00221 /*          COMP describes which input tests to perform: */
00222 /*            = .FALSE. if the computed condition numbers are not to */
00223 /*                      be tested against RCDVIN and RCDEIN */
00224 /*            = .TRUE.  if they are to be compared */
00225 
00226 /*  JTYPE   (input) INTEGER */
00227 /*          Type of input matrix. Used to label output if error occurs. */
00228 
00229 /*  ISEED   (input) INTEGER array, dimension (4) */
00230 /*          If COMP = .FALSE., the random number generator seed */
00231 /*          used to produce matrix. */
00232 /*          If COMP = .TRUE., ISEED(1) = the number of the example. */
00233 /*          Used to label output if error occurs. */
00234 
00235 /*  THRESH  (input) DOUBLE PRECISION */
00236 /*          A test will count as "failed" if the "error", computed as */
00237 /*          described above, exceeds THRESH.  Note that the error */
00238 /*          is scaled to be O(1), so THRESH should be a reasonably */
00239 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00240 /*          it should not depend on the precision (single vs. double) */
00241 /*          or the size of the matrix.  It must be at least zero. */
00242 
00243 /*  NOUNIT  (input) INTEGER */
00244 /*          The FORTRAN unit number for printing out error messages */
00245 /*          (e.g., if a routine returns INFO not equal to 0.) */
00246 
00247 /*  N       (input) INTEGER */
00248 /*          The dimension of A. N must be at least 0. */
00249 
00250 /*  A       (input/output) COMPLEX*16 array, dimension (LDA, N) */
00251 /*          Used to hold the matrix whose eigenvalues are to be */
00252 /*          computed. */
00253 
00254 /*  LDA     (input) INTEGER */
00255 /*          The leading dimension of A, and H. LDA must be at */
00256 /*          least 1 and at least N. */
00257 
00258 /*  H       (workspace) COMPLEX*16 array, dimension (LDA, N) */
00259 /*          Another copy of the test matrix A, modified by ZGEESX. */
00260 
00261 /*  HT      (workspace) COMPLEX*16 array, dimension (LDA, N) */
00262 /*          Yet another copy of the test matrix A, modified by ZGEESX. */
00263 
00264 /*  W       (workspace) COMPLEX*16 array, dimension (N) */
00265 /*          The computed eigenvalues of A. */
00266 
00267 /*  WT      (workspace) COMPLEX*16 array, dimension (N) */
00268 /*          Like W, this array contains the eigenvalues of A, */
00269 /*          but those computed when ZGEESX only computes a partial */
00270 /*          eigendecomposition, i.e. not Schur vectors */
00271 
00272 /*  WTMP    (workspace) COMPLEX*16 array, dimension (N) */
00273 /*          Like W, this array contains the eigenvalues of A, */
00274 /*          but sorted by increasing real or imaginary part. */
00275 
00276 /*  VS      (workspace) COMPLEX*16 array, dimension (LDVS, N) */
00277 /*          VS holds the computed Schur vectors. */
00278 
00279 /*  LDVS    (input) INTEGER */
00280 /*          Leading dimension of VS. Must be at least max(1, N). */
00281 
00282 /*  VS1     (workspace) COMPLEX*16 array, dimension (LDVS, N) */
00283 /*          VS1 holds another copy of the computed Schur vectors. */
00284 
00285 /*  RCDEIN  (input) DOUBLE PRECISION */
00286 /*          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal */
00287 /*          condition number for the average of selected eigenvalues. */
00288 
00289 /*  RCDVIN  (input) DOUBLE PRECISION */
00290 /*          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal */
00291 /*          condition number for the selected right invariant subspace. */
00292 
00293 /*  NSLCT   (input) INTEGER */
00294 /*          When COMP = .TRUE. the number of selected eigenvalues */
00295 /*          corresponding to the precomputed values RCDEIN and RCDVIN. */
00296 
00297 /*  ISLCT   (input) INTEGER array, dimension (NSLCT) */
00298 /*          When COMP = .TRUE. ISLCT selects the eigenvalues of the */
00299 /*          input matrix corresponding to the precomputed values RCDEIN */
00300 /*          and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the */
00301 /*          eigenvalue with the J-th largest real or imaginary part is */
00302 /*          selected. The real part is used if ISRT = 0, and the */
00303 /*          imaginary part if ISRT = 1. */
00304 /*          Not referenced if COMP = .FALSE. */
00305 
00306 /*  ISRT    (input) INTEGER */
00307 /*          When COMP = .TRUE., ISRT describes how ISLCT is used to */
00308 /*          choose a subset of the spectrum. */
00309 /*          Not referenced if COMP = .FALSE. */
00310 
00311 /*  RESULT  (output) DOUBLE PRECISION array, dimension (17) */
00312 /*          The values computed by the 17 tests described above. */
00313 /*          The values are currently limited to 1/ulp, to avoid */
00314 /*          overflow. */
00315 
00316 /*  WORK    (workspace) COMPLEX*16 array, dimension (2*N*N) */
00317 
00318 /*  LWORK   (input) INTEGER */
00319 /*          The number of entries in WORK to be passed to ZGEESX. This */
00320 /*          must be at least 2*N, and N*(N+1)/2 if tests 14--16 are to */
00321 /*          be performed. */
00322 
00323 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N) */
00324 
00325 /*  BWORK   (workspace) LOGICAL array, dimension (N) */
00326 
00327 /*  INFO    (output) INTEGER */
00328 /*          If 0,  successful exit. */
00329 /*          If <0, input parameter -INFO had an incorrect value. */
00330 /*          If >0, ZGEESX returned an error code, the absolute */
00331 /*                 value of which is returned. */
00332 
00333 /*  ===================================================================== */
00334 
00335 /*     .. Parameters .. */
00336 /*     .. */
00337 /*     .. Local Scalars .. */
00338 /*     .. */
00339 /*     .. Local Arrays .. */
00340 /*     .. */
00341 /*     .. External Functions .. */
00342 /*     .. */
00343 /*     .. External Subroutines .. */
00344 /*     .. */
00345 /*     .. Intrinsic Functions .. */
00346 /*     .. */
00347 /*     .. Arrays in Common .. */
00348 /*     .. */
00349 /*     .. Scalars in Common .. */
00350 /*     .. */
00351 /*     .. Common blocks .. */
00352 /*     .. */
00353 /*     .. Executable Statements .. */
00354 
00355 /*     Check for errors */
00356 
00357     /* Parameter adjustments */
00358     --iseed;
00359     ht_dim1 = *lda;
00360     ht_offset = 1 + ht_dim1;
00361     ht -= ht_offset;
00362     h_dim1 = *lda;
00363     h_offset = 1 + h_dim1;
00364     h__ -= h_offset;
00365     a_dim1 = *lda;
00366     a_offset = 1 + a_dim1;
00367     a -= a_offset;
00368     --w;
00369     --wt;
00370     --wtmp;
00371     vs1_dim1 = *ldvs;
00372     vs1_offset = 1 + vs1_dim1;
00373     vs1 -= vs1_offset;
00374     vs_dim1 = *ldvs;
00375     vs_offset = 1 + vs_dim1;
00376     vs -= vs_offset;
00377     --islct;
00378     --result;
00379     --work;
00380     --rwork;
00381     --bwork;
00382 
00383     /* Function Body */
00384     *info = 0;
00385     if (*thresh < 0.) {
00386         *info = -3;
00387     } else if (*nounit <= 0) {
00388         *info = -5;
00389     } else if (*n < 0) {
00390         *info = -6;
00391     } else if (*lda < 1 || *lda < *n) {
00392         *info = -8;
00393     } else if (*ldvs < 1 || *ldvs < *n) {
00394         *info = -15;
00395     } else if (*lwork < *n << 1) {
00396         *info = -24;
00397     }
00398 
00399     if (*info != 0) {
00400         i__1 = -(*info);
00401         xerbla_("ZGET24", &i__1);
00402         return 0;
00403     }
00404 
00405 /*     Quick return if nothing to do */
00406 
00407     for (i__ = 1; i__ <= 17; ++i__) {
00408         result[i__] = -1.;
00409 /* L10: */
00410     }
00411 
00412     if (*n == 0) {
00413         return 0;
00414     }
00415 
00416 /*     Important constants */
00417 
00418     smlnum = dlamch_("Safe minimum");
00419     ulp = dlamch_("Precision");
00420     ulpinv = 1. / ulp;
00421 
00422 /*     Perform tests (1)-(13) */
00423 
00424     sslct_1.selopt = 0;
00425     for (isort = 0; isort <= 1; ++isort) {
00426         if (isort == 0) {
00427             *(unsigned char *)sort = 'N';
00428             rsub = 0;
00429         } else {
00430             *(unsigned char *)sort = 'S';
00431             rsub = 6;
00432         }
00433 
00434 /*        Compute Schur form and Schur vectors, and test them */
00435 
00436         zlacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00437         zgeesx_("V", sort, (L_fp)zslect_, "N", n, &h__[h_offset], lda, &sdim, 
00438                 &w[1], &vs[vs_offset], ldvs, &rconde, &rcondv, &work[1], 
00439                 lwork, &rwork[1], &bwork[1], &iinfo);
00440         if (iinfo != 0) {
00441             result[rsub + 1] = ulpinv;
00442             if (*jtype != 22) {
00443                 io___12.ciunit = *nounit;
00444                 s_wsfe(&io___12);
00445                 do_fio(&c__1, "ZGEESX1", (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__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00450                 e_wsfe();
00451             } else {
00452                 io___13.ciunit = *nounit;
00453                 s_wsfe(&io___13);
00454                 do_fio(&c__1, "ZGEESX1", (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         if (isort == 0) {
00464             zcopy_(n, &w[1], &c__1, &wtmp[1], &c__1);
00465         }
00466 
00467 /*        Do Test (1) or Test (7) */
00468 
00469         result[rsub + 1] = 0.;
00470         i__1 = *n - 1;
00471         for (j = 1; j <= i__1; ++j) {
00472             i__2 = *n;
00473             for (i__ = j + 1; i__ <= i__2; ++i__) {
00474                 i__3 = i__ + j * h_dim1;
00475                 if (h__[i__3].r != 0. || h__[i__3].i != 0.) {
00476                     result[rsub + 1] = ulpinv;
00477                 }
00478 /* L20: */
00479             }
00480 /* L30: */
00481         }
00482 
00483 /*        Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP) */
00484 
00485 /*        Copy A to VS1, used as workspace */
00486 
00487         zlacpy_(" ", n, n, &a[a_offset], lda, &vs1[vs1_offset], ldvs);
00488 
00489 /*        Compute Q*H and store in HT. */
00490 
00491         zgemm_("No transpose", "No transpose", n, n, n, &c_b2, &vs[vs_offset], 
00492                  ldvs, &h__[h_offset], lda, &c_b1, &ht[ht_offset], lda);
00493 
00494 /*        Compute A - Q*H*Q' */
00495 
00496         z__1.r = -1., z__1.i = -0.;
00497         zgemm_("No transpose", "Conjugate transpose", n, n, n, &z__1, &ht[
00498                 ht_offset], lda, &vs[vs_offset], ldvs, &c_b2, &vs1[vs1_offset]
00499 , ldvs);
00500 
00501 /* Computing MAX */
00502         d__1 = zlange_("1", n, n, &a[a_offset], lda, &rwork[1]);
00503         anorm = max(d__1,smlnum);
00504         wnorm = zlange_("1", n, n, &vs1[vs1_offset], ldvs, &rwork[1]);
00505 
00506         if (anorm > wnorm) {
00507             result[rsub + 2] = wnorm / anorm / (*n * ulp);
00508         } else {
00509             if (anorm < 1.) {
00510 /* Computing MIN */
00511                 d__1 = wnorm, d__2 = *n * anorm;
00512                 result[rsub + 2] = min(d__1,d__2) / anorm / (*n * ulp);
00513             } else {
00514 /* Computing MIN */
00515                 d__1 = wnorm / anorm, d__2 = (doublereal) (*n);
00516                 result[rsub + 2] = min(d__1,d__2) / (*n * ulp);
00517             }
00518         }
00519 
00520 /*        Test (3) or (9):  Compute norm( I - Q'*Q ) / ( N * ULP ) */
00521 
00522         zunt01_("Columns", n, n, &vs[vs_offset], ldvs, &work[1], lwork, &
00523                 rwork[1], &result[rsub + 3]);
00524 
00525 /*        Do Test (4) or Test (10) */
00526 
00527         result[rsub + 4] = 0.;
00528         i__1 = *n;
00529         for (i__ = 1; i__ <= i__1; ++i__) {
00530             i__2 = i__ + i__ * h_dim1;
00531             i__3 = i__;
00532             if (h__[i__2].r != w[i__3].r || h__[i__2].i != w[i__3].i) {
00533                 result[rsub + 4] = ulpinv;
00534             }
00535 /* L40: */
00536         }
00537 
00538 /*        Do Test (5) or Test (11) */
00539 
00540         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
00541         zgeesx_("N", sort, (L_fp)zslect_, "N", n, &ht[ht_offset], lda, &sdim, 
00542                 &wt[1], &vs[vs_offset], ldvs, &rconde, &rcondv, &work[1], 
00543                 lwork, &rwork[1], &bwork[1], &iinfo);
00544         if (iinfo != 0) {
00545             result[rsub + 5] = ulpinv;
00546             if (*jtype != 22) {
00547                 io___17.ciunit = *nounit;
00548                 s_wsfe(&io___17);
00549                 do_fio(&c__1, "ZGEESX2", (ftnlen)7);
00550                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00551                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00552                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00553                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00554                 e_wsfe();
00555             } else {
00556                 io___18.ciunit = *nounit;
00557                 s_wsfe(&io___18);
00558                 do_fio(&c__1, "ZGEESX2", (ftnlen)7);
00559                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00560                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00561                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00562                 e_wsfe();
00563             }
00564             *info = abs(iinfo);
00565             goto L220;
00566         }
00567 
00568         result[rsub + 5] = 0.;
00569         i__1 = *n;
00570         for (j = 1; j <= i__1; ++j) {
00571             i__2 = *n;
00572             for (i__ = 1; i__ <= i__2; ++i__) {
00573                 i__3 = i__ + j * h_dim1;
00574                 i__4 = i__ + j * ht_dim1;
00575                 if (h__[i__3].r != ht[i__4].r || h__[i__3].i != ht[i__4].i) {
00576                     result[rsub + 5] = ulpinv;
00577                 }
00578 /* L50: */
00579             }
00580 /* L60: */
00581         }
00582 
00583 /*        Do Test (6) or Test (12) */
00584 
00585         result[rsub + 6] = 0.;
00586         i__1 = *n;
00587         for (i__ = 1; i__ <= i__1; ++i__) {
00588             i__2 = i__;
00589             i__3 = i__;
00590             if (w[i__2].r != wt[i__3].r || w[i__2].i != wt[i__3].i) {
00591                 result[rsub + 6] = ulpinv;
00592             }
00593 /* L70: */
00594         }
00595 
00596 /*        Do Test (13) */
00597 
00598         if (isort == 1) {
00599             result[13] = 0.;
00600             knteig = 0;
00601             i__1 = *n;
00602             for (i__ = 1; i__ <= i__1; ++i__) {
00603                 if (zslect_(&w[i__])) {
00604                     ++knteig;
00605                 }
00606                 if (i__ < *n) {
00607                     if (zslect_(&w[i__ + 1]) && ! zslect_(&w[i__])) {
00608                         result[13] = ulpinv;
00609                     }
00610                 }
00611 /* L80: */
00612             }
00613             if (sdim != knteig) {
00614                 result[13] = ulpinv;
00615             }
00616         }
00617 
00618 /* L90: */
00619     }
00620 
00621 /*     If there is enough workspace, perform tests (14) and (15) */
00622 /*     as well as (10) through (13) */
00623 
00624     if (*lwork >= *n * (*n + 1) / 2) {
00625 
00626 /*        Compute both RCONDE and RCONDV with VS */
00627 
00628         *(unsigned char *)sort = 'S';
00629         result[14] = 0.;
00630         result[15] = 0.;
00631         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
00632         zgeesx_("V", sort, (L_fp)zslect_, "B", n, &ht[ht_offset], lda, &sdim1, 
00633                  &wt[1], &vs1[vs1_offset], ldvs, &rconde, &rcondv, &work[1], 
00634                 lwork, &rwork[1], &bwork[1], &iinfo);
00635         if (iinfo != 0) {
00636             result[14] = ulpinv;
00637             result[15] = ulpinv;
00638             if (*jtype != 22) {
00639                 io___21.ciunit = *nounit;
00640                 s_wsfe(&io___21);
00641                 do_fio(&c__1, "ZGEESX3", (ftnlen)7);
00642                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00643                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00644                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00645                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00646                 e_wsfe();
00647             } else {
00648                 io___22.ciunit = *nounit;
00649                 s_wsfe(&io___22);
00650                 do_fio(&c__1, "ZGEESX3", (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 *)&iseed[1], (ftnlen)sizeof(integer));
00654                 e_wsfe();
00655             }
00656             *info = abs(iinfo);
00657             goto L220;
00658         }
00659 
00660 /*        Perform tests (10), (11), (12), and (13) */
00661 
00662         i__1 = *n;
00663         for (i__ = 1; i__ <= i__1; ++i__) {
00664             i__2 = i__;
00665             i__3 = i__;
00666             if (w[i__2].r != wt[i__3].r || w[i__2].i != wt[i__3].i) {
00667                 result[10] = ulpinv;
00668             }
00669             i__2 = *n;
00670             for (j = 1; j <= i__2; ++j) {
00671                 i__3 = i__ + j * h_dim1;
00672                 i__4 = i__ + j * ht_dim1;
00673                 if (h__[i__3].r != ht[i__4].r || h__[i__3].i != ht[i__4].i) {
00674                     result[11] = ulpinv;
00675                 }
00676                 i__3 = i__ + j * vs_dim1;
00677                 i__4 = i__ + j * vs1_dim1;
00678                 if (vs[i__3].r != vs1[i__4].r || vs[i__3].i != vs1[i__4].i) {
00679                     result[12] = ulpinv;
00680                 }
00681 /* L100: */
00682             }
00683 /* L110: */
00684         }
00685         if (sdim != sdim1) {
00686             result[13] = ulpinv;
00687         }
00688 
00689 /*        Compute both RCONDE and RCONDV without VS, and compare */
00690 
00691         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
00692         zgeesx_("N", sort, (L_fp)zslect_, "B", n, &ht[ht_offset], lda, &sdim1, 
00693                  &wt[1], &vs1[vs1_offset], ldvs, &rcnde1, &rcndv1, &work[1], 
00694                 lwork, &rwork[1], &bwork[1], &iinfo);
00695         if (iinfo != 0) {
00696             result[14] = ulpinv;
00697             result[15] = ulpinv;
00698             if (*jtype != 22) {
00699                 io___25.ciunit = *nounit;
00700                 s_wsfe(&io___25);
00701                 do_fio(&c__1, "ZGEESX4", (ftnlen)7);
00702                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00703                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00704                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00705                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00706                 e_wsfe();
00707             } else {
00708                 io___26.ciunit = *nounit;
00709                 s_wsfe(&io___26);
00710                 do_fio(&c__1, "ZGEESX4", (ftnlen)7);
00711                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00712                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00713                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00714                 e_wsfe();
00715             }
00716             *info = abs(iinfo);
00717             goto L220;
00718         }
00719 
00720 /*        Perform tests (14) and (15) */
00721 
00722         if (rcnde1 != rconde) {
00723             result[14] = ulpinv;
00724         }
00725         if (rcndv1 != rcondv) {
00726             result[15] = ulpinv;
00727         }
00728 
00729 /*        Perform tests (10), (11), (12), and (13) */
00730 
00731         i__1 = *n;
00732         for (i__ = 1; i__ <= i__1; ++i__) {
00733             i__2 = i__;
00734             i__3 = i__;
00735             if (w[i__2].r != wt[i__3].r || w[i__2].i != wt[i__3].i) {
00736                 result[10] = ulpinv;
00737             }
00738             i__2 = *n;
00739             for (j = 1; j <= i__2; ++j) {
00740                 i__3 = i__ + j * h_dim1;
00741                 i__4 = i__ + j * ht_dim1;
00742                 if (h__[i__3].r != ht[i__4].r || h__[i__3].i != ht[i__4].i) {
00743                     result[11] = ulpinv;
00744                 }
00745                 i__3 = i__ + j * vs_dim1;
00746                 i__4 = i__ + j * vs1_dim1;
00747                 if (vs[i__3].r != vs1[i__4].r || vs[i__3].i != vs1[i__4].i) {
00748                     result[12] = ulpinv;
00749                 }
00750 /* L120: */
00751             }
00752 /* L130: */
00753         }
00754         if (sdim != sdim1) {
00755             result[13] = ulpinv;
00756         }
00757 
00758 /*        Compute RCONDE with VS, and compare */
00759 
00760         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
00761         zgeesx_("V", sort, (L_fp)zslect_, "E", n, &ht[ht_offset], lda, &sdim1, 
00762                  &wt[1], &vs1[vs1_offset], ldvs, &rcnde1, &rcndv1, &work[1], 
00763                 lwork, &rwork[1], &bwork[1], &iinfo);
00764         if (iinfo != 0) {
00765             result[14] = ulpinv;
00766             if (*jtype != 22) {
00767                 io___27.ciunit = *nounit;
00768                 s_wsfe(&io___27);
00769                 do_fio(&c__1, "ZGEESX5", (ftnlen)7);
00770                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00771                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00772                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00773                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00774                 e_wsfe();
00775             } else {
00776                 io___28.ciunit = *nounit;
00777                 s_wsfe(&io___28);
00778                 do_fio(&c__1, "ZGEESX5", (ftnlen)7);
00779                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00780                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00781                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00782                 e_wsfe();
00783             }
00784             *info = abs(iinfo);
00785             goto L220;
00786         }
00787 
00788 /*        Perform test (14) */
00789 
00790         if (rcnde1 != rconde) {
00791             result[14] = ulpinv;
00792         }
00793 
00794 /*        Perform tests (10), (11), (12), and (13) */
00795 
00796         i__1 = *n;
00797         for (i__ = 1; i__ <= i__1; ++i__) {
00798             i__2 = i__;
00799             i__3 = i__;
00800             if (w[i__2].r != wt[i__3].r || w[i__2].i != wt[i__3].i) {
00801                 result[10] = ulpinv;
00802             }
00803             i__2 = *n;
00804             for (j = 1; j <= i__2; ++j) {
00805                 i__3 = i__ + j * h_dim1;
00806                 i__4 = i__ + j * ht_dim1;
00807                 if (h__[i__3].r != ht[i__4].r || h__[i__3].i != ht[i__4].i) {
00808                     result[11] = ulpinv;
00809                 }
00810                 i__3 = i__ + j * vs_dim1;
00811                 i__4 = i__ + j * vs1_dim1;
00812                 if (vs[i__3].r != vs1[i__4].r || vs[i__3].i != vs1[i__4].i) {
00813                     result[12] = ulpinv;
00814                 }
00815 /* L140: */
00816             }
00817 /* L150: */
00818         }
00819         if (sdim != sdim1) {
00820             result[13] = ulpinv;
00821         }
00822 
00823 /*        Compute RCONDE without VS, and compare */
00824 
00825         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
00826         zgeesx_("N", sort, (L_fp)zslect_, "E", n, &ht[ht_offset], lda, &sdim1, 
00827                  &wt[1], &vs1[vs1_offset], ldvs, &rcnde1, &rcndv1, &work[1], 
00828                 lwork, &rwork[1], &bwork[1], &iinfo);
00829         if (iinfo != 0) {
00830             result[14] = ulpinv;
00831             if (*jtype != 22) {
00832                 io___29.ciunit = *nounit;
00833                 s_wsfe(&io___29);
00834                 do_fio(&c__1, "ZGEESX6", (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 *)&(*jtype), (ftnlen)sizeof(integer));
00838                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00839                 e_wsfe();
00840             } else {
00841                 io___30.ciunit = *nounit;
00842                 s_wsfe(&io___30);
00843                 do_fio(&c__1, "ZGEESX6", (ftnlen)7);
00844                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00845                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00846                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00847                 e_wsfe();
00848             }
00849             *info = abs(iinfo);
00850             goto L220;
00851         }
00852 
00853 /*        Perform test (14) */
00854 
00855         if (rcnde1 != rconde) {
00856             result[14] = ulpinv;
00857         }
00858 
00859 /*        Perform tests (10), (11), (12), and (13) */
00860 
00861         i__1 = *n;
00862         for (i__ = 1; i__ <= i__1; ++i__) {
00863             i__2 = i__;
00864             i__3 = i__;
00865             if (w[i__2].r != wt[i__3].r || w[i__2].i != wt[i__3].i) {
00866                 result[10] = ulpinv;
00867             }
00868             i__2 = *n;
00869             for (j = 1; j <= i__2; ++j) {
00870                 i__3 = i__ + j * h_dim1;
00871                 i__4 = i__ + j * ht_dim1;
00872                 if (h__[i__3].r != ht[i__4].r || h__[i__3].i != ht[i__4].i) {
00873                     result[11] = ulpinv;
00874                 }
00875                 i__3 = i__ + j * vs_dim1;
00876                 i__4 = i__ + j * vs1_dim1;
00877                 if (vs[i__3].r != vs1[i__4].r || vs[i__3].i != vs1[i__4].i) {
00878                     result[12] = ulpinv;
00879                 }
00880 /* L160: */
00881             }
00882 /* L170: */
00883         }
00884         if (sdim != sdim1) {
00885             result[13] = ulpinv;
00886         }
00887 
00888 /*        Compute RCONDV with VS, and compare */
00889 
00890         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
00891         zgeesx_("V", sort, (L_fp)zslect_, "V", n, &ht[ht_offset], lda, &sdim1, 
00892                  &wt[1], &vs1[vs1_offset], ldvs, &rcnde1, &rcndv1, &work[1], 
00893                 lwork, &rwork[1], &bwork[1], &iinfo);
00894         if (iinfo != 0) {
00895             result[15] = ulpinv;
00896             if (*jtype != 22) {
00897                 io___31.ciunit = *nounit;
00898                 s_wsfe(&io___31);
00899                 do_fio(&c__1, "ZGEESX7", (ftnlen)7);
00900                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00901                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00902                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00903                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00904                 e_wsfe();
00905             } else {
00906                 io___32.ciunit = *nounit;
00907                 s_wsfe(&io___32);
00908                 do_fio(&c__1, "ZGEESX7", (ftnlen)7);
00909                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00910                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00911                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00912                 e_wsfe();
00913             }
00914             *info = abs(iinfo);
00915             goto L220;
00916         }
00917 
00918 /*        Perform test (15) */
00919 
00920         if (rcndv1 != rcondv) {
00921             result[15] = ulpinv;
00922         }
00923 
00924 /*        Perform tests (10), (11), (12), and (13) */
00925 
00926         i__1 = *n;
00927         for (i__ = 1; i__ <= i__1; ++i__) {
00928             i__2 = i__;
00929             i__3 = i__;
00930             if (w[i__2].r != wt[i__3].r || w[i__2].i != wt[i__3].i) {
00931                 result[10] = ulpinv;
00932             }
00933             i__2 = *n;
00934             for (j = 1; j <= i__2; ++j) {
00935                 i__3 = i__ + j * h_dim1;
00936                 i__4 = i__ + j * ht_dim1;
00937                 if (h__[i__3].r != ht[i__4].r || h__[i__3].i != ht[i__4].i) {
00938                     result[11] = ulpinv;
00939                 }
00940                 i__3 = i__ + j * vs_dim1;
00941                 i__4 = i__ + j * vs1_dim1;
00942                 if (vs[i__3].r != vs1[i__4].r || vs[i__3].i != vs1[i__4].i) {
00943                     result[12] = ulpinv;
00944                 }
00945 /* L180: */
00946             }
00947 /* L190: */
00948         }
00949         if (sdim != sdim1) {
00950             result[13] = ulpinv;
00951         }
00952 
00953 /*        Compute RCONDV without VS, and compare */
00954 
00955         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
00956         zgeesx_("N", sort, (L_fp)zslect_, "V", n, &ht[ht_offset], lda, &sdim1, 
00957                  &wt[1], &vs1[vs1_offset], ldvs, &rcnde1, &rcndv1, &work[1], 
00958                 lwork, &rwork[1], &bwork[1], &iinfo);
00959         if (iinfo != 0) {
00960             result[15] = ulpinv;
00961             if (*jtype != 22) {
00962                 io___33.ciunit = *nounit;
00963                 s_wsfe(&io___33);
00964                 do_fio(&c__1, "ZGEESX8", (ftnlen)7);
00965                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00966                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00967                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00968                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00969                 e_wsfe();
00970             } else {
00971                 io___34.ciunit = *nounit;
00972                 s_wsfe(&io___34);
00973                 do_fio(&c__1, "ZGEESX8", (ftnlen)7);
00974                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00975                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00976                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00977                 e_wsfe();
00978             }
00979             *info = abs(iinfo);
00980             goto L220;
00981         }
00982 
00983 /*        Perform test (15) */
00984 
00985         if (rcndv1 != rcondv) {
00986             result[15] = ulpinv;
00987         }
00988 
00989 /*        Perform tests (10), (11), (12), and (13) */
00990 
00991         i__1 = *n;
00992         for (i__ = 1; i__ <= i__1; ++i__) {
00993             i__2 = i__;
00994             i__3 = i__;
00995             if (w[i__2].r != wt[i__3].r || w[i__2].i != wt[i__3].i) {
00996                 result[10] = ulpinv;
00997             }
00998             i__2 = *n;
00999             for (j = 1; j <= i__2; ++j) {
01000                 i__3 = i__ + j * h_dim1;
01001                 i__4 = i__ + j * ht_dim1;
01002                 if (h__[i__3].r != ht[i__4].r || h__[i__3].i != ht[i__4].i) {
01003                     result[11] = ulpinv;
01004                 }
01005                 i__3 = i__ + j * vs_dim1;
01006                 i__4 = i__ + j * vs1_dim1;
01007                 if (vs[i__3].r != vs1[i__4].r || vs[i__3].i != vs1[i__4].i) {
01008                     result[12] = ulpinv;
01009                 }
01010 /* L200: */
01011             }
01012 /* L210: */
01013         }
01014         if (sdim != sdim1) {
01015             result[13] = ulpinv;
01016         }
01017 
01018     }
01019 
01020 L220:
01021 
01022 /*     If there are precomputed reciprocal condition numbers, compare */
01023 /*     computed values with them. */
01024 
01025     if (*comp) {
01026 
01027 /*        First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that */
01028 /*        the logical function ZSLECT selects the eigenvalues specified */
01029 /*        by NSLCT, ISLCT and ISRT. */
01030 
01031         sslct_1.seldim = *n;
01032         sslct_1.selopt = 1;
01033         eps = max(ulp,5.9605e-8);
01034         i__1 = *n;
01035         for (i__ = 1; i__ <= i__1; ++i__) {
01036             ipnt[i__ - 1] = i__;
01037             sslct_1.selval[i__ - 1] = FALSE_;
01038             i__2 = i__;
01039             sslct_1.selwr[i__ - 1] = wtmp[i__2].r;
01040             sslct_1.selwi[i__ - 1] = d_imag(&wtmp[i__]);
01041 /* L230: */
01042         }
01043         i__1 = *n - 1;
01044         for (i__ = 1; i__ <= i__1; ++i__) {
01045             kmin = i__;
01046             if (*isrt == 0) {
01047                 i__2 = i__;
01048                 vrimin = wtmp[i__2].r;
01049             } else {
01050                 vrimin = d_imag(&wtmp[i__]);
01051             }
01052             i__2 = *n;
01053             for (j = i__ + 1; j <= i__2; ++j) {
01054                 if (*isrt == 0) {
01055                     i__3 = j;
01056                     vricmp = wtmp[i__3].r;
01057                 } else {
01058                     vricmp = d_imag(&wtmp[j]);
01059                 }
01060                 if (vricmp < vrimin) {
01061                     kmin = j;
01062                     vrimin = vricmp;
01063                 }
01064 /* L240: */
01065             }
01066             i__2 = kmin;
01067             ctmp.r = wtmp[i__2].r, ctmp.i = wtmp[i__2].i;
01068             i__2 = kmin;
01069             i__3 = i__;
01070             wtmp[i__2].r = wtmp[i__3].r, wtmp[i__2].i = wtmp[i__3].i;
01071             i__2 = i__;
01072             wtmp[i__2].r = ctmp.r, wtmp[i__2].i = ctmp.i;
01073             itmp = ipnt[i__ - 1];
01074             ipnt[i__ - 1] = ipnt[kmin - 1];
01075             ipnt[kmin - 1] = itmp;
01076 /* L250: */
01077         }
01078         i__1 = *nslct;
01079         for (i__ = 1; i__ <= i__1; ++i__) {
01080             sslct_1.selval[ipnt[islct[i__] - 1] - 1] = TRUE_;
01081 /* L260: */
01082         }
01083 
01084 /*        Compute condition numbers */
01085 
01086         zlacpy_("F", n, n, &a[a_offset], lda, &ht[ht_offset], lda);
01087         zgeesx_("N", "S", (L_fp)zslect_, "B", n, &ht[ht_offset], lda, &sdim1, 
01088                 &wt[1], &vs1[vs1_offset], ldvs, &rconde, &rcondv, &work[1], 
01089                 lwork, &rwork[1], &bwork[1], &iinfo);
01090         if (iinfo != 0) {
01091             result[16] = ulpinv;
01092             result[17] = ulpinv;
01093             io___42.ciunit = *nounit;
01094             s_wsfe(&io___42);
01095             do_fio(&c__1, "ZGEESX9", (ftnlen)7);
01096             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01097             do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
01098             do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
01099             e_wsfe();
01100             *info = abs(iinfo);
01101             goto L270;
01102         }
01103 
01104 /*        Compare condition number for average of selected eigenvalues */
01105 /*        taking its condition number into account */
01106 
01107         anorm = zlange_("1", n, n, &a[a_offset], lda, &rwork[1]);
01108 /* Computing MAX */
01109         d__1 = (doublereal) (*n) * eps * anorm;
01110         v = max(d__1,smlnum);
01111         if (anorm == 0.) {
01112             v = 1.;
01113         }
01114         if (v > rcondv) {
01115             tol = 1.;
01116         } else {
01117             tol = v / rcondv;
01118         }
01119         if (v > *rcdvin) {
01120             tolin = 1.;
01121         } else {
01122             tolin = v / *rcdvin;
01123         }
01124 /* Computing MAX */
01125         d__1 = tol, d__2 = smlnum / eps;
01126         tol = max(d__1,d__2);
01127 /* Computing MAX */
01128         d__1 = tolin, d__2 = smlnum / eps;
01129         tolin = max(d__1,d__2);
01130         if (eps * (*rcdein - tolin) > rconde + tol) {
01131             result[16] = ulpinv;
01132         } else if (*rcdein - tolin > rconde + tol) {
01133             result[16] = (*rcdein - tolin) / (rconde + tol);
01134         } else if (*rcdein + tolin < eps * (rconde - tol)) {
01135             result[16] = ulpinv;
01136         } else if (*rcdein + tolin < rconde - tol) {
01137             result[16] = (rconde - tol) / (*rcdein + tolin);
01138         } else {
01139             result[16] = 1.;
01140         }
01141 
01142 /*        Compare condition numbers for right invariant subspace */
01143 /*        taking its condition number into account */
01144 
01145         if (v > rcondv * rconde) {
01146             tol = rcondv;
01147         } else {
01148             tol = v / rconde;
01149         }
01150         if (v > *rcdvin * *rcdein) {
01151             tolin = *rcdvin;
01152         } else {
01153             tolin = v / *rcdein;
01154         }
01155 /* Computing MAX */
01156         d__1 = tol, d__2 = smlnum / eps;
01157         tol = max(d__1,d__2);
01158 /* Computing MAX */
01159         d__1 = tolin, d__2 = smlnum / eps;
01160         tolin = max(d__1,d__2);
01161         if (eps * (*rcdvin - tolin) > rcondv + tol) {
01162             result[17] = ulpinv;
01163         } else if (*rcdvin - tolin > rcondv + tol) {
01164             result[17] = (*rcdvin - tolin) / (rcondv + tol);
01165         } else if (*rcdvin + tolin < eps * (rcondv - tol)) {
01166             result[17] = ulpinv;
01167         } else if (*rcdvin + tolin < rcondv - tol) {
01168             result[17] = (rcondv - tol) / (*rcdvin + tolin);
01169         } else {
01170             result[17] = 1.;
01171         }
01172 
01173 L270:
01174 
01175         ;
01176     }
01177 
01178 
01179     return 0;
01180 
01181 /*     End of ZGET24 */
01182 
01183 } /* zget24_ */


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