#include "blaswrap.h" /* sget34.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__16 = 16; static integer c__0 = 0; static integer c__1 = 1; static integer c__4 = 4; static integer c__5 = 5; static logical c_true = TRUE_; static integer c__2 = 2; static integer c__32 = 32; static integer c__3 = 3; static real c_b64 = 1.f; /* Subroutine */ int sget34_(real *rmax, integer *lmax, integer *ninfo, integer *knt) { /* System generated locals */ real r__1, r__2, r__3; /* Builtin functions */ double sqrt(doublereal), r_sign(real *, real *); /* Local variables */ static integer i__, j; static real q[16] /* was [4][4] */, t[16] /* was [4][4] */, t1[16] /* was [4][4] */; static integer ia, ib, ic; static real vm[2]; static integer ia11, ia12, ia21, ia22, ic11, ic12, ic21, ic22, iam, icm; static real val[9], eps, res; static integer info; static real tnrm, work[32]; extern /* Subroutine */ int shst01_(integer *, integer *, integer *, real *, integer *, real *, integer *, real *, integer *, real *, integer *, real *), scopy_(integer *, real *, integer *, real *, integer *), slabad_(real *, real *); extern doublereal slamch_(char *); extern /* Subroutine */ int slaexc_(logical *, integer *, real *, integer *, real *, integer *, integer *, integer *, integer *, real *, integer *); static real bignum, smlnum, result[2]; /* -- LAPACK test routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= SGET34 tests SLAEXC, a routine for swapping adjacent blocks (either 1 by 1 or 2 by 2) on the diagonal of a matrix in real Schur form. Thus, SLAEXC computes an orthogonal matrix Q such that Q' * [ A B ] * Q = [ C1 B1 ] [ 0 C ] [ 0 A1 ] where C1 is similar to C and A1 is similar to A. Both A and C are assumed to be in standard form (equal diagonal entries and offdiagonal with differing signs) and A1 and C1 are returned with the same properties. The test code verifies these last last assertions, as well as that the residual in the above equation 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 array, dimension (2) NINFO(J) is the number of examples where INFO=J occurred. KNT (output) INTEGER Total number of examples tested. ===================================================================== Get machine parameters Parameter adjustments */ --ninfo; /* Function Body */ eps = slamch_("P"); smlnum = slamch_("S") / eps; bignum = 1.f / smlnum; slabad_(&smlnum, &bignum); /* Set up test case parameters */ val[0] = 0.f; val[1] = sqrt(smlnum); val[2] = 1.f; val[3] = 2.f; val[4] = sqrt(bignum); val[5] = -sqrt(smlnum); val[6] = -1.f; val[7] = -2.f; val[8] = -sqrt(bignum); vm[0] = 1.f; vm[1] = eps * 2.f + 1.f; scopy_(&c__16, &val[3], &c__0, t, &c__1); ninfo[1] = 0; ninfo[2] = 0; *knt = 0; *lmax = 0; *rmax = 0.f; /* Begin test loop */ for (ia = 1; ia <= 9; ++ia) { for (iam = 1; iam <= 2; ++iam) { for (ib = 1; ib <= 9; ++ib) { for (ic = 1; ic <= 9; ++ic) { t[0] = val[ia - 1] * vm[iam - 1]; t[5] = val[ic - 1]; t[4] = val[ib - 1]; t[1] = 0.f; /* Computing MAX */ r__1 = dabs(t[0]), r__2 = dabs(t[5]), r__1 = max(r__1, r__2), r__2 = dabs(t[4]); tnrm = dmax(r__1,r__2); scopy_(&c__16, t, &c__1, t1, &c__1); scopy_(&c__16, val, &c__0, q, &c__1); scopy_(&c__4, &val[2], &c__0, q, &c__5); slaexc_(&c_true, &c__2, t, &c__4, q, &c__4, &c__1, &c__1, &c__1, work, &info); if (info != 0) { ++ninfo[info]; } shst01_(&c__2, &c__1, &c__2, t1, &c__4, t, &c__4, q, & c__4, work, &c__32, result); res = result[0] + result[1]; if (info != 0) { res += 1.f / eps; } if (t[0] != t1[5]) { res += 1.f / eps; } if (t[5] != t1[0]) { res += 1.f / eps; } if (t[1] != 0.f) { res += 1.f / eps; } ++(*knt); if (res > *rmax) { *lmax = *knt; *rmax = res; } /* L10: */ } /* L20: */ } /* L30: */ } /* L40: */ } for (ia = 1; ia <= 5; ++ia) { for (iam = 1; iam <= 2; ++iam) { for (ib = 1; ib <= 5; ++ib) { for (ic11 = 1; ic11 <= 5; ++ic11) { for (ic12 = 2; ic12 <= 5; ++ic12) { for (ic21 = 2; ic21 <= 4; ++ic21) { for (ic22 = -1; ic22 <= 1; ic22 += 2) { t[0] = val[ia - 1] * vm[iam - 1]; t[4] = val[ib - 1]; t[8] = val[ib - 1] * -2.f; t[1] = 0.f; t[5] = val[ic11 - 1]; t[9] = val[ic12 - 1]; t[2] = 0.f; t[6] = -val[ic21 - 1]; t[10] = val[ic11 - 1] * (real) ic22; /* Computing MAX */ r__1 = dabs(t[0]), r__2 = dabs(t[4]), r__1 = max(r__1,r__2), r__2 = dabs(t[8]), r__1 = max(r__1,r__2), r__2 = dabs(t[ 5]), r__1 = max(r__1,r__2), r__2 = dabs(t[9]), r__1 = max(r__1,r__2), r__2 = dabs(t[6]), r__1 = max(r__1, r__2), r__2 = dabs(t[10]); tnrm = dmax(r__1,r__2); scopy_(&c__16, t, &c__1, t1, &c__1); scopy_(&c__16, val, &c__0, q, &c__1); scopy_(&c__4, &val[2], &c__0, q, &c__5); slaexc_(&c_true, &c__3, t, &c__4, q, &c__4, & c__1, &c__1, &c__2, work, &info); if (info != 0) { ++ninfo[info]; } shst01_(&c__3, &c__1, &c__3, t1, &c__4, t, & c__4, q, &c__4, work, &c__32, result); res = result[0] + result[1]; if (info == 0) { if (t1[0] != t[10]) { res += 1.f / eps; } if (t[2] != 0.f) { res += 1.f / eps; } if (t[6] != 0.f) { res += 1.f / eps; } if (t[1] != 0.f && (t[0] != t[5] || r_sign(&c_b64, &t[4]) == r_sign(& c_b64, &t[1]))) { res += 1.f / eps; } } ++(*knt); if (res > *rmax) { *lmax = *knt; *rmax = res; } /* L50: */ } /* L60: */ } /* L70: */ } /* L80: */ } /* L90: */ } /* L100: */ } /* L110: */ } for (ia11 = 1; ia11 <= 5; ++ia11) { for (ia12 = 2; ia12 <= 5; ++ia12) { for (ia21 = 2; ia21 <= 4; ++ia21) { for (ia22 = -1; ia22 <= 1; ia22 += 2) { for (icm = 1; icm <= 2; ++icm) { for (ib = 1; ib <= 5; ++ib) { for (ic = 1; ic <= 5; ++ic) { t[0] = val[ia11 - 1]; t[4] = val[ia12 - 1]; t[8] = val[ib - 1] * -2.f; t[1] = -val[ia21 - 1]; t[5] = val[ia11 - 1] * (real) ia22; t[9] = val[ib - 1]; t[2] = 0.f; t[6] = 0.f; t[10] = val[ic - 1] * vm[icm - 1]; /* Computing MAX */ r__1 = dabs(t[0]), r__2 = dabs(t[4]), r__1 = max(r__1,r__2), r__2 = dabs(t[8]), r__1 = max(r__1,r__2), r__2 = dabs(t[ 5]), r__1 = max(r__1,r__2), r__2 = dabs(t[9]), r__1 = max(r__1,r__2), r__2 = dabs(t[6]), r__1 = max(r__1, r__2), r__2 = dabs(t[10]); tnrm = dmax(r__1,r__2); scopy_(&c__16, t, &c__1, t1, &c__1); scopy_(&c__16, val, &c__0, q, &c__1); scopy_(&c__4, &val[2], &c__0, q, &c__5); slaexc_(&c_true, &c__3, t, &c__4, q, &c__4, & c__1, &c__2, &c__1, work, &info); if (info != 0) { ++ninfo[info]; } shst01_(&c__3, &c__1, &c__3, t1, &c__4, t, & c__4, q, &c__4, work, &c__32, result); res = result[0] + result[1]; if (info == 0) { if (t1[10] != t[0]) { res += 1.f / eps; } if (t[1] != 0.f) { res += 1.f / eps; } if (t[2] != 0.f) { res += 1.f / eps; } if (t[6] != 0.f && (t[5] != t[10] || r_sign(&c_b64, &t[9]) == r_sign(& c_b64, &t[6]))) { res += 1.f / eps; } } ++(*knt); if (res > *rmax) { *lmax = *knt; *rmax = res; } /* L120: */ } /* L130: */ } /* L140: */ } /* L150: */ } /* L160: */ } /* L170: */ } /* L180: */ } for (ia11 = 1; ia11 <= 5; ++ia11) { for (ia12 = 2; ia12 <= 5; ++ia12) { for (ia21 = 2; ia21 <= 4; ++ia21) { for (ia22 = -1; ia22 <= 1; ia22 += 2) { for (ib = 1; ib <= 5; ++ib) { for (ic11 = 3; ic11 <= 4; ++ic11) { for (ic12 = 3; ic12 <= 4; ++ic12) { for (ic21 = 3; ic21 <= 4; ++ic21) { for (ic22 = -1; ic22 <= 1; ic22 += 2) { for (icm = 5; icm <= 7; ++icm) { iam = 1; t[0] = val[ia11 - 1] * vm[iam - 1] ; t[4] = val[ia12 - 1] * vm[iam - 1] ; t[8] = val[ib - 1] * -2.f; t[12] = val[ib - 1] * .5f; t[1] = -t[4] * val[ia21 - 1]; t[5] = val[ia11 - 1] * (real) ia22 * vm[iam - 1]; t[9] = val[ib - 1]; t[13] = val[ib - 1] * 3.f; t[2] = 0.f; t[6] = 0.f; t[10] = val[ic11 - 1] * (r__1 = val[icm - 1], dabs(r__1)); t[14] = val[ic12 - 1] * (r__1 = val[icm - 1], dabs(r__1)); t[3] = 0.f; t[7] = 0.f; t[11] = -t[14] * val[ic21 - 1] * ( r__1 = val[icm - 1], dabs( r__1)); t[15] = val[ic11 - 1] * (real) ic22 * (r__1 = val[icm - 1], dabs(r__1)); tnrm = 0.f; for (i__ = 1; i__ <= 4; ++i__) { for (j = 1; j <= 4; ++j) { /* Computing MAX */ r__2 = tnrm, r__3 = (r__1 = t[i__ + (j << 2) - 5], dabs(r__1)); tnrm = dmax(r__2,r__3); /* L190: */ } /* L200: */ } scopy_(&c__16, t, &c__1, t1, & c__1); scopy_(&c__16, val, &c__0, q, & c__1); scopy_(&c__4, &val[2], &c__0, q, & c__5); slaexc_(&c_true, &c__4, t, &c__4, q, &c__4, &c__1, &c__2, & c__2, work, &info); if (info != 0) { ++ninfo[info]; } shst01_(&c__4, &c__1, &c__4, t1, & c__4, t, &c__4, q, &c__4, work, &c__32, result); res = result[0] + result[1]; if (info == 0) { if (t[2] != 0.f) { res += 1.f / eps; } if (t[3] != 0.f) { res += 1.f / eps; } if (t[6] != 0.f) { res += 1.f / eps; } if (t[7] != 0.f) { res += 1.f / eps; } if (t[1] != 0.f && (t[0] != t[5] || r_sign(&c_b64, & t[4]) == r_sign(&c_b64, &t[1]))) { res += 1.f / eps; } if (t[11] != 0.f && (t[10] != t[15] || r_sign(& c_b64, &t[14]) == r_sign(&c_b64, &t[11]))) { res += 1.f / eps; } } ++(*knt); if (res > *rmax) { *lmax = *knt; *rmax = res; } /* L210: */ } /* L220: */ } /* L230: */ } /* L240: */ } /* L250: */ } /* L260: */ } /* L270: */ } /* L280: */ } /* L290: */ } /* L300: */ } return 0; /* End of SGET34 */ } /* sget34_ */