ddrgvx.c

Go to the documentation of this file.

/* ddrgvx.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"
#include "blaswrap.h"

/* Table of constant values */

static integer c__1 = 1;
static integer c__0 = 0;
static integer c__5 = 5;
static logical c_true = TRUE_;
static logical c_false = FALSE_;
static integer c__3 = 3;

/* Subroutine */ int ddrgvx_(integer *nsize, doublereal *thresh, integer *nin, 
         integer *nout, doublereal *a, integer *lda, doublereal *b, 
        doublereal *ai, doublereal *bi, doublereal *alphar, doublereal *
        alphai, doublereal *beta, doublereal *vl, doublereal *vr, integer *
        ilo, integer *ihi, doublereal *lscale, doublereal *rscale, doublereal 
        *s, doublereal *dtru, doublereal *dif, doublereal *diftru, doublereal 
        *work, integer *lwork, integer *iwork, integer *liwork, doublereal *
        result, logical *bwork, integer *info)
{
    /* Format strings */
    static char fmt_9999[] = "(\002 DDRGVX: \002,a,\002 returned INFO=\002,i"
            "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002)\002)";
    static char fmt_9998[] = "(\002 DDRGVX: \002,a,\002 Eigenvectors from"
            " \002,a,\002 incorrectly \002,\002normalized.\002,/\002 Bits of "
            "error=\002,0p,g10.3,\002,\002,9x,\002N=\002,i6,\002, JTYPE=\002,"
            "i6,\002, IWA=\002,i5,\002, IWB=\002,i5,\002, IWX=\002,i5,\002, I"
            "WY=\002,i5)";
    static char fmt_9997[] = "(/1x,a3,\002 -- Real Expert Eigenvalue/vecto"
            "r\002,\002 problem driver\002)";
    static char fmt_9995[] = "(\002 Matrix types: \002,/)";
    static char fmt_9994[] = "(\002 TYPE 1: Da is diagonal, Db is identity,"
            " \002,/\002     A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) \002,/"
            "\002     YH and X are left and right eigenvectors. \002,/)";
    static char fmt_9993[] = "(\002 TYPE 2: Da is quasi-diagonal, Db is iden"
            "tity, \002,/\002     A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1)"
            " \002,/\002     YH and X are left and right eigenvectors. \002,/)"
            ;
    static char fmt_9992[] = "(/\002 Tests performed:  \002,/4x,\002 a is al"
            "pha, b is beta, l is a left eigenvector, \002,/4x,\002 r is a ri"
            "ght eigenvector and \002,a,\002 means \002,a,\002.\002,/\002 1 ="
            " max | ( b A - a B )\002,a,\002 l | / const.\002,/\002 2 = max |"
            " ( b A - a B ) r | / const.\002,/\002 3 = max ( Sest/Stru, Stru/"
            "Sest ) \002,\002 over all eigenvalues\002,/\002 4 = max( DIFest/"
            "DIFtru, DIFtru/DIFest ) \002,\002 over the 1st and 5th eigenvect"
            "ors\002,/)";
    static char fmt_9991[] = "(\002 Type=\002,i2,\002,\002,\002 IWA=\002,i2"
            ",\002, IWB=\002,i2,\002, IWX=\002,i2,\002, IWY=\002,i2,\002, res"
            "ult \002,i2,\002 is\002,0p,f8.2)";
    static char fmt_9990[] = "(\002 Type=\002,i2,\002,\002,\002 IWA=\002,i2"
            ",\002, IWB=\002,i2,\002, IWX=\002,i2,\002, IWY=\002,i2,\002, res"
            "ult \002,i2,\002 is\002,1p,d10.3)";
    static char fmt_9987[] = "(\002 DDRGVX: \002,a,\002 returned INFO=\002,i"
            "6,\002.\002,/9x,\002N=\002,i6,\002, Input example #\002,i2,\002"
            ")\002)";
    static char fmt_9986[] = "(\002 DDRGVX: \002,a,\002 Eigenvectors from"
            " \002,a,\002 incorrectly \002,\002normalized.\002,/\002 Bits of "
            "error=\002,0p,g10.3,\002,\002,9x,\002N=\002,i6,\002, Input Examp"
            "le #\002,i2,\002)\002)";
    static char fmt_9996[] = "(\002 Input Example\002)";
    static char fmt_9989[] = "(\002 Input example #\002,i2,\002, matrix orde"
            "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,0p,f8.2)";
    static char fmt_9988[] = "(\002 Input example #\002,i2,\002, matrix orde"
            "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,1p,d10.3)";

    /* System generated locals */
    integer a_dim1, a_offset, ai_dim1, ai_offset, b_dim1, b_offset, bi_dim1, 
            bi_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, i__2;
    doublereal d__1, d__2, d__3, d__4;

    /* Builtin functions */
    double sqrt(doublereal);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
             s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
            e_rsle(void);

    /* Local variables */
    integer i__, j, n, iwa, iwb;
    doublereal ulp;
    integer iwx, iwy, nmax;
    extern /* Subroutine */ int dget52_(logical *, integer *, doublereal *, 
            integer *, doublereal *, integer *, doublereal *, integer *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *);
    integer linfo;
    doublereal anorm, bnorm;
    integer nerrs;
    extern /* Subroutine */ int dlatm6_(integer *, integer *, doublereal *, 
            integer *, doublereal *, doublereal *, integer *, doublereal *, 
            integer *, doublereal *, doublereal *, doublereal *, doublereal *, 
             doublereal *, doublereal *);
    doublereal ratio1, ratio2, thrsh2;
    extern doublereal dlamch_(char *), dlange_(char *, integer *, 
            integer *, doublereal *, integer *, doublereal *);
    extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 
            doublereal *, integer *, doublereal *, integer *), 
            xerbla_(char *, integer *);
    doublereal abnorm;
    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
            integer *, integer *);
    extern /* Subroutine */ int alasvm_(char *, integer *, integer *, integer 
            *, integer *), dggevx_(char *, char *, char *, char *, 
            integer *, doublereal *, integer *, doublereal *, integer *, 
            doublereal *, doublereal *, doublereal *, doublereal *, integer *, 
             doublereal *, integer *, integer *, integer *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *, doublereal *, integer *, integer *, logical *, 
            integer *);
    doublereal weight[5];
    integer minwrk, maxwrk, iptype;
    doublereal ulpinv;
    integer nptknt, ntestt;

    /* Fortran I/O blocks */
    static cilist io___20 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___22 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___23 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___28 = { 0, 0, 0, fmt_9997, 0 };
    static cilist io___29 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___30 = { 0, 0, 0, fmt_9994, 0 };
    static cilist io___31 = { 0, 0, 0, fmt_9993, 0 };
    static cilist io___32 = { 0, 0, 0, fmt_9992, 0 };
    static cilist io___33 = { 0, 0, 0, fmt_9991, 0 };
    static cilist io___34 = { 0, 0, 0, fmt_9990, 0 };
    static cilist io___35 = { 0, 0, 1, 0, 0 };
    static cilist io___36 = { 0, 0, 0, 0, 0 };
    static cilist io___37 = { 0, 0, 0, 0, 0 };
    static cilist io___38 = { 0, 0, 0, 0, 0 };
    static cilist io___39 = { 0, 0, 0, 0, 0 };
    static cilist io___40 = { 0, 0, 0, fmt_9987, 0 };
    static cilist io___41 = { 0, 0, 0, fmt_9986, 0 };
    static cilist io___42 = { 0, 0, 0, fmt_9986, 0 };
    static cilist io___43 = { 0, 0, 0, fmt_9997, 0 };
    static cilist io___44 = { 0, 0, 0, fmt_9996, 0 };
    static cilist io___45 = { 0, 0, 0, fmt_9992, 0 };
    static cilist io___46 = { 0, 0, 0, fmt_9989, 0 };
    static cilist io___47 = { 0, 0, 0, fmt_9988, 0 };



