#include "blaswrap.h"
/* sget39.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__10 = 10;
static integer c__1 = 1;
static logical c_false = FALSE_;
static logical c_true = TRUE_;
static real c_b25 = 1.f;
static real c_b59 = -1.f;

/* Subroutine */ int sget39_(real *rmax, integer *lmax, integer *ninfo, 
	integer *knt)
{
    /* Initialized data */

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

    /* System generated locals */
    integer i__1, i__2;
    real r__1, r__2;

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

    /* Local variables */
    static real b[10], d__[20];
    static integer i__, j, k, n;
    static real t[100]	/* was [10][10] */, w, x[20], y[20], vm1[5], vm2[5], 
	    vm3[5], vm4[5], vm5[3], dum[1], eps;
    static integer ivm1, ivm2, ivm3, ivm4, ivm5, ndim, info;
    static real dumm;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    static real norm, work[10], scale, domin, resid;
    extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 
	    real *, integer *, real *, integer *, real *, real *, integer *);
    extern doublereal sasum_(integer *, real *, integer *);
    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
	    integer *);
    static real xnorm;
    extern /* Subroutine */ int slabad_(real *, real *);
    extern doublereal slamch_(char *), slange_(char *, integer *, 
	    integer *, real *, integer *, real *);
    static real bignum;
    extern integer isamax_(integer *, real *, integer *);
    static real normtb;
    extern /* Subroutine */ int slaqtr_(logical *, logical *, integer *, real 
	    *, integer *, real *, real *, real *, real *, real *, integer *);
    static real smlnum;


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


    Purpose   
    =======   

    SGET39 tests SLAQTR, a routine for solving the real or   
    special complex quasi upper triangular system   

         op(T)*p = scale*c,   
    or   
         op(T + iB)*(p+iq) = scale*(c+id),   

    in real arithmetic. T is upper quasi-triangular.   
    If it is complex, then the first diagonal block of T must be   
    1 by 1, B has the special structure   

                   B = [ b(1) b(2) ... b(n) ]   
                       [       w            ]   
                       [           w        ]   
                       [              .     ]   
                       [                 w  ]   

    op(A) = A or A', where A' denotes the conjugate transpose of   
    the matrix A.   

    On input, X = [ c ].  On output, X = [ p ].   
                  [ d ]                  [ q ]   

    Scale is an output less than or equal to 1, chosen to avoid   
    overflow in X.   
    This subroutine is specially designed for the condition number   
    estimation in the eigenproblem routine STRSNA.   

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

         ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||   
       -----------------------------------------   
           max(ulp*(||T||+||B||)*(||x1||+||x2||),   
               (||T||+||B||)*smlnum/ulp,   
               smlnum)   

    (The (||T||+||B||)*smlnum/ulp term accounts for possible   
     (gradual or nongradual) underflow in x1 and x2.)   

    Arguments   
    ==========   

    RMAX    (output) REAL   
            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 = slamch_("P");
    smlnum = slamch_("S");
    bignum = 1.f / smlnum;
    slabad_(&smlnum, &bignum);

/*     Set up test case parameters */

    vm1[0] = 1.f;
    vm1[1] = sqrt(smlnum);
    vm1[2] = sqrt(vm1[1]);
    vm1[3] = sqrt(bignum);
    vm1[4] = sqrt(vm1[3]);

    vm2[0] = 1.f;
    vm2[1] = sqrt(smlnum);
    vm2[2] = sqrt(vm2[1]);
    vm2[3] = sqrt(bignum);
    vm2[4] = sqrt(vm2[3]);

    vm3[0] = 1.f;
    vm3[1] = sqrt(smlnum);
    vm3[2] = sqrt(vm3[1]);
    vm3[3] = sqrt(bignum);
    vm3[4] = sqrt(vm3[3]);

    vm4[0] = 1.f;
    vm4[1] = sqrt(smlnum);
    vm4[2] = sqrt(vm4[1]);
    vm4[3] = sqrt(bignum);
    vm4[4] = sqrt(vm4[3]);

    vm5[0] = 1.f;
    vm5[1] = eps;
    vm5[2] = sqrt(smlnum);

