00001 SUBROUTINE DSVDCH( N, S, E, SVD, TOL, INFO ) 00002 * 00003 * -- LAPACK test routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Scalar Arguments .. 00008 INTEGER INFO, N 00009 DOUBLE PRECISION TOL 00010 * .. 00011 * .. Array Arguments .. 00012 DOUBLE PRECISION E( * ), S( * ), SVD( * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * DSVDCH checks to see if SVD(1) ,..., SVD(N) are accurate singular 00019 * values of the bidiagonal matrix B with diagonal entries 00020 * S(1) ,..., S(N) and superdiagonal entries E(1) ,..., E(N-1)). 00021 * It does this by expanding each SVD(I) into an interval 00022 * [SVD(I) * (1-EPS) , SVD(I) * (1+EPS)], merging overlapping intervals 00023 * if any, and using Sturm sequences to count and verify whether each 00024 * resulting interval has the correct number of singular values (using 00025 * DSVDCT). Here EPS=TOL*MAX(N/10,1)*MAZHEP, where MACHEP is the 00026 * machine precision. The routine assumes the singular values are sorted 00027 * with SVD(1) the largest and SVD(N) smallest. If each interval 00028 * contains the correct number of singular values, INFO = 0 is returned, 00029 * otherwise INFO is the index of the first singular value in the first 00030 * bad interval. 00031 * 00032 * Arguments 00033 * ========== 00034 * 00035 * N (input) INTEGER 00036 * The dimension of the bidiagonal matrix B. 00037 * 00038 * S (input) DOUBLE PRECISION array, dimension (N) 00039 * The diagonal entries of the bidiagonal matrix B. 00040 * 00041 * E (input) DOUBLE PRECISION array, dimension (N-1) 00042 * The superdiagonal entries of the bidiagonal matrix B. 00043 * 00044 * SVD (input) DOUBLE PRECISION array, dimension (N) 00045 * The computed singular values to be checked. 00046 * 00047 * TOL (input) DOUBLE PRECISION 00048 * Error tolerance for checking, a multiplier of the 00049 * machine precision. 00050 * 00051 * INFO (output) INTEGER 00052 * =0 if the singular values are all correct (to within 00053 * 1 +- TOL*MAZHEPS) 00054 * >0 if the interval containing the INFO-th singular value 00055 * contains the incorrect number of singular values. 00056 * 00057 * ===================================================================== 00058 * 00059 * .. Parameters .. 00060 DOUBLE PRECISION ONE 00061 PARAMETER ( ONE = 1.0D0 ) 00062 DOUBLE PRECISION ZERO 00063 PARAMETER ( ZERO = 0.0D0 ) 00064 * .. 00065 * .. Local Scalars .. 00066 INTEGER BPNT, COUNT, NUML, NUMU, TPNT 00067 DOUBLE PRECISION EPS, LOWER, OVFL, TUPPR, UNFL, UNFLEP, UPPER 00068 * .. 00069 * .. External Functions .. 00070 DOUBLE PRECISION DLAMCH 00071 EXTERNAL DLAMCH 00072 * .. 00073 * .. External Subroutines .. 00074 EXTERNAL DSVDCT 00075 * .. 00076 * .. Intrinsic Functions .. 00077 INTRINSIC MAX, SQRT 00078 * .. 00079 * .. Executable Statements .. 00080 * 00081 * Get machine constants 00082 * 00083 INFO = 0 00084 IF( N.LE.0 ) 00085 $ RETURN 00086 UNFL = DLAMCH( 'Safe minimum' ) 00087 OVFL = DLAMCH( 'Overflow' ) 00088 EPS = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00089 * 00090 * UNFLEP is chosen so that when an eigenvalue is multiplied by the 00091 * scale factor sqrt(OVFL)*sqrt(sqrt(UNFL))/MX in DSVDCT, it exceeds 00092 * sqrt(UNFL), which is the lower limit for DSVDCT. 00093 * 00094 UNFLEP = ( SQRT( SQRT( UNFL ) ) / SQRT( OVFL ) )*SVD( 1 ) + 00095 $ UNFL / EPS 00096 * 00097 * The value of EPS works best when TOL .GE. 10. 00098 * 00099 EPS = TOL*MAX( N / 10, 1 )*EPS 00100 * 00101 * TPNT points to singular value at right endpoint of interval 00102 * BPNT points to singular value at left endpoint of interval 00103 * 00104 TPNT = 1 00105 BPNT = 1 00106 * 00107 * Begin loop over all intervals 00108 * 00109 10 CONTINUE 00110 UPPER = ( ONE+EPS )*SVD( TPNT ) + UNFLEP 00111 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP 00112 IF( LOWER.LE.UNFLEP ) 00113 $ LOWER = -UPPER 00114 * 00115 * Begin loop merging overlapping intervals 00116 * 00117 20 CONTINUE 00118 IF( BPNT.EQ.N ) 00119 $ GO TO 30 00120 TUPPR = ( ONE+EPS )*SVD( BPNT+1 ) + UNFLEP 00121 IF( TUPPR.LT.LOWER ) 00122 $ GO TO 30 00123 * 00124 * Merge 00125 * 00126 BPNT = BPNT + 1 00127 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP 00128 IF( LOWER.LE.UNFLEP ) 00129 $ LOWER = -UPPER 00130 GO TO 20 00131 30 CONTINUE 00132 * 00133 * Count singular values in interval [ LOWER, UPPER ] 00134 * 00135 CALL DSVDCT( N, S, E, LOWER, NUML ) 00136 CALL DSVDCT( N, S, E, UPPER, NUMU ) 00137 COUNT = NUMU - NUML 00138 IF( LOWER.LT.ZERO ) 00139 $ COUNT = COUNT / 2 00140 IF( COUNT.NE.BPNT-TPNT+1 ) THEN 00141 * 00142 * Wrong number of singular values in interval 00143 * 00144 INFO = TPNT 00145 GO TO 40 00146 END IF 00147 TPNT = BPNT + 1 00148 BPNT = TPNT 00149 IF( TPNT.LE.N ) 00150 $ GO TO 10 00151 40 CONTINUE 00152 RETURN 00153 * 00154 * End of DSVDCH 00155 * 00156 END
1.7.2