#include "blaswrap.h"
/* dget35.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"

/* Table of constant values */

static integer c__6 = 6;
static doublereal c_b35 = 1.;

/* Subroutine */ int dget35_(doublereal *rmax, integer *lmax, integer *ninfo, 
	integer *knt)
{
    /* Initialized data */

    static integer idim[8] = { 1,2,3,4,3,3,6,4 };
    static integer ival[288]	/* was [6][6][8] */ = { 1,0,0,0,0,0,0,0,0,0,0,
	    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,2,0,0,0,0,-2,
	    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
	    0,0,5,1,2,0,0,0,-8,-2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
	    3,4,0,0,0,0,-5,3,0,0,0,0,1,2,1,4,0,0,-3,-9,-1,1,0,0,0,0,0,0,0,0,0,
	    0,0,0,0,0,1,0,0,0,0,0,2,3,0,0,0,0,5,6,7,0,0,0,0,0,0,0,0,0,0,0,0,0,
	    0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,3,-4,0,0,0,2,5,2,0,0,0,0,0,0,0,0,0,
	    0,0,0,0,0,0,0,0,0,0,0,0,1,2,0,0,0,0,-2,0,0,0,0,0,5,6,3,4,0,0,-1,
	    -9,-5,2,0,0,8,8,8,8,5,6,9,9,9,9,-7,5,1,0,0,0,0,0,1,5,2,0,0,0,2,
	    -21,5,0,0,0,1,2,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };

    /* System generated locals */
    integer i__1, i__2, i__3;
    doublereal d__1, d__2, d__3;

    /* Builtin functions */
    double sqrt(doublereal), sin(doublereal);

    /* Local variables */
    static doublereal a[36]	/* was [6][6] */, b[36]	/* was [6][6] */, c__[
	    36]	/* was [6][6] */;
    static integer i__, j, m, n;
    static doublereal cc[36]	/* was [6][6] */, vm1[3], vm2[3];
    static integer ima, imb;
    static doublereal dum[1], eps, res, res1;
    static integer info;
    static doublereal cnrm;
    static integer isgn;
    static doublereal rmul, tnrm, xnrm, scale;
    extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
	    integer *, doublereal *, doublereal *, integer *, doublereal *, 
	    integer *, doublereal *, doublereal *, integer *);
    static char trana[1], tranb[1];
    static integer imlda1, imlda2, imldb1;
    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
    extern doublereal dlamch_(char *), dlange_(char *, integer *, 
	    integer *, doublereal *, integer *, doublereal *);
    static integer imloff, itrana, itranb;
    static doublereal bignum, smlnum;
    extern /* Subroutine */ int dtrsyl_(char *, char *, integer *, integer *, 
	    integer *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, integer *, doublereal *, integer *);


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


    Purpose   
    =======   

    DGET35 tests DTRSYL, a routine for solving the Sylvester matrix   
    equation   

       op(A)*X + ISGN*X*op(B) = scale*C,   

    A and B are assumed to be in Schur canonical form, op() represents an   
    optional transpose, and ISGN can be -1 or +1.  Scale is an output   
    less than or equal to 1, chosen to avoid overflow in X.   

    The test code verifies that the following residual is order 1:   

       norm(op(A)*X + ISGN*X*op(B) - scale*C) /   
           (EPS*max(norm(A),norm(B))*norm(X))   

    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.   

    =====================================================================   


       Get machine parameters */

    eps = dlamch_("P");
    smlnum = dlamch_("S") * 4. / eps;
    bignum = 1. / smlnum;
    dlabad_(&smlnum, &bignum);

/*     Set up test case parameters */

    vm1[0] = sqrt(smlnum);
    vm1[1] = 1.;
    vm1[2] = sqrt(bignum);
    vm2[0] = 1.;
    vm2[1] = eps * 2. + 1.;
    vm2[2] = 2.;

    *knt = 0;
    *ninfo = 0;
    *lmax = 0;
    *rmax = 0.;