/*     Initalization */

    *knt = 0;
    *rmax = 0.f;
    *ninfo = 0;
    smlnum /= eps;

/*     Begin test loop */

    for (ivm5 = 1; ivm5 <= 3; ++ivm5) {
	for (ivm4 = 1; ivm4 <= 5; ++ivm4) {
	    for (ivm3 = 1; ivm3 <= 5; ++ivm3) {
		for (ivm2 = 1; ivm2 <= 5; ++ivm2) {
		    for (ivm1 = 1; ivm1 <= 5; ++ivm1) {
			for (ndim = 1; ndim <= 6; ++ndim) {

			    n = idim[ndim - 1];
			    i__1 = n;
			    for (i__ = 1; i__ <= i__1; ++i__) {
				i__2 = n;
				for (j = 1; j <= i__2; ++j) {
				    t[i__ + j * 10 - 11] = (real) ival[i__ + (
					    j + ndim * 5) * 5 - 31] * vm1[
					    ivm1 - 1];
				    if (i__ >= j) {
					t[i__ + j * 10 - 11] *= vm5[ivm5 - 1];
				    }
/* L10: */
				}
/* L20: */
			    }

			    w = vm2[ivm2 - 1] * 1.f;

			    i__1 = n;
			    for (i__ = 1; i__ <= i__1; ++i__) {
				b[i__ - 1] = cos((real) i__) * vm3[ivm3 - 1];
/* L30: */
			    }

			    i__1 = n << 1;
			    for (i__ = 1; i__ <= i__1; ++i__) {
				d__[i__ - 1] = sin((real) i__) * vm4[ivm4 - 1]
					;
/* L40: */
			    }

			    norm = slange_("1", &n, &n, t, &c__10, work);
			    k = isamax_(&n, b, &c__1);
			    normtb = norm + (r__1 = b[k - 1], dabs(r__1)) + 
				    dabs(w);

			    scopy_(&n, d__, &c__1, x, &c__1);
			    ++(*knt);
			    slaqtr_(&c_false, &c_true, &n, t, &c__10, dum, &
				    dumm, &scale, x, work, &info);
			    if (info != 0) {
				++(*ninfo);
			    }

/*                       || T*x - scale*d || /   
                           max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) */

			    scopy_(&n, d__, &c__1, y, &c__1);
			    r__1 = -scale;
			    sgemv_("No transpose", &n, &n, &c_b25, t, &c__10, 
				    x, &c__1, &r__1, y, &c__1);
			    xnorm = sasum_(&n, x, &c__1);
			    resid = sasum_(&n, y, &c__1);
/* Computing MAX */
			    r__1 = smlnum, r__2 = smlnum / eps * norm, r__1 = 
				    max(r__1,r__2), r__2 = norm * eps * xnorm;
			    domin = dmax(r__1,r__2);
			    resid /= domin;
			    if (resid > *rmax) {
				*rmax = resid;
				*lmax = *knt;
			    }

			    scopy_(&n, d__, &c__1, x, &c__1);
			    ++(*knt);
			    slaqtr_(&c_true, &c_true, &n, t, &c__10, dum, &
				    dumm, &scale, x, work, &info);
			    if (info != 0) {
				++(*ninfo);
			    }

/*                       || T*x - scale*d || /   
                           max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum) */

			    scopy_(&n, d__, &c__1, y, &c__1);
			    r__1 = -scale;
			    sgemv_("Transpose", &n, &n, &c_b25, t, &c__10, x, 
				    &c__1, &r__1, y, &c__1);
			    xnorm = sasum_(&n, x, &c__1);
			    resid = sasum_(&n, y, &c__1);
/* Computing MAX */
			    r__1 = smlnum, r__2 = smlnum / eps * norm, r__1 = 
				    max(r__1,r__2), r__2 = norm * eps * xnorm;
			    domin = dmax(r__1,r__2);
			    resid /= domin;
			    if (resid > *rmax) {
				*rmax = resid;
				*lmax = *knt;
			    }

			    i__1 = n << 1;
			    scopy_(&i__1, d__, &c__1, x, &c__1);
			    ++(*knt);
			    slaqtr_(&c_false, &c_false, &n, t, &c__10, b, &w, 
				    &scale, x, work, &info);
			    if (info != 0) {
				++(*ninfo);
			    }

/*                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /   
                            max(ulp*(||T||+||B||)*(||x1||+||x2||),   
                                    smlnum/ulp * (||T||+||B||), smlnum ) */


			    i__1 = n << 1;
			    scopy_(&i__1, d__, &c__1, y, &c__1);
			    y[0] = sdot_(&n, b, &c__1, &x[n], &c__1) + scale *
				     y[0];
			    i__1 = n;
			    for (i__ = 2; i__ <= i__1; ++i__) {
				y[i__ - 1] = w * x[i__ + n - 1] + scale * y[
					i__ - 1];
/* L50: */
			    }
			    sgemv_("No transpose", &n, &n, &c_b25, t, &c__10, 
				    x, &c__1, &c_b59, y, &c__1);

			    y[n] = sdot_(&n, b, &c__1, x, &c__1) - scale * y[
				    n];
			    i__1 = n;
			    for (i__ = 2; i__ <= i__1; ++i__) {
				y[i__ + n - 1] = w * x[i__ - 1] - scale * y[
					i__ + n - 1];
/* L60: */
			    }
			    sgemv_("No transpose", &n, &n, &c_b25, t, &c__10, 
				    &x[n], &c__1, &c_b25, &y[n], &c__1);

			    i__1 = n << 1;
			    resid = sasum_(&i__1, y, &c__1);
/* Computing MAX */
			    i__1 = n << 1;
			    r__1 = smlnum, r__2 = smlnum / eps * normtb, r__1 
				    = max(r__1,r__2), r__2 = eps * (normtb * 
				    sasum_(&i__1, x, &c__1));
			    domin = dmax(r__1,r__2);
			    resid /= domin;
			    if (resid > *rmax) {
				*rmax = resid;
				*lmax = *knt;
			    }

			    i__1 = n << 1;
			    scopy_(&i__1, d__, &c__1, x, &c__1);
			    ++(*knt);
			    slaqtr_(&c_true, &c_false, &n, t, &c__10, b, &w, &
				    scale, x, work, &info);
			    if (info != 0) {
				++(*ninfo);
			    }

/*                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /   
                            max(ulp*(||T||+||B||)*(||x1||+||x2||),   
                                    smlnum/ulp * (||T||+||B||), smlnum ) */

			    i__1 = n << 1;
			    scopy_(&i__1, d__, &c__1, y, &c__1);
			    y[0] = b[0] * x[n] - scale * y[0];
			    i__1 = n;
			    for (i__ = 2; i__ <= i__1; ++i__) {
				y[i__ - 1] = b[i__ - 1] * x[n] + w * x[i__ + 
					n - 1] - scale * y[i__ - 1];
/* L70: */
			    }
			    sgemv_("Transpose", &n, &n, &c_b25, t, &c__10, x, 
				    &c__1, &c_b25, y, &c__1);

			    y[n] = b[0] * x[0] + scale * y[n];
			    i__1 = n;
			    for (i__ = 2; i__ <= i__1; ++i__) {
				y[i__ + n - 1] = b[i__ - 1] * x[0] + w * x[
					i__ - 1] + scale * y[i__ + n - 1];
/* L80: */
			    }
			    sgemv_("Transpose", &n, &n, &c_b25, t, &c__10, &x[
				    n], &c__1, &c_b59, &y[n], &c__1);

			    i__1 = n << 1;
			    resid = sasum_(&i__1, y, &c__1);
/* Computing MAX */
			    i__1 = n << 1;
			    r__1 = smlnum, r__2 = smlnum / eps * normtb, r__1 
				    = max(r__1,r__2), r__2 = eps * (normtb * 
				    sasum_(&i__1, x, &c__1));
			    domin = dmax(r__1,r__2);
			    resid /= domin;
			    if (resid > *rmax) {
				*rmax = resid;
				*lmax = *knt;
			    }

/* L90: */
			}
/* L100: */
		    }
/* L110: */
		}
/* L120: */
	    }
/* L130: */
	}
/* L140: */
    }

    return 0;

/*     End of SGET39 */

} /* sget39_ */