#include "blaswrap.h" /* ssvdct.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" /* Subroutine */ int ssvdct_(integer *n, real *s, real *e, real *shift, integer *num) { /* System generated locals */ integer i__1; real r__1, r__2, r__3, r__4; /* Builtin functions */ double sqrt(doublereal); /* Local variables */ static integer i__; static real u, m1, m2, mx, tmp, tom, sun, sov, unfl, ovfl, ssun; extern doublereal slamch_(char *); static real sshift; /* -- LAPACK test routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= SSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N tridiagonal matrix T which are less than or equal to SHIFT. T is formed by putting zeros on the diagonal and making the off-diagonals equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N). If SHIFT is positive, NUM is equal to N plus the number of singular values of a bidiagonal matrix B less than or equal to SHIFT. Here B has diagonal entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1). If SHIFT is negative, NUM is equal to the number of singular values of B greater than or equal to -SHIFT. See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966 Arguments ========= N (input) INTEGER The dimension of the bidiagonal matrix B. S (input) REAL array, dimension (N) The diagonal entries of the bidiagonal matrix B. E (input) REAL array of dimension (N-1) The superdiagonal entries of the bidiagonal matrix B. SHIFT (input) REAL The shift, used as described under Purpose. NUM (output) INTEGER The number of eigenvalues of T less than or equal to SHIFT. ===================================================================== Get machine constants Parameter adjustments */ --e; --s; /* Function Body */ unfl = slamch_("Safe minimum") * 2; ovfl = 1.f / unfl; /* Find largest entry */ mx = dabs(s[1]); i__1 = *n - 1; for (i__ = 1; i__ <= i__1; ++i__) { /* Computing MAX */ r__3 = mx, r__4 = (r__1 = s[i__ + 1], dabs(r__1)), r__3 = max(r__3, r__4), r__4 = (r__2 = e[i__], dabs(r__2)); mx = dmax(r__3,r__4); /* L10: */ } if (mx == 0.f) { if (*shift < 0.f) { *num = 0; } else { *num = *n << 1; } return 0; } /* Compute scale factors as in Kahan's report */ sun = sqrt(unfl); ssun = sqrt(sun); sov = sqrt(ovfl); tom = ssun * sov; if (mx <= 1.f) { m1 = 1.f / mx; m2 = tom; } else { m1 = 1.f; m2 = tom / mx; } /* Begin counting */ u = 1.f; *num = 0; sshift = *shift * m1 * m2; u = -sshift; if (u <= sun) { if (u <= 0.f) { ++(*num); if (u > -sun) { u = -sun; } } else { u = sun; } } tmp = s[1] * m1 * m2; u = -tmp * (tmp / u) - sshift; if (u <= sun) { if (u <= 0.f) { ++(*num); if (u > -sun) { u = -sun; } } else { u = sun; } } i__1 = *n - 1; for (i__ = 1; i__ <= i__1; ++i__) { tmp = e[i__] * m1 * m2; u = -tmp * (tmp / u) - sshift; if (u <= sun) { if (u <= 0.f) { ++(*num); if (u > -sun) { u = -sun; } } else { u = sun; } } tmp = s[i__ + 1] * m1 * m2; u = -tmp * (tmp / u) - sshift; if (u <= sun) { if (u <= 0.f) { ++(*num); if (u > -sun) { u = -sun; } } else { u = sun; } } /* L20: */ } return 0; /* End of SSVDCT */ } /* ssvdct_ */