dget32.c

Go to the documentation of this file.

/* 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"
#include "blaswrap.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 */
    doublereal b[4]     /* was [2][2] */, x[4]  /* was [2][2] */;
    integer n1, n2, ib;
    doublereal tl[4]    /* was [2][2] */, tr[4] /* was [2][2] */;
    integer ib1, ib2, ib3;
    doublereal den, val[3], eps;
    integer itl;
    doublereal res, sgn;
    integer itr;
    doublereal tmp;
    integer info, isgn;
    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 *);
    doublereal bignum;
    integer itranl, itlscl;
    logical ltranl;
    integer itranr, itrscl;
    logical ltranr;
    doublereal smlnum;


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

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

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

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

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

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