sstect.c
Go to the documentation of this file.
00001 /* sstect.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int sstect_(integer *n, real *a, real *b, real *shift, 
00017         integer *num)
00018 {
00019     /* System generated locals */
00020     integer i__1;
00021     real r__1, r__2, r__3, r__4;
00022 
00023     /* Builtin functions */
00024     double sqrt(doublereal);
00025 
00026     /* Local variables */
00027     integer i__;
00028     real u, m1, m2, mx, tmp, tom, sun, sov, unfl, ovfl, ssun;
00029     extern doublereal slamch_(char *);
00030     real sshift;
00031 
00032 
00033 /*  -- LAPACK test routine (version 3.1) -- */
00034 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00035 /*     November 2006 */
00036 
00037 /*     .. Scalar Arguments .. */
00038 /*     .. */
00039 /*     .. Array Arguments .. */
00040 /*     .. */
00041 
00042 /*  Purpose */
00043 /*  ======= */
00044 
00045 /*     SSTECT counts the number NUM of eigenvalues of a tridiagonal */
00046 /*     matrix T which are less than or equal to SHIFT. T has */
00047 /*     diagonal entries A(1), ... , A(N), and offdiagonal entries */
00048 /*     B(1), ..., B(N-1). */
00049 /*     See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
00050 /*     Matrix", Report CS41, Computer Science Dept., Stanford */
00051 /*     University, July 21, 1966 */
00052 
00053 /*  Arguments */
00054 /*  ========= */
00055 
00056 /*  N       (input) INTEGER */
00057 /*          The dimension of the tridiagonal matrix T. */
00058 
00059 /*  A       (input) REAL array, dimension (N) */
00060 /*          The diagonal entries of the tridiagonal matrix T. */
00061 
00062 /*  B       (input) REAL array, dimension (N-1) */
00063 /*          The offdiagonal entries of the tridiagonal matrix T. */
00064 
00065 /*  SHIFT   (input) REAL */
00066 /*          The shift, used as described under Purpose. */
00067 
00068 /*  NUM     (output) INTEGER */
00069 /*          The number of eigenvalues of T less than or equal */
00070 /*          to SHIFT. */
00071 
00072 /*  ===================================================================== */
00073 
00074 /*     .. Parameters .. */
00075 /*     .. */
00076 /*     .. Local Scalars .. */
00077 /*     .. */
00078 /*     .. External Functions .. */
00079 /*     .. */
00080 /*     .. Intrinsic Functions .. */
00081 /*     .. */
00082 /*     .. Executable Statements .. */
00083 
00084 /*     Get machine constants */
00085 
00086     /* Parameter adjustments */
00087     --b;
00088     --a;
00089 
00090     /* Function Body */
00091     unfl = slamch_("Safe minimum");
00092     ovfl = slamch_("Overflow");
00093 
00094 /*     Find largest entry */
00095 
00096     mx = dabs(a[1]);
00097     i__1 = *n - 1;
00098     for (i__ = 1; i__ <= i__1; ++i__) {
00099 /* Computing MAX */
00100         r__3 = mx, r__4 = (r__1 = a[i__ + 1], dabs(r__1)), r__3 = max(r__3,
00101                 r__4), r__4 = (r__2 = b[i__], dabs(r__2));
00102         mx = dmax(r__3,r__4);
00103 /* L10: */
00104     }
00105 
00106 /*     Handle easy cases, including zero matrix */
00107 
00108     if (*shift >= mx * 3.f) {
00109         *num = *n;
00110         return 0;
00111     }
00112     if (*shift < mx * -3.f) {
00113         *num = 0;
00114         return 0;
00115     }
00116 
00117 /*     Compute scale factors as in Kahan's report */
00118 /*     At this point, MX .NE. 0 so we can divide by it */
00119 
00120     sun = sqrt(unfl);
00121     ssun = sqrt(sun);
00122     sov = sqrt(ovfl);
00123     tom = ssun * sov;
00124     if (mx <= 1.f) {
00125         m1 = 1.f / mx;
00126         m2 = tom;
00127     } else {
00128         m1 = 1.f;
00129         m2 = tom / mx;
00130     }
00131 
00132 /*     Begin counting */
00133 
00134     *num = 0;
00135     sshift = *shift * m1 * m2;
00136     u = a[1] * m1 * m2 - sshift;
00137     if (u <= sun) {
00138         if (u <= 0.f) {
00139             ++(*num);
00140             if (u > -sun) {
00141                 u = -sun;
00142             }
00143         } else {
00144             u = sun;
00145         }
00146     }
00147     i__1 = *n;
00148     for (i__ = 2; i__ <= i__1; ++i__) {
00149         tmp = b[i__ - 1] * m1 * m2;
00150         u = a[i__] * m1 * m2 - tmp * (tmp / u) - sshift;
00151         if (u <= sun) {
00152             if (u <= 0.f) {
00153                 ++(*num);
00154                 if (u > -sun) {
00155                     u = -sun;
00156                 }
00157             } else {
00158                 u = sun;
00159             }
00160         }
00161 /* L20: */
00162     }
00163     return 0;
00164 
00165 /*     End of SSTECT */
00166 
00167 } /* sstect_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:14