ssvdch.c
Go to the documentation of this file.
00001 /* ssvdch.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 ssvdch_(integer *n, real *s, real *e, real *svd, real *
00017         tol, integer *info)
00018 {
00019     /* System generated locals */
00020     integer i__1;
00021 
00022     /* Builtin functions */
00023     double sqrt(doublereal);
00024 
00025     /* Local variables */
00026     real eps;
00027     integer bpnt;
00028     real unfl, ovfl;
00029     integer numl, numu, tpnt, count;
00030     real lower, upper, tuppr;
00031     extern doublereal slamch_(char *);
00032     real unflep;
00033     extern /* Subroutine */ int ssvdct_(integer *, real *, real *, real *, 
00034             integer *);
00035 
00036 
00037 /*  -- LAPACK test routine (version 3.1) -- */
00038 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00039 /*     November 2006 */
00040 
00041 /*     .. Scalar Arguments .. */
00042 /*     .. */
00043 /*     .. Array Arguments .. */
00044 /*     .. */
00045 
00046 /*  Purpose */
00047 /*  ======= */
00048 
00049 /*  SSVDCH checks to see if SVD(1) ,..., SVD(N) are accurate singular */
00050 /*  values of the bidiagonal matrix B with diagonal entries */
00051 /*  S(1) ,..., S(N) and superdiagonal entries E(1) ,..., E(N-1)). */
00052 /*  It does this by expanding each SVD(I) into an interval */
00053 /*  [SVD(I) * (1-EPS) , SVD(I) * (1+EPS)], merging overlapping intervals */
00054 /*  if any, and using Sturm sequences to count and verify whether each */
00055 /*  resulting interval has the correct number of singular values (using */
00056 /*  SSVDCT). Here EPS=TOL*MAX(N/10,1)*MACHEP, where MACHEP is the */
00057 /*  machine precision. The routine assumes the singular values are sorted */
00058 /*  with SVD(1) the largest and SVD(N) smallest.  If each interval */
00059 /*  contains the correct number of singular values, INFO = 0 is returned, */
00060 /*  otherwise INFO is the index of the first singular value in the first */
00061 /*  bad interval. */
00062 
00063 /*  Arguments */
00064 /*  ========== */
00065 
00066 /*  N       (input) INTEGER */
00067 /*          The dimension of the bidiagonal matrix B. */
00068 
00069 /*  S       (input) REAL array, dimension (N) */
00070 /*          The diagonal entries of the bidiagonal matrix B. */
00071 
00072 /*  E       (input) REAL array, dimension (N-1) */
00073 /*          The superdiagonal entries of the bidiagonal matrix B. */
00074 
00075 /*  SVD     (input) REAL array, dimension (N) */
00076 /*          The computed singular values to be checked. */
00077 
00078 /*  TOL     (input) REAL */
00079 /*          Error tolerance for checking, a multiplier of the */
00080 /*          machine precision. */
00081 
00082 /*  INFO    (output) INTEGER */
00083 /*          =0 if the singular values are all correct (to within */
00084 /*             1 +- TOL*MACHEPS) */
00085 /*          >0 if the interval containing the INFO-th singular value */
00086 /*             contains the incorrect number of singular values. */
00087 
00088 /*  ===================================================================== */
00089 
00090 /*     .. Parameters .. */
00091 /*     .. */
00092 /*     .. Local Scalars .. */
00093 /*     .. */
00094 /*     .. External Functions .. */
00095 /*     .. */
00096 /*     .. External Subroutines .. */
00097 /*     .. */
00098 /*     .. Intrinsic Functions .. */
00099 /*     .. */
00100 /*     .. Executable Statements .. */
00101 
00102 /*     Get machine constants */
00103 
00104     /* Parameter adjustments */
00105     --svd;
00106     --e;
00107     --s;
00108 
00109     /* Function Body */
00110     *info = 0;
00111     if (*n <= 0) {
00112         return 0;
00113     }
00114     unfl = slamch_("Safe minimum");
00115     ovfl = slamch_("Overflow");
00116     eps = slamch_("Epsilon") * slamch_("Base");
00117 
00118 /*     UNFLEP is chosen so that when an eigenvalue is multiplied by the */
00119 /*     scale factor sqrt(OVFL)*sqrt(sqrt(UNFL))/MX in SSVDCT, it exceeds */
00120 /*     sqrt(UNFL), which is the lower limit for SSVDCT. */
00121 
00122     unflep = sqrt(sqrt(unfl)) / sqrt(ovfl) * svd[1] + unfl / eps;
00123 
00124 /*     The value of EPS works best when TOL .GE. 10. */
00125 
00126 /* Computing MAX */
00127     i__1 = *n / 10;
00128     eps = *tol * max(i__1,1) * eps;
00129 
00130 /*     TPNT points to singular value at right endpoint of interval */
00131 /*     BPNT points to singular value at left  endpoint of interval */
00132 
00133     tpnt = 1;
00134     bpnt = 1;
00135 
00136 /*     Begin loop over all intervals */
00137 
00138 L10:
00139     upper = (eps + 1.f) * svd[tpnt] + unflep;
00140     lower = (1.f - eps) * svd[bpnt] - unflep;
00141     if (lower <= unflep) {
00142         lower = -upper;
00143     }
00144 
00145 /*     Begin loop merging overlapping intervals */
00146 
00147 L20:
00148     if (bpnt == *n) {
00149         goto L30;
00150     }
00151     tuppr = (eps + 1.f) * svd[bpnt + 1] + unflep;
00152     if (tuppr < lower) {
00153         goto L30;
00154     }
00155 
00156 /*     Merge */
00157 
00158     ++bpnt;
00159     lower = (1.f - eps) * svd[bpnt] - unflep;
00160     if (lower <= unflep) {
00161         lower = -upper;
00162     }
00163     goto L20;
00164 L30:
00165 
00166 /*     Count singular values in interval [ LOWER, UPPER ] */
00167 
00168     ssvdct_(n, &s[1], &e[1], &lower, &numl);
00169     ssvdct_(n, &s[1], &e[1], &upper, &numu);
00170     count = numu - numl;
00171     if (lower < 0.f) {
00172         count /= 2;
00173     }
00174     if (count != bpnt - tpnt + 1) {
00175 
00176 /*        Wrong number of singular values in interval */
00177 
00178         *info = tpnt;
00179         goto L40;
00180     }
00181     tpnt = bpnt + 1;
00182     bpnt = tpnt;
00183     if (tpnt <= *n) {
00184         goto L10;
00185     }
00186 L40:
00187     return 0;
00188 
00189 /*     End of SSVDCH */
00190 
00191 } /* ssvdch_ */


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