#include "blaswrap.h"
/* dget32.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__2 = 2;

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

    static integer itval[32]	/* was [2][2][8] */ = { 8,4,2,1,4,8,1,2,2,1,8,
	    4,1,2,4,8,9,4,2,1,4,9,1,2,2,1,9,4,1,2,4,9 };

    /* System generated locals */
    doublereal d__1, d__2;

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

    /* Local variables */
    static doublereal b[4]	/* was [2][2] */, x[4]	/* was [2][2] */;
    static integer n1, n2, ib;
    static doublereal tl[4]	/* was [2][2] */, tr[4]	/* was [2][2] */;
    static integer ib1, ib2, ib3;
    static doublereal den, val[3], eps;
    static integer itl;
    static doublereal res, sgn;
    static integer itr;
    static doublereal tmp;
    static integer info, isgn;
    static doublereal tnrm, xnrm, scale, xnorm;
    extern /* Subroutine */ int dlasy2_(logical *, logical *, integer *, 
	    integer *, integer *, doublereal *, integer *, doublereal *, 
	    integer *, doublereal *, integer *, doublereal *, doublereal *, 
	    integer *, doublereal *, integer *), dlabad_(doublereal *, 
	    doublereal *);
    extern doublereal dlamch_(char *);
    static doublereal bignum;
    static integer itranl, itlscl;
    static logical ltranl;
    static integer itranr, itrscl;
    static logical ltranr;
    static doublereal smlnum;


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


    Purpose   
    =======   

    DGET32 tests DLASY2, a routine for solving   

            op(TL)*X + ISGN*X*op(TR) = SCALE*B   

    where TL is N1 by N1, TR is N2 by N2, and N1,N2 =1 or 2 only.   
    X and B are N1 by N2, op() is an optional transpose, an   
    ISGN = 1 or -1. SCALE is chosen less than or equal to 1 to   
    avoid overflow in X.   

    The test condition is that the scaled residual   

    norm( op(TL)*X + ISGN*X*op(TR) = SCALE*B )   
         / ( max( ulp*norm(TL), ulp*norm(TR)) * norm(X), SMLNUM )   

    should be on the order of 1. Here, ulp is the machine precision.   
    Also, it is verified that SCALE is less than or equal to 1, and   
    that XNORM = infinity-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 returned with INFO.NE.0.   

    KNT     (output) INTEGER   
            Total number of examples tested.   

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


       Get machine parameters */

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

/*     Set up test case parameters */

    val[0] = sqrt(smlnum);
    val[1] = 1.;
    val[2] = sqrt(bignum);

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