/*  -- LAPACK test routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  DDRGVX checks the nonsymmetric generalized eigenvalue problem */
/*  expert driver DGGEVX. */

/*  DGGEVX computes the generalized eigenvalues, (optionally) the left */
/*  and/or right eigenvectors, (optionally) computes a balancing */
/*  transformation to improve the conditioning, and (optionally) */
/*  reciprocal condition numbers for the eigenvalues and eigenvectors. */

/*  When DDRGVX is called with NSIZE > 0, two types of test matrix pairs */
/*  are generated by the subroutine DLATM6 and test the driver DGGEVX. */
/*  The test matrices have the known exact condition numbers for */
/*  eigenvalues. For the condition numbers of the eigenvectors */
/*  corresponding the first and last eigenvalues are also know */
/*  ``exactly'' (see DLATM6). */

/*  For each matrix pair, the following tests will be performed and */
/*  compared with the threshhold THRESH. */

/*  (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of */

/*     | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) ) */

/*      where l**H is the conjugate tranpose of l. */

/*  (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of */

/*        | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) ) */

/*  (3) The condition number S(i) of eigenvalues computed by DGGEVX */
/*      differs less than a factor THRESH from the exact S(i) (see */
/*      DLATM6). */

/*  (4) DIF(i) computed by DTGSNA differs less than a factor 10*THRESH */
/*      from the exact value (for the 1st and 5th vectors only). */

