ssvdct.c
Go to the documentation of this file.
00001 /* ssvdct.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 ssvdct_(integer *n, real *s, real *e, 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 /*  SSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N */
00046 /*  tridiagonal matrix T which are less than or equal to SHIFT.  T is */
00047 /*  formed by putting zeros on the diagonal and making the off-diagonals */
00048 /*  equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N).  If SHIFT is */
00049 /*  positive, NUM is equal to N plus the number of singular values of a */
00050 /*  bidiagonal matrix B less than or equal to SHIFT.  Here B has diagonal */
00051 /*  entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1). */
00052 /*  If SHIFT is negative, NUM is equal to the number of singular values */
00053 /*  of B greater than or equal to -SHIFT. */
00054 
00055 /*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
00056 /*  Matrix", Report CS41, Computer Science Dept., Stanford University, */
00057 /*  July 21, 1966 */
00058 
00059 /*  Arguments */
00060 /*  ========= */
00061 
00062 /*  N       (input) INTEGER */
00063 /*          The dimension of the bidiagonal matrix B. */
00064 
00065 /*  S       (input) REAL array, dimension (N) */
00066 /*          The diagonal entries of the bidiagonal matrix B. */
00067 
00068 /*  E       (input) REAL array of dimension (N-1) */
00069 /*          The superdiagonal entries of the bidiagonal matrix B. */
00070 
00071 /*  SHIFT   (input) REAL */
00072 /*          The shift, used as described under Purpose. */
00073 
00074 /*  NUM     (output) INTEGER */
00075 /*          The number of eigenvalues of T less than or equal to SHIFT. */
00076 
00077 /*  ===================================================================== */
00078 
00079 /*     .. Parameters .. */
00080 /*     .. */
00081 /*     .. Local Scalars .. */
00082 /*     .. */
00083 /*     .. External Functions .. */
00084 /*     .. */
00085 /*     .. Intrinsic Functions .. */
00086 /*     .. */
00087 /*     .. Executable Statements .. */
00088 
00089 /*     Get machine constants */
00090 
00091     /* Parameter adjustments */
00092     --e;
00093     --s;
00094 
00095     /* Function Body */
00096     unfl = slamch_("Safe minimum") * 2;
00097     ovfl = 1.f / unfl;
00098 
00099 /*     Find largest entry */
00100 
00101     mx = dabs(s[1]);
00102     i__1 = *n - 1;
00103     for (i__ = 1; i__ <= i__1; ++i__) {
00104 /* Computing MAX */
00105         r__3 = mx, r__4 = (r__1 = s[i__ + 1], dabs(r__1)), r__3 = max(r__3,
00106                 r__4), r__4 = (r__2 = e[i__], dabs(r__2));
00107         mx = dmax(r__3,r__4);
00108 /* L10: */
00109     }
00110 
00111     if (mx == 0.f) {
00112         if (*shift < 0.f) {
00113             *num = 0;
00114         } else {
00115             *num = *n << 1;
00116         }
00117         return 0;
00118     }
00119 
00120 /*     Compute scale factors as in Kahan's report */
00121 
00122     sun = sqrt(unfl);
00123     ssun = sqrt(sun);
00124     sov = sqrt(ovfl);
00125     tom = ssun * sov;
00126     if (mx <= 1.f) {
00127         m1 = 1.f / mx;
00128         m2 = tom;
00129     } else {
00130         m1 = 1.f;
00131         m2 = tom / mx;
00132     }
00133 
00134 /*     Begin counting */
00135 
00136     u = 1.f;
00137     *num = 0;
00138     sshift = *shift * m1 * m2;
00139     u = -sshift;
00140     if (u <= sun) {
00141         if (u <= 0.f) {
00142             ++(*num);
00143             if (u > -sun) {
00144                 u = -sun;
00145             }
00146         } else {
00147             u = sun;
00148         }
00149     }
00150     tmp = s[1] * m1 * m2;
00151     u = -tmp * (tmp / u) - sshift;
00152     if (u <= sun) {
00153         if (u <= 0.f) {
00154             ++(*num);
00155             if (u > -sun) {
00156                 u = -sun;
00157             }
00158         } else {
00159             u = sun;
00160         }
00161     }
00162     i__1 = *n - 1;
00163     for (i__ = 1; i__ <= i__1; ++i__) {
00164         tmp = e[i__] * m1 * m2;
00165         u = -tmp * (tmp / u) - sshift;
00166         if (u <= sun) {
00167             if (u <= 0.f) {
00168                 ++(*num);
00169                 if (u > -sun) {
00170                     u = -sun;
00171                 }
00172             } else {
00173                 u = sun;
00174             }
00175         }
00176         tmp = s[i__ + 1] * m1 * m2;
00177         u = -tmp * (tmp / u) - sshift;
00178         if (u <= sun) {
00179             if (u <= 0.f) {
00180                 ++(*num);
00181                 if (u > -sun) {
00182                     u = -sun;
00183                 }
00184             } else {
00185                 u = sun;
00186             }
00187         }
00188 /* L20: */
00189     }
00190     return 0;
00191 
00192 /*     End of SSVDCT */
00193 
00194 } /* ssvdct_ */


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