/*     Begin test loop */

    for (itranl = 0; itranl <= 1; ++itranl) {
	for (itranr = 0; itranr <= 1; ++itranr) {
	    for (isgn = -1; isgn <= 1; isgn += 2) {
		sgn = (doublereal) isgn;
		ltranl = itranl == 1;
		ltranr = itranr == 1;

		n1 = 1;
		n2 = 1;
		for (itl = 1; itl <= 3; ++itl) {
		    for (itr = 1; itr <= 3; ++itr) {
			for (ib = 1; ib <= 3; ++ib) {
			    tl[0] = val[itl - 1];
			    tr[0] = val[itr - 1];
			    b[0] = val[ib - 1];
			    ++(*knt);
			    dlasy2_(&ltranl, &ltranr, &isgn, &n1, &n2, tl, &
				    c__2, tr, &c__2, b, &c__2, &scale, x, &
				    c__2, &xnorm, &info);
			    if (info != 0) {
				++(*ninfo);
			    }
			    res = (d__1 = (tl[0] + sgn * tr[0]) * x[0] - 
				    scale * b[0], abs(d__1));
			    if (info == 0) {
/* Computing MAX */
				d__1 = eps * ((abs(tr[0]) + abs(tl[0])) * abs(
					x[0]));
				den = max(d__1,smlnum);
			    } else {
/* Computing MAX */
				d__1 = abs(x[0]);
				den = smlnum * max(d__1,1.);
			    }
			    res /= den;
			    if (scale > 1.) {
				res += 1. / eps;
			    }
			    res += (d__1 = xnorm - abs(x[0]), abs(d__1)) / 
				    max(smlnum,xnorm) / eps;
			    if (info != 0 && info != 1) {
				res += 1. / eps;
			    }
			    if (res > *rmax) {
				*lmax = *knt;
				*rmax = res;
			    }
/* L10: */
			}
/* L20: */
		    }
/* L30: */
		}

		n1 = 2;
		n2 = 1;
		for (itl = 1; itl <= 8; ++itl) {
		    for (itlscl = 1; itlscl <= 3; ++itlscl) {
			for (itr = 1; itr <= 3; ++itr) {
			    for (ib1 = 1; ib1 <= 3; ++ib1) {
				for (ib2 = 1; ib2 <= 3; ++ib2) {
				    b[0] = val[ib1 - 1];
				    b[1] = val[ib2 - 1] * -4.;
				    tl[0] = itval[((itl << 1) + 1 << 1) - 6] *
					     val[itlscl - 1];
				    tl[1] = itval[((itl << 1) + 1 << 1) - 5] *
					     val[itlscl - 1];
				    tl[2] = itval[((itl << 1) + 2 << 1) - 6] *
					     val[itlscl - 1];
				    tl[3] = itval[((itl << 1) + 2 << 1) - 5] *
					     val[itlscl - 1];
				    tr[0] = val[itr - 1];
				    ++(*knt);
				    dlasy2_(&ltranl, &ltranr, &isgn, &n1, &n2,
					     tl, &c__2, tr, &c__2, b, &c__2, &
					    scale, x, &c__2, &xnorm, &info);
				    if (info != 0) {
					++(*ninfo);
				    }
				    if (ltranl) {
					tmp = tl[2];
					tl[2] = tl[1];
					tl[1] = tmp;
				    }
				    res = (d__1 = (tl[0] + sgn * tr[0]) * x[0]
					     + tl[2] * x[1] - scale * b[0], 
					    abs(d__1));
				    res += (d__1 = (tl[3] + sgn * tr[0]) * x[
					    1] + tl[1] * x[0] - scale * b[1], 
					    abs(d__1)); 				    tnrm = abs(tr[0]) + abs(tl[0]) + abs(tl[2] 					    ) + abs(tl[1]) + abs(tl[3]);
/* Computing MAX */
				    d__1 = abs(x[0]), d__2 = abs(x[1]);
				    xnrm = max(d__1,d__2);
/* Computing MAX */
				    d__1 = smlnum, d__2 = smlnum * xnrm, d__1 
					    = max(d__1,d__2), d__2 = tnrm * 
					    eps * xnrm;
				    den = max(d__1,d__2);
				    res /= den;
				    if (scale > 1.) {
					res += 1. / eps;
				    }
				    res += (d__1 = xnorm - xnrm, abs(d__1)) / 
					    max(smlnum,xnorm) / eps;
				    if (res > *rmax) {
					*lmax = *knt;
					*rmax = res;
				    }
/* L40: */
				}
/* L50: */
			    }
/* L60: */
			}
/* L70: */
		    }
/* L80: */
		}

		n1 = 1;
		n2 = 2;
		for (itr = 1; itr <= 8; ++itr) {
		    for (itrscl = 1; itrscl <= 3; ++itrscl) {
			for (itl = 1; itl <= 3; ++itl) {
			    for (ib1 = 1; ib1 <= 3; ++ib1) {
				for (ib2 = 1; ib2 <= 3; ++ib2) {
				    b[0] = val[ib1 - 1];
				    b[2] = val[ib2 - 1] * -2.;
				    tr[0] = itval[((itr << 1) + 1 << 1) - 6] *
					     val[itrscl - 1];
				    tr[1] = itval[((itr << 1) + 1 << 1) - 5] *
					     val[itrscl - 1];
				    tr[2] = itval[((itr << 1) + 2 << 1) - 6] *
					     val[itrscl - 1];
				    tr[3] = itval[((itr << 1) + 2 << 1) - 5] *
					     val[itrscl - 1];
				    tl[0] = val[itl - 1];
				    ++(*knt);
				    dlasy2_(&ltranl, &ltranr, &isgn, &n1, &n2,
					     tl, &c__2, tr, &c__2, b, &c__2, &
					    scale, x, &c__2, &xnorm, &info);
				    if (info != 0) {
					++(*ninfo);
				    }
				    if (ltranr) {
					tmp = tr[2];
					tr[2] = tr[1];
					tr[1] = tmp;
				    } 				    tnrm = abs(tl[0]) + abs(tr[0]) + abs(tr[2] 					    ) + abs(tr[3]) + abs(tr[1]);
				    xnrm = abs(x[0]) + abs(x[2]);
				    res = (d__1 = (tl[0] + sgn * tr[0]) * x[0]
					     + sgn * tr[1] * x[2] - scale * b[
					    0], abs(d__1));
				    res += (d__1 = (tl[0] + sgn * tr[3]) * x[
					    2] + sgn * tr[2] * x[0] - scale * 
					    b[2], abs(d__1));
/* Computing MAX */
				    d__1 = smlnum, d__2 = smlnum * xnrm, d__1 
					    = max(d__1,d__2), d__2 = tnrm * 
					    eps * xnrm;
				    den = max(d__1,d__2);
				    res /= den;
				    if (scale > 1.) {
					res += 1. / eps;
				    }
				    res += (d__1 = xnorm - xnrm, abs(d__1)) / 
					    max(smlnum,xnorm) / eps;
				    if (res > *rmax) {
					*lmax = *knt;
					*rmax = res;
				    }
/* L90: */
				}
/* L100: */
			    }
/* L110: */
			}
/* L120: */
		    }
/* L130: */
		}

		n1 = 2;
		n2 = 2;
		for (itr = 1; itr <= 8; ++itr) {
		    for (itrscl = 1; itrscl <= 3; ++itrscl) {
			for (itl = 1; itl <= 8; ++itl) {
			    for (itlscl = 1; itlscl <= 3; ++itlscl) {
				for (ib1 = 1; ib1 <= 3; ++ib1) {
				    for (ib2 = 1; ib2 <= 3; ++ib2) {
					for (ib3 = 1; ib3 <= 3; ++ib3) {
					    b[0] = val[ib1 - 1];
					    b[1] = val[ib2 - 1] * -4.;
					    b[2] = val[ib3 - 1] * -2.;
/* Computing MIN */
					    d__1 = val[ib1 - 1], d__2 = val[
						    ib2 - 1], d__1 = min(d__1,
						    d__2), d__2 = val[ib3 - 1]
						    ;
					    b[3] = min(d__1,d__2) * 8.;
					    tr[0] = itval[((itr << 1) + 1 << 
						    1) - 6] * val[itrscl - 1];
					    tr[1] = itval[((itr << 1) + 1 << 
						    1) - 5] * val[itrscl - 1];
					    tr[2] = itval[((itr << 1) + 2 << 
						    1) - 6] * val[itrscl - 1];
					    tr[3] = itval[((itr << 1) + 2 << 
						    1) - 5] * val[itrscl - 1];
					    tl[0] = itval[((itl << 1) + 1 << 
						    1) - 6] * val[itlscl - 1];
					    tl[1] = itval[((itl << 1) + 1 << 
						    1) - 5] * val[itlscl - 1];
					    tl[2] = itval[((itl << 1) + 2 << 
						    1) - 6] * val[itlscl - 1];
					    tl[3] = itval[((itl << 1) + 2 << 
						    1) - 5] * val[itlscl - 1];
					    ++(*knt);
					    dlasy2_(&ltranl, &ltranr, &isgn, &
						    n1, &n2, tl, &c__2, tr, &
						    c__2, b, &c__2, &scale, x,
						     &c__2, &xnorm, &info);
					    if (info != 0) {
			  ++(*ninfo);
					    }
					    if (ltranr) {
			  tmp = tr[2];
			  tr[2] = tr[1];
			  tr[1] = tmp;
					    }
					    if (ltranl) {
			  tmp = tl[2];
			  tl[2] = tl[1];
			  tl[1] = tmp;
					    }
					    tnrm = abs(tr[0]) + abs(tr[1]) + 
						    abs(tr[2]) + abs(tr[3]) + 
						    abs(tl[0]) + abs(tl[1]) + 
						    abs(tl[2]) + abs(tl[3]);
/* Computing MAX */
					    d__1 = abs(x[0]) + abs(x[2]), 
						    d__2 = abs(x[1]) + abs(x[
						    3]);
					    xnrm = max(d__1,d__2); 					    res = (d__1 = (tl[0] + sgn * tr[0] 						    ) * x[0] + sgn * tr[1] * 
						    x[2] + tl[2] * x[1] - 
						    scale * b[0], abs(d__1));
					    res += (d__1 = tl[0] * x[2] + sgn 
						    * tr[2] * x[0] + sgn * tr[
						    3] * x[2] + tl[2] * x[3] 
						    - scale * b[2], abs(d__1))
						    ;
					    res += (d__1 = tl[1] * x[0] + sgn 
						    * tr[0] * x[1] + sgn * tr[
						    1] * x[3] + tl[3] * x[1] 
						    - scale * b[1], abs(d__1))
						    ;
					    res += (d__1 = (tl[3] + sgn * tr[
						    3]) * x[3] + sgn * tr[2] *
						     x[1] + tl[1] * x[2] - 
						    scale * b[3], abs(d__1));
/* Computing MAX */
					    d__1 = smlnum, d__2 = smlnum * 
						    xnrm, d__1 = max(d__1,
						    d__2), d__2 = tnrm * eps *
						     xnrm;
					    den = max(d__1,d__2);
					    res /= den;
					    if (scale > 1.) {
			  res += 1. / eps;
					    }
					    res += (d__1 = xnorm - xnrm, abs(
						    d__1)) / max(smlnum,xnorm)
						     / eps;
					    if (res > *rmax) {
			  *lmax = *knt;
			  *rmax = res;
					    }
/* L140: */
					}
/* L150: */
				    }
/* L160: */
				}
/* L170: */
			    }
/* L180: */
			}
/* L190: */
		    }
/* L200: */
		}
/* L210: */
	    }
/* L220: */
	}
/* L230: */
    }

    return 0;

/*     End of DGET32 */

} /* dget32_ */