sget33.c

Go to the documentation of this file.

/* sget33.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 real c_b19 = 1.f;

/* Subroutine */ int sget33_(real *rmax, integer *lmax, integer *ninfo, 
        integer *knt)
{
    /* System generated locals */
    real r__1, r__2, r__3;

    /* Builtin functions */
    double r_sign(real *, real *);

    /* Local variables */
    real q[4]   /* was [2][2] */, t[4]  /* was [2][2] */;
    integer i1, i2, i3, i4, j1, j2, j3;
    real t1[4]  /* was [2][2] */, t2[4] /* was [2][2] */, cs, sn, vm[3];
    integer im1, im2, im3, im4;
    real wi1, wi2, wr1, wr2, val[4], eps, res, sum, tnrm;
    extern /* Subroutine */ int slanv2_(real *, real *, real *, real *, real *
, real *, real *, real *, real *, real *), slabad_(real *, real *)
            ;
    extern doublereal slamch_(char *);
    real bignum, smlnum;


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

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

/*  Purpose */
/*  ======= */

/*  SGET33 tests SLANV2, a routine for putting 2 by 2 blocks into */
/*  standard form.  In other words, it computes a two by two rotation */
/*  [[C,S];[-S,C]] where in */

/*     [ C S ][T(1,1) T(1,2)][ C -S ] = [ T11 T12 ] */
/*     [-S C ][T(2,1) T(2,2)][ S  C ]   [ T21 T22 ] */

/*  either */
/*     1) T21=0 (real eigenvalues), or */
/*     2) T11=T22 and T21*T12<0 (complex conjugate eigenvalues). */
/*  We also  verify that the residual is small. */

/*  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 returned with INFO .NE. 0. */

/*  KNT     (output) INTEGER */
/*          Total number of examples tested. */

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

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

/*     Get machine parameters */

    eps = slamch_("P");
    smlnum = slamch_("S") / eps;
    bignum = 1.f / smlnum;
    slabad_(&smlnum, &bignum);

/*     Set up test case parameters */

    val[0] = 1.f;
    val[1] = eps * 2.f + 1.f;
    val[2] = 2.f;
    val[3] = 2.f - eps * 4.f;
    vm[0] = smlnum;
    vm[1] = 1.f;
    vm[2] = bignum;

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

/*     Begin test loop */

    for (i1 = 1; i1 <= 4; ++i1) {
        for (i2 = 1; i2 <= 4; ++i2) {
            for (i3 = 1; i3 <= 4; ++i3) {
                for (i4 = 1; i4 <= 4; ++i4) {
                    for (im1 = 1; im1 <= 3; ++im1) {
                        for (im2 = 1; im2 <= 3; ++im2) {
                            for (im3 = 1; im3 <= 3; ++im3) {
                                for (im4 = 1; im4 <= 3; ++im4) {
                                    t[0] = val[i1 - 1] * vm[im1 - 1];
                                    t[2] = val[i2 - 1] * vm[im2 - 1];
                                    t[1] = -val[i3 - 1] * vm[im3 - 1];
                                    t[3] = val[i4 - 1] * vm[im4 - 1];
/* Computing MAX */
                                    r__1 = dabs(t[0]), r__2 = dabs(t[2]), 
                                            r__1 = max(r__1,r__2), r__2 = 
                                            dabs(t[1]), r__1 = max(r__1,r__2),
                                             r__2 = dabs(t[3]);
                                    tnrm = dmax(r__1,r__2);
                                    t1[0] = t[0];
                                    t1[2] = t[2];
                                    t1[1] = t[1];
                                    t1[3] = t[3];
                                    q[0] = 1.f;
                                    q[2] = 0.f;
                                    q[1] = 0.f;
                                    q[3] = 1.f;

                                    slanv2_(t, &t[2], &t[1], &t[3], &wr1, &
                                            wi1, &wr2, &wi2, &cs, &sn);
                                    for (j1 = 1; j1 <= 2; ++j1) {
                                        res = q[j1 - 1] * cs + q[j1 + 1] * sn;
                                        q[j1 + 1] = -q[j1 - 1] * sn + q[j1 + 
                                                1] * cs;
                                        q[j1 - 1] = res;
/* L10: */
                                    }

                                    res = 0.f;
/* Computing 2nd power */
                                    r__2 = q[0];
/* Computing 2nd power */
                                    r__3 = q[2];
                                    res += (r__1 = r__2 * r__2 + r__3 * r__3 
                                            - 1.f, dabs(r__1)) / eps;
/* Computing 2nd power */
                                    r__2 = q[3];
/* Computing 2nd power */
                                    r__3 = q[1];
                                    res += (r__1 = r__2 * r__2 + r__3 * r__3 
                                            - 1.f, dabs(r__1)) / eps;
                                    res += (r__1 = q[0] * q[1] + q[2] * q[3], 
                                            dabs(r__1)) / eps;
                                    for (j1 = 1; j1 <= 2; ++j1) {
                                        for (j2 = 1; j2 <= 2; ++j2) {
                                            t2[j1 + (j2 << 1) - 3] = 0.f;
                                            for (j3 = 1; j3 <= 2; ++j3) {
                          t2[j1 + (j2 << 1) - 3] += t1[j1 + (j3 << 1) - 3] * 
                                  q[j3 + (j2 << 1) - 3];
/* L20: */
                                            }
/* L30: */
                                        }
/* L40: */
                                    }
                                    for (j1 = 1; j1 <= 2; ++j1) {
                                        for (j2 = 1; j2 <= 2; ++j2) {
                                            sum = t[j1 + (j2 << 1) - 3];
                                            for (j3 = 1; j3 <= 2; ++j3) {
                          sum -= q[j3 + (j1 << 1) - 3] * t2[j3 + (j2 << 1) - 
                                  3];
/* L50: */
                                            }
                                            res += dabs(sum) / eps / tnrm;
/* L60: */
                                        }
/* L70: */
                                    }
                                    if (t[1] != 0.f && (t[0] != t[3] || 
                                            r_sign(&c_b19, &t[2]) * r_sign(&
                                            c_b19, &t[1]) > 0.f)) {
                                        res += 1.f / eps;
                                    }
                                    ++(*knt);
                                    if (res > *rmax) {
                                        *lmax = *knt;
                                        *rmax = res;
                                    }
/* L80: */
                                }
/* L90: */
                            }
/* L100: */
                        }
/* L110: */
                    }
/* L120: */
                }
/* L130: */
            }
/* L140: */
        }
/* L150: */
    }

    return 0;

/*     End of SGET33 */

} /* sget33_ */