/*  Test Matrices */
/*  ============= */

/*  Two kinds of test matrix pairs */

/*           (A, B) = inverse(YH) * (Da, Db) * inverse(X) */

/*  are used in the tests: */

/*  1: Da = 1+a   0    0    0    0    Db = 1   0   0   0   0 */
/*           0   2+a   0    0    0         0   1   0   0   0 */
/*           0    0   3+a   0    0         0   0   1   0   0 */
/*           0    0    0   4+a   0         0   0   0   1   0 */
/*           0    0    0    0   5+a ,      0   0   0   0   1 , and */

/*  2: Da =  1   -1    0    0    0    Db = 1   0   0   0   0 */
/*           1    1    0    0    0         0   1   0   0   0 */
/*           0    0    1    0    0         0   0   1   0   0 */
/*           0    0    0   1+a  1+b        0   0   0   1   0 */
/*           0    0    0  -1-b  1+a ,      0   0   0   0   1 . */

/*  In both cases the same inverse(YH) and inverse(X) are used to compute */
/*  (A, B), giving the exact eigenvectors to (A,B) as (YH, X): */

/*  YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x */
/*          0    1   -y    y   -y         0   1   x  -x  -x */
/*          0    0    1    0    0         0   0   1   0   0 */
/*          0    0    0    1    0         0   0   0   1   0 */
/*          0    0    0    0    1,        0   0   0   0   1 , where */

/*  a, b, x and y will have all values independently of each other from */
/*  { sqrt(sqrt(ULP)),  0.1,  1,  10,  1/sqrt(sqrt(ULP)) }. */

/*  Arguments */
/*  ========= */

/*  NSIZE   (input) INTEGER */
/*          The number of sizes of matrices to use.  NSIZE must be at */
/*          least zero. If it is zero, no randomly generated matrices */
/*          are tested, but any test matrices read from NIN will be */
/*          tested. */

/*  THRESH  (input) DOUBLE PRECISION */
/*          A test will count as "failed" if the "error", computed as */
/*          described above, exceeds THRESH.  Note that the error */
/*          is scaled to be O(1), so THRESH should be a reasonably */
/*          small multiple of 1, e.g., 10 or 100.  In particular, */
/*          it should not depend on the precision (single vs. double) */
/*          or the size of the matrix.  It must be at least zero. */

/*  NIN     (input) INTEGER */
/*          The FORTRAN unit number for reading in the data file of */
/*          problems to solve. */

/*  NOUT    (input) INTEGER */
/*          The FORTRAN unit number for printing out error messages */
/*          (e.g., if a routine returns IINFO not equal to 0.) */

/*  A       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE) */
/*          Used to hold the matrix whose eigenvalues are to be */
/*          computed.  On exit, A contains the last matrix actually used. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of A, B, AI, BI, Ao, and Bo. */
/*          It must be at least 1 and at least NSIZE. */

/*  B       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE) */
/*          Used to hold the matrix whose eigenvalues are to be */
/*          computed.  On exit, B contains the last matrix actually used. */

/*  AI      (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE) */
/*          Copy of A, modified by DGGEVX. */

/*  BI      (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE) */
/*          Copy of B, modified by DGGEVX. */

