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_ */