/*     Begin test loop */

    for (itrana = 1; itrana <= 2; ++itrana) {
	for (itranb = 1; itranb <= 2; ++itranb) {
	    for (isgn = -1; isgn <= 1; isgn += 2) {
		for (ima = 1; ima <= 8; ++ima) {
		    for (imlda1 = 1; imlda1 <= 3; ++imlda1) {
			for (imlda2 = 1; imlda2 <= 3; ++imlda2) {
			    for (imloff = 1; imloff <= 2; ++imloff) {
				for (imb = 1; imb <= 8; ++imb) {
				    for (imldb1 = 1; imldb1 <= 3; ++imldb1) {
					if (itrana == 1) {
					    *(unsigned char *)trana = 'N';
					}
					if (itrana == 2) {
					    *(unsigned char *)trana = 'T';
					}
					if (itranb == 1) {
					    *(unsigned char *)tranb = 'N';
					}
					if (itranb == 2) {
					    *(unsigned char *)tranb = 'T';
					}
					m = idim[ima - 1];
					n = idim[imb - 1];
					tnrm = 0.;
					i__1 = m;
					for (i__ = 1; i__ <= i__1; ++i__) {
					    i__2 = m;
					    for (j = 1; j <= i__2; ++j) {
			  a[i__ + j * 6 - 7] = (doublereal) ival[i__ + (j + 
				  ima * 6) * 6 - 43];
			  if ((i__3 = i__ - j, abs(i__3)) <= 1) {
			      a[i__ + j * 6 - 7] *= vm1[imlda1 - 1];
			      a[i__ + j * 6 - 7] *= vm2[imlda2 - 1];
			  } else {
			      a[i__ + j * 6 - 7] *= vm1[imloff - 1];
			  }
/* Computing MAX */
			  d__2 = tnrm, d__3 = (d__1 = a[i__ + j * 6 - 7], abs(
				  d__1));
			  tnrm = max(d__2,d__3);
/* L10: */
					    }
/* L20: */
					}
					i__1 = n;
					for (i__ = 1; i__ <= i__1; ++i__) {
					    i__2 = n;
					    for (j = 1; j <= i__2; ++j) {
			  b[i__ + j * 6 - 7] = (doublereal) ival[i__ + (j + 
				  imb * 6) * 6 - 43];
			  if ((i__3 = i__ - j, abs(i__3)) <= 1) {
			      b[i__ + j * 6 - 7] *= vm1[imldb1 - 1];
			  } else {
			      b[i__ + j * 6 - 7] *= vm1[imloff - 1];
			  }
/* Computing MAX */
			  d__2 = tnrm, d__3 = (d__1 = b[i__ + j * 6 - 7], abs(
				  d__1));
			  tnrm = max(d__2,d__3);
/* L30: */
					    }
/* L40: */
					}
					cnrm = 0.;
					i__1 = m;
					for (i__ = 1; i__ <= i__1; ++i__) {
					    i__2 = n;
					    for (j = 1; j <= i__2; ++j) {
			  c__[i__ + j * 6 - 7] = sin((doublereal) (i__ * j));
/* Computing MAX */
			  d__1 = cnrm, d__2 = c__[i__ + j * 6 - 7];
			  cnrm = max(d__1,d__2);
			  cc[i__ + j * 6 - 7] = c__[i__ + j * 6 - 7];
/* L50: */
					    }
/* L60: */
					}
					++(*knt);
					dtrsyl_(trana, tranb, &isgn, &m, &n, 
						a, &c__6, b, &c__6, c__, &
						c__6, &scale, &info);
					if (info != 0) {
					    ++(*ninfo);
					}
					xnrm = dlange_("M", &m, &n, c__, &
						c__6, dum);
					rmul = 1.;
					if (xnrm > 1. && tnrm > 1.) {
					    if (xnrm > bignum / tnrm) {
			  rmul = 1. / max(xnrm,tnrm);
					    }
					}
					d__1 = -scale * rmul;
					dgemm_(trana, "N", &m, &n, &m, &rmul, 
						a, &c__6, c__, &c__6, &d__1, 
						cc, &c__6);
					d__1 = (doublereal) isgn * rmul;
					dgemm_("N", tranb, &m, &n, &n, &d__1, 
						c__, &c__6, b, &c__6, &c_b35, 
						cc, &c__6);
					res1 = dlange_("M", &m, &n, cc, &c__6,
						 dum);
/* Computing MAX */
					d__1 = smlnum, d__2 = smlnum * xnrm, 
						d__1 = max(d__1,d__2), d__2 = 
						rmul * tnrm * eps * xnrm;
					res = res1 / max(d__1,d__2);
					if (res > *rmax) {
					    *lmax = *knt;
					    *rmax = res;
					}
/* L70: */
				    }
/* L80: */
				}
/* L90: */
			    }
/* L100: */
			}
/* L110: */
		    }
/* L120: */
		}
/* L130: */
	    }
/* L140: */
	}
/* L150: */
    }

    return 0;

/*     End of DGET35 */

} /* dget35_ */