#include "blaswrap.h"
/*  -- translated by f2c (version 19990503).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#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 integer ndim, info;
    static real dumm;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    static real norm, work[10], b[10], d__[20];
    static integer i__, j, k, n;
    static real t[100]	/* was [10][10] */, scale, w, x[20], y[20], 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, vm1[5], vm2[5], vm3[5], vm4[5], vm5[3], dum[1], eps;
    static integer ivm1, ivm2, ivm3, ivm4, ivm5;


#define ival_ref(a_1,a_2,a_3) ival[((a_3)*5 + (a_2))*5 + a_1 - 31]
#define t_ref(a_1,a_2) t[(a_2)*10 + a_1 - 11]


/*  -- LAPACK test routine (version 3.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       September 30, 1994   


    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_ref(i__, j) = (real) ival_ref(i__, j, 
					    ndim) * vm1[ivm1 - 1];
				    if (i__ >= j) {
					t_ref(i__, j) = t_ref(i__, j) * 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_ */

#undef t_ref
#undef ival_ref