/*  ALPHAR  (workspace) DOUBLE PRECISION array, dimension (NSIZE) */
/*  ALPHAI  (workspace) DOUBLE PRECISION array, dimension (NSIZE) */
/*  BETA    (workspace) DOUBLE PRECISION array, dimension (NSIZE) */
/*          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues. */

/*  VL      (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE) */
/*          VL holds the left eigenvectors computed by DGGEVX. */

/*  VR      (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE) */
/*          VR holds the right eigenvectors computed by DGGEVX. */

/*  ILO     (output/workspace) INTEGER */

/*  IHI     (output/workspace) INTEGER */

/*  LSCALE  (output/workspace) DOUBLE PRECISION array, dimension (N) */

/*  RSCALE  (output/workspace) DOUBLE PRECISION array, dimension (N) */

/*  S       (output/workspace) DOUBLE PRECISION array, dimension (N) */

/*  DTRU    (output/workspace) DOUBLE PRECISION array, dimension (N) */

/*  DIF     (output/workspace) DOUBLE PRECISION array, dimension (N) */

/*  DIFTRU  (output/workspace) DOUBLE PRECISION array, dimension (N) */

/*  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK) */

/*  LWORK   (input) INTEGER */
/*          Leading dimension of WORK.  LWORK >= 2*N*N+12*N+16. */

/*  IWORK   (workspace) INTEGER array, dimension (LIWORK) */

/*  LIWORK  (input) INTEGER */
/*          Leading dimension of IWORK.  Must be at least N+6. */

/*  RESULT  (output/workspace) DOUBLE PRECISION array, dimension (4) */

/*  BWORK   (workspace) LOGICAL array, dimension (N) */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
/*          > 0:  A routine returned an error code. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Local Arrays .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Check for errors */

    /* Parameter adjustments */
    vr_dim1 = *lda;
    vr_offset = 1 + vr_dim1;
    vr -= vr_offset;
    vl_dim1 = *lda;
    vl_offset = 1 + vl_dim1;
    vl -= vl_offset;
    bi_dim1 = *lda;
    bi_offset = 1 + bi_dim1;
    bi -= bi_offset;
    ai_dim1 = *lda;
    ai_offset = 1 + ai_dim1;
    ai -= ai_offset;
    b_dim1 = *lda;
    b_offset = 1 + b_dim1;
    b -= b_offset;
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --alphar;
    --alphai;
    --beta;
    --lscale;
    --rscale;
    --s;
    --dtru;
    --dif;
    --diftru;
    --work;
    --iwork;
    --result;
    --bwork;

    /* Function Body */
    *info = 0;

    nmax = 5;

    if (*nsize < 0) {
        *info = -1;
    } else if (*thresh < 0.) {
        *info = -2;
    } else if (*nin <= 0) {
        *info = -3;
    } else if (*nout <= 0) {
        *info = -4;
    } else if (*lda < 1 || *lda < nmax) {
        *info = -6;
    } else if (*liwork < nmax + 6) {
        *info = -26;
    }

/*     Compute workspace */
/*      (Note: Comments in the code beginning "Workspace:" describe the */
/*       minimal amount of workspace needed at that point in the code, */
/*       as well as the preferred amount for good performance. */
/*       NB refers to the optimal block size for the immediately */
/*       following subroutine, as returned by ILAENV.) */

    minwrk = 1;
    if (*info == 0 && *lwork >= 1) {
        minwrk = (nmax << 1) * nmax + nmax * 12 + 16;
        maxwrk = nmax * 6 + nmax * ilaenv_(&c__1, "DGEQRF", " ", &nmax, &c__1, 
                 &nmax, &c__0);
/* Computing MAX */
        i__1 = maxwrk, i__2 = (nmax << 1) * nmax + nmax * 12 + 16;
        maxwrk = max(i__1,i__2);
        work[1] = (doublereal) maxwrk;
    }

    if (*lwork < minwrk) {
        *info = -24;
    }

    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("DDRGVX", &i__1);
        return 0;
    }

    n = 5;
    ulp = dlamch_("P");
    ulpinv = 1. / ulp;
    thrsh2 = *thresh * 10.;
    nerrs = 0;
    nptknt = 0;
    ntestt = 0;

    if (*nsize == 0) {
        goto L90;
    }

/*     Parameters used for generating test matrices. */

    weight[0] = sqrt(sqrt(ulp));
    weight[1] = .1;
    weight[2] = 1.;
    weight[3] = 1. / weight[1];
    weight[4] = 1. / weight[0];

    for (iptype = 1; iptype <= 2; ++iptype) {
        for (iwa = 1; iwa <= 5; ++iwa) {
            for (iwb = 1; iwb <= 5; ++iwb) {
                for (iwx = 1; iwx <= 5; ++iwx) {
                    for (iwy = 1; iwy <= 5; ++iwy) {

/*                    generated a test matrix pair */

                        dlatm6_(&iptype, &c__5, &a[a_offset], lda, &b[
                                b_offset], &vr[vr_offset], lda, &vl[vl_offset]
, lda, &weight[iwa - 1], &weight[iwb - 1], &
                                weight[iwx - 1], &weight[iwy - 1], &dtru[1], &
                                diftru[1]);

/*                    Compute eigenvalues/eigenvectors of (A, B). */
/*                    Compute eigenvalue/eigenvector condition numbers */
/*                    using computed eigenvectors. */

                        dlacpy_("F", &n, &n, &a[a_offset], lda, &ai[ai_offset]
, lda);
                        dlacpy_("F", &n, &n, &b[b_offset], lda, &bi[bi_offset]
, lda);

                        dggevx_("N", "V", "V", "B", &n, &ai[ai_offset], lda, &
                                bi[bi_offset], lda, &alphar[1], &alphai[1], &
                                beta[1], &vl[vl_offset], lda, &vr[vr_offset], 
                                lda, ilo, ihi, &lscale[1], &rscale[1], &anorm, 
                                 &bnorm, &s[1], &dif[1], &work[1], lwork, &
                                iwork[1], &bwork[1], &linfo);
                        if (linfo != 0) {
                            result[1] = ulpinv;
                            io___20.ciunit = *nout;
                            s_wsfe(&io___20);
                            do_fio(&c__1, "DGGEVX", (ftnlen)6);
                            do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
                                    ;
                            do_fio(&c__1, (char *)&iptype, (ftnlen)sizeof(
                                    integer));
                            e_wsfe();
                            goto L30;
                        }

/*                    Compute the norm(A, B) */

                        dlacpy_("Full", &n, &n, &ai[ai_offset], lda, &work[1], 
                                 &n);
                        dlacpy_("Full", &n, &n, &bi[bi_offset], lda, &work[n *
                                 n + 1], &n);
                        i__1 = n << 1;
                        abnorm = dlange_("Fro", &n, &i__1, &work[1], &n, &
                                work[1]);

/*                    Tests (1) and (2) */

                        result[1] = 0.;
                        dget52_(&c_true, &n, &a[a_offset], lda, &b[b_offset], 
                                lda, &vl[vl_offset], lda, &alphar[1], &alphai[
                                1], &beta[1], &work[1], &result[1]);
                        if (result[2] > *thresh) {
                            io___22.ciunit = *nout;
                            s_wsfe(&io___22);
                            do_fio(&c__1, "Left", (ftnlen)4);
                            do_fio(&c__1, "DGGEVX", (ftnlen)6);
                            do_fio(&c__1, (char *)&result[2], (ftnlen)sizeof(
                                    doublereal));
                            do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
                                    ;
                            do_fio(&c__1, (char *)&iptype, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwa, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwb, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwx, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwy, (ftnlen)sizeof(
                                    integer));
                            e_wsfe();
                        }

                        result[2] = 0.;
                        dget52_(&c_false, &n, &a[a_offset], lda, &b[b_offset], 
                                 lda, &vr[vr_offset], lda, &alphar[1], &
                                alphai[1], &beta[1], &work[1], &result[2]);
                        if (result[3] > *thresh) {
                            io___23.ciunit = *nout;
                            s_wsfe(&io___23);
                            do_fio(&c__1, "Right", (ftnlen)5);
                            do_fio(&c__1, "DGGEVX", (ftnlen)6);
                            do_fio(&c__1, (char *)&result[3], (ftnlen)sizeof(
                                    doublereal));
                            do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
                                    ;
                            do_fio(&c__1, (char *)&iptype, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwa, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwb, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwx, (ftnlen)sizeof(
                                    integer));
                            do_fio(&c__1, (char *)&iwy, (ftnlen)sizeof(
                                    integer));
                            e_wsfe();
                        }

/*                    Test (3) */

                        result[3] = 0.;
                        i__1 = n;
                        for (i__ = 1; i__ <= i__1; ++i__) {
                            if (s[i__] == 0.) {
                                if (dtru[i__] > abnorm * ulp) {
                                    result[3] = ulpinv;
                                }
                            } else if (dtru[i__] == 0.) {
                                if (s[i__] > abnorm * ulp) {
                                    result[3] = ulpinv;
                                }
                            } else {
/* Computing MAX */
                                d__3 = (d__1 = dtru[i__] / s[i__], abs(d__1)),
                                         d__4 = (d__2 = s[i__] / dtru[i__], 
                                        abs(d__2));
                                work[i__] = max(d__3,d__4);
/* Computing MAX */
                                d__1 = result[3], d__2 = work[i__];
                                result[3] = max(d__1,d__2);
                            }
/* L10: */
                        }

/*                    Test (4) */

                        result[4] = 0.;
                        if (dif[1] == 0.) {
                            if (diftru[1] > abnorm * ulp) {
                                result[4] = ulpinv;
                            }
                        } else if (diftru[1] == 0.) {
                            if (dif[1] > abnorm * ulp) {
                                result[4] = ulpinv;
                            }
                        } else if (dif[5] == 0.) {
                            if (diftru[5] > abnorm * ulp) {
                                result[4] = ulpinv;
                            }
                        } else if (diftru[5] == 0.) {
                            if (dif[5] > abnorm * ulp) {
                                result[4] = ulpinv;
                            }
                        } else {
/* Computing MAX */
                            d__3 = (d__1 = diftru[1] / dif[1], abs(d__1)), 
                                    d__4 = (d__2 = dif[1] / diftru[1], abs(
                                    d__2));
                            ratio1 = max(d__3,d__4);
/* Computing MAX */
                            d__3 = (d__1 = diftru[5] / dif[5], abs(d__1)), 
                                    d__4 = (d__2 = dif[5] / diftru[5], abs(
                                    d__2));
                            ratio2 = max(d__3,d__4);
                            result[4] = max(ratio1,ratio2);
                        }

                        ntestt += 4;

/*                    Print out tests which fail. */

                        for (j = 1; j <= 4; ++j) {
                            if (result[j] >= thrsh2 && j >= 4 || result[j] >= 
                                    *thresh && j <= 3) {

/*                       If this is the first test to fail, */
/*                       print a header to the data file. */

                                if (nerrs == 0) {
                                    io___28.ciunit = *nout;
                                    s_wsfe(&io___28);
                                    do_fio(&c__1, "DXV", (ftnlen)3);
                                    e_wsfe();

/*                          Print out messages for built-in examples */

/*                          Matrix types */

                                    io___29.ciunit = *nout;
                                    s_wsfe(&io___29);
                                    e_wsfe();
                                    io___30.ciunit = *nout;
                                    s_wsfe(&io___30);
                                    e_wsfe();
                                    io___31.ciunit = *nout;
                                    s_wsfe(&io___31);
                                    e_wsfe();

/*                          Tests performed */

                                    io___32.ciunit = *nout;
                                    s_wsfe(&io___32);
                                    do_fio(&c__1, "'", (ftnlen)1);
                                    do_fio(&c__1, "transpose", (ftnlen)9);
                                    do_fio(&c__1, "'", (ftnlen)1);
                                    e_wsfe();

                                }
                                ++nerrs;
                                if (result[j] < 1e4) {
                                    io___33.ciunit = *nout;
                                    s_wsfe(&io___33);
                                    do_fio(&c__1, (char *)&iptype, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwa, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwb, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwx, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwy, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
                                            integer));
                                    do_fio(&c__1, (char *)&result[j], (ftnlen)
                                            sizeof(doublereal));
                                    e_wsfe();
                                } else {
                                    io___34.ciunit = *nout;
                                    s_wsfe(&io___34);
                                    do_fio(&c__1, (char *)&iptype, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwa, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwb, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwx, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&iwy, (ftnlen)
                                            sizeof(integer));
                                    do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
                                            integer));
                                    do_fio(&c__1, (char *)&result[j], (ftnlen)
                                            sizeof(doublereal));
                                    e_wsfe();
                                }
                            }
/* L20: */
                        }

L30:

/* L40: */
                        ;
                    }
/* L50: */
                }
