sget39.c

Go to the documentation of this file.

/* 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"
#include "blaswrap.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 */
    real b[10], d__[20];
    integer i__, j, k, n;
    real t[100] /* was [10][10] */, w, x[20], y[20], vm1[5], vm2[5], vm3[5], 
            vm4[5], vm5[3], dum[1], eps;
    integer ivm1, ivm2, ivm3, ivm4, ivm5, ndim, info;
    real dumm;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    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 *);
    real xnorm;
    extern /* Subroutine */ int slabad_(real *, real *);
    extern doublereal slamch_(char *), slange_(char *, integer *, 
            integer *, real *, integer *, real *);
    real bignum;
    extern integer isamax_(integer *, real *, integer *);
    real normtb;
    extern /* Subroutine */ int slaqtr_(logical *, logical *, integer *, real 
            *, integer *, real *, real *, real *, real *, real *, integer *);
    real smlnum;


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

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

/*  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. */

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

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

/*     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_ */