zget36.c

Go to the documentation of this file.

/* zget36.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 doublecomplex c_b1 = {0.,0.};
static doublecomplex c_b2 = {1.,0.};
static integer c__3 = 3;
static integer c__1 = 1;
static integer c__7 = 7;
static integer c__10 = 10;
static integer c__11 = 11;
static integer c__200 = 200;

/* Subroutine */ int zget36_(doublereal *rmax, integer *lmax, integer *ninfo, 
        integer *knt, integer *nin)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4;

    /* Builtin functions */
    integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
            e_rsle(void);

    /* Local variables */
    integer i__, j, n;
    doublecomplex q[100]        /* was [10][10] */, t1[100]     /* was [10][
            10] */, t2[100]     /* was [10][10] */;
    doublereal eps, res;
    doublecomplex tmp[100]      /* was [10][10] */, diag[10];
    integer ifst, ilst;
    doublecomplex work[200];
    integer info1, info2;
    doublecomplex ctemp;
    extern /* Subroutine */ int zhst01_(integer *, integer *, integer *, 
            doublecomplex *, integer *, doublecomplex *, integer *, 
            doublecomplex *, integer *, doublecomplex *, integer *, 
            doublereal *, doublereal *);
    doublereal rwork[10];
    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
            doublecomplex *, integer *);
    extern doublereal dlamch_(char *);
    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 
            doublecomplex *, integer *, doublecomplex *, integer *), 
            zlaset_(char *, integer *, integer *, doublecomplex *, 
            doublecomplex *, doublecomplex *, integer *);
    doublereal result[2];
    extern /* Subroutine */ int ztrexc_(char *, integer *, doublecomplex *, 
            integer *, doublecomplex *, integer *, integer *, integer *, 
            integer *);

    /* Fortran I/O blocks */
    static cilist io___2 = { 0, 0, 0, 0, 0 };
    static cilist io___7 = { 0, 0, 0, 0, 0 };



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

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

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

/*  ZGET36 tests ZTREXC, a routine for reordering diagonal entries of a */
/*  matrix in complex Schur form. Thus, ZLAEXC computes a unitary matrix */
/*  Q such that */

/*     Q' * T1 * Q  = T2 */

/*  and where one of the diagonal blocks of T1 (the one at row IFST) has */
/*  been moved to position ILST. */

/*  The test code verifies that the residual Q'*T1*Q-T2 is small, that T2 */
/*  is in Schur form, and that the final position of the IFST block is */
/*  ILST. */

/*  The test matrices are read from a file with logical unit number NIN. */

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

/*  RMAX    (output) DOUBLE PRECISION */
/*          Value of the largest test ratio. */

/*  LMAX    (output) INTEGER */
/*          Example number where largest test ratio achieved. */

/*  NINFO   (output) INTEGER */
/*          Number of examples where INFO is nonzero. */

/*  KNT     (output) INTEGER */
/*          Total number of examples tested. */

/*  NIN     (input) INTEGER */
/*          Input logical unit number. */

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

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

    eps = dlamch_("P");
    *rmax = 0.;
    *lmax = 0;
    *knt = 0;
    *ninfo = 0;

/*     Read input data until N=0 */

L10:
    io___2.ciunit = *nin;
    s_rsle(&io___2);
    do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&ifst, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&ilst, (ftnlen)sizeof(integer));
    e_rsle();
    if (n == 0) {
        return 0;
    }
    ++(*knt);
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        io___7.ciunit = *nin;
        s_rsle(&io___7);
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            do_lio(&c__7, &c__1, (char *)&tmp[i__ + j * 10 - 11], (ftnlen)
                    sizeof(doublecomplex));
        }
        e_rsle();
/* L20: */
    }
    zlacpy_("F", &n, &n, tmp, &c__10, t1, &c__10);
    zlacpy_("F", &n, &n, tmp, &c__10, t2, &c__10);
    res = 0.;

/*     Test without accumulating Q */

    zlaset_("Full", &n, &n, &c_b1, &c_b2, q, &c__10);
    ztrexc_("N", &n, t1, &c__10, q, &c__10, &ifst, &ilst, &info1);
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            i__3 = i__ + j * 10 - 11;
            if (i__ == j && (q[i__3].r != 1. || q[i__3].i != 0.)) {
                res += 1. / eps;
            }
            i__3 = i__ + j * 10 - 11;
            if (i__ != j && (q[i__3].r != 0. || q[i__3].i != 0.)) {
                res += 1. / eps;
            }
/* L30: */
        }
/* L40: */
    }

/*     Test with accumulating Q */

    zlaset_("Full", &n, &n, &c_b1, &c_b2, q, &c__10);
    ztrexc_("V", &n, t2, &c__10, q, &c__10, &ifst, &ilst, &info2);

/*     Compare T1 with T2 */

    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            i__3 = i__ + j * 10 - 11;
            i__4 = i__ + j * 10 - 11;
            if (t1[i__3].r != t2[i__4].r || t1[i__3].i != t2[i__4].i) {
                res += 1. / eps;
            }
/* L50: */
        }
/* L60: */
    }
    if (info1 != 0 || info2 != 0) {
        ++(*ninfo);
    }
    if (info1 != info2) {
        res += 1. / eps;
    }

/*     Test for successful reordering of T2 */

    zcopy_(&n, tmp, &c__11, diag, &c__1);
    if (ifst < ilst) {
        i__1 = ilst;
        for (i__ = ifst + 1; i__ <= i__1; ++i__) {
            i__2 = i__ - 1;
            ctemp.r = diag[i__2].r, ctemp.i = diag[i__2].i;
            i__2 = i__ - 1;
            i__3 = i__ - 2;
            diag[i__2].r = diag[i__3].r, diag[i__2].i = diag[i__3].i;
            i__2 = i__ - 2;
            diag[i__2].r = ctemp.r, diag[i__2].i = ctemp.i;
/* L70: */
        }
    } else if (ifst > ilst) {
        i__1 = ilst;
        for (i__ = ifst - 1; i__ >= i__1; --i__) {
            i__2 = i__;
            ctemp.r = diag[i__2].r, ctemp.i = diag[i__2].i;
            i__2 = i__;
            i__3 = i__ - 1;
            diag[i__2].r = diag[i__3].r, diag[i__2].i = diag[i__3].i;
            i__2 = i__ - 1;
            diag[i__2].r = ctemp.r, diag[i__2].i = ctemp.i;
/* L80: */
        }
    }
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = i__ + i__ * 10 - 11;
        i__3 = i__ - 1;
        if (t2[i__2].r != diag[i__3].r || t2[i__2].i != diag[i__3].i) {
            res += 1. / eps;
        }
/* L90: */
    }

/*     Test for small residual, and orthogonality of Q */

    zhst01_(&n, &c__1, &n, tmp, &c__10, t2, &c__10, q, &c__10, work, &c__200, 
            rwork, result);
    res = res + result[0] + result[1];

/*     Test for T2 being in Schur form */

    i__1 = n - 1;
    for (j = 1; j <= i__1; ++j) {
        i__2 = n;
        for (i__ = j + 1; i__ <= i__2; ++i__) {
            i__3 = i__ + j * 10 - 11;
            if (t2[i__3].r != 0. || t2[i__3].i != 0.) {
                res += 1. / eps;
            }
/* L100: */
        }
/* L110: */
    }
    if (res > *rmax) {
        *rmax = res;
        *lmax = *knt;
    }
    goto L10;

/*     End of ZGET36 */

} /* zget36_ */