/* L60: */
            }
/* L70: */
        }
/* L80: */
    }

    goto L150;

L90:

/*     Read in data from file to check accuracy of condition estimation */
/*     Read input data until N=0 */

    io___35.ciunit = *nin;
    i__1 = s_rsle(&io___35);
    if (i__1 != 0) {
        goto L150;
    }
    i__1 = do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
    if (i__1 != 0) {
        goto L150;
    }
    i__1 = e_rsle();
    if (i__1 != 0) {
        goto L150;
    }
    if (n == 0) {
        goto L150;
    }
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        io___36.ciunit = *nin;
        s_rsle(&io___36);
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            do_lio(&c__5, &c__1, (char *)&a[i__ + j * a_dim1], (ftnlen)sizeof(
                    doublereal));
        }
        e_rsle();
/* L100: */
    }
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        io___37.ciunit = *nin;
        s_rsle(&io___37);
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            do_lio(&c__5, &c__1, (char *)&b[i__ + j * b_dim1], (ftnlen)sizeof(
                    doublereal));
        }
        e_rsle();
/* L110: */
    }
    io___38.ciunit = *nin;
    s_rsle(&io___38);
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        do_lio(&c__5, &c__1, (char *)&dtru[i__], (ftnlen)sizeof(doublereal));
    }
    e_rsle();
    io___39.ciunit = *nin;
    s_rsle(&io___39);
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        do_lio(&c__5, &c__1, (char *)&diftru[i__], (ftnlen)sizeof(doublereal))
                ;
    }
    e_rsle();

    ++nptknt;

/*     Compute eigenvalues/eigenvectors of (A, B). */
/*     Compute eigenvalue/eigenvector condition numbers */
/*     using computed eigenvectors. */

    dlacpy_("F", &n, &n, &a[a_offset], lda, &ai[ai_offset], lda);
    dlacpy_("F", &n, &n, &b[b_offset], lda, &bi[bi_offset], lda);

    dggevx_("N", "V", "V", "B", &n, &ai[ai_offset], lda, &bi[bi_offset], lda, 
            &alphar[1], &alphai[1], &beta[1], &vl[vl_offset], lda, &vr[
            vr_offset], lda, ilo, ihi, &lscale[1], &rscale[1], &anorm, &bnorm, 
             &s[1], &dif[1], &work[1], lwork, &iwork[1], &bwork[1], &linfo);

    if (linfo != 0) {
        result[1] = ulpinv;
        io___40.ciunit = *nout;
        s_wsfe(&io___40);
        do_fio(&c__1, "DGGEVX", (ftnlen)6);
        do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(integer));
        do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
        do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
        e_wsfe();
        goto L140;
    }

/*     Compute the norm(A, B) */

    dlacpy_("Full", &n, &n, &ai[ai_offset], lda, &work[1], &n);
    dlacpy_("Full", &n, &n, &bi[bi_offset], lda, &work[n * n + 1], &n);
    i__1 = n << 1;
    abnorm = dlange_("Fro", &n, &i__1, &work[1], &n, &work[1]);

/*     Tests (1) and (2) */

    result[1] = 0.;
    dget52_(&c_true, &n, &a[a_offset], lda, &b[b_offset], lda, &vl[vl_offset], 
             lda, &alphar[1], &alphai[1], &beta[1], &work[1], &result[1]);
    if (result[2] > *thresh) {
        io___41.ciunit = *nout;
        s_wsfe(&io___41);
        do_fio(&c__1, "Left", (ftnlen)4);
        do_fio(&c__1, "DGGEVX", (ftnlen)6);
        do_fio(&c__1, (char *)&result[2], (ftnlen)sizeof(doublereal));
        do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
        do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
        e_wsfe();
    }

    result[2] = 0.;
    dget52_(&c_false, &n, &a[a_offset], lda, &b[b_offset], lda, &vr[vr_offset]
, lda, &alphar[1], &alphai[1], &beta[1], &work[1], &result[2]);
    if (result[3] > *thresh) {
        io___42.ciunit = *nout;
        s_wsfe(&io___42);
        do_fio(&c__1, "Right", (ftnlen)5);
        do_fio(&c__1, "DGGEVX", (ftnlen)6);
        do_fio(&c__1, (char *)&result[3], (ftnlen)sizeof(doublereal));
        do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
        do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
        e_wsfe();
    }

/*     Test (3) */

    result[3] = 0.;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        if (s[i__] == 0.) {
            if (dtru[i__] > abnorm * ulp) {
                result[3] = ulpinv;
            }
        } else if (dtru[i__] == 0.) {
            if (s[i__] > abnorm * ulp) {
                result[3] = ulpinv;
            }
        } else {
/* Computing MAX */
            d__3 = (d__1 = dtru[i__] / s[i__], abs(d__1)), d__4 = (d__2 = s[
                    i__] / dtru[i__], abs(d__2));
            work[i__] = max(d__3,d__4);
/* Computing MAX */
            d__1 = result[3], d__2 = work[i__];
            result[3] = max(d__1,d__2);
        }
/* L120: */
    }

/*     Test (4) */

    result[4] = 0.;
    if (dif[1] == 0.) {
        if (diftru[1] > abnorm * ulp) {
            result[4] = ulpinv;
        }
    } else if (diftru[1] == 0.) {
        if (dif[1] > abnorm * ulp) {
            result[4] = ulpinv;
        }
    } else if (dif[5] == 0.) {
        if (diftru[5] > abnorm * ulp) {
            result[4] = ulpinv;
        }
    } else if (diftru[5] == 0.) {
        if (dif[5] > abnorm * ulp) {
            result[4] = ulpinv;
        }
    } else {
/* Computing MAX */
        d__3 = (d__1 = diftru[1] / dif[1], abs(d__1)), d__4 = (d__2 = dif[1] /
                 diftru[1], abs(d__2));
        ratio1 = max(d__3,d__4);
/* Computing MAX */
        d__3 = (d__1 = diftru[5] / dif[5], abs(d__1)), d__4 = (d__2 = dif[5] /
                 diftru[5], abs(d__2));
        ratio2 = max(d__3,d__4);
        result[4] = max(ratio1,ratio2);
    }

    ntestt += 4;

/*     Print out tests which fail. */

    for (j = 1; j <= 4; ++j) {
        if (result[j] >= thrsh2) {

/*           If this is the first test to fail, */
/*           print a header to the data file. */

            if (nerrs == 0) {
                io___43.ciunit = *nout;
                s_wsfe(&io___43);
                do_fio(&c__1, "DXV", (ftnlen)3);
                e_wsfe();

/*              Print out messages for built-in examples */

/*              Matrix types */

                io___44.ciunit = *nout;
                s_wsfe(&io___44);
                e_wsfe();

/*              Tests performed */

                io___45.ciunit = *nout;
                s_wsfe(&io___45);
                do_fio(&c__1, "'", (ftnlen)1);
                do_fio(&c__1, "transpose", (ftnlen)9);
                do_fio(&c__1, "'", (ftnlen)1);
                e_wsfe();

            }
            ++nerrs;
            if (result[j] < 1e4) {
                io___46.ciunit = *nout;
                s_wsfe(&io___46);
                do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
                do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
                do_fio(&c__1, (char *)&result[j], (ftnlen)sizeof(doublereal));
                e_wsfe();
            } else {
                io___47.ciunit = *nout;
                s_wsfe(&io___47);
                do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
                do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
                do_fio(&c__1, (char *)&result[j], (ftnlen)sizeof(doublereal));
                e_wsfe();
            }
        }
/* L130: */
    }

L140:

    goto L90;
L150:

/*     Summary */

    alasvm_("DXV", nout, &nerrs, &ntestt, &c__0);

    work[1] = (doublereal) maxwrk;

    return 0;












/*     End of DDRGVX */

} /* ddrgvx_ */