![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSVDCH 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE DSVDCH( N, S, E, SVD, TOL, INFO ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER INFO, N 00015 * DOUBLE PRECISION TOL 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION E( * ), S( * ), SVD( * ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> DSVDCH checks to see if SVD(1) ,..., SVD(N) are accurate singular 00028 *> values of the bidiagonal matrix B with diagonal entries 00029 *> S(1) ,..., S(N) and superdiagonal entries E(1) ,..., E(N-1)). 00030 *> It does this by expanding each SVD(I) into an interval 00031 *> [SVD(I) * (1-EPS) , SVD(I) * (1+EPS)], merging overlapping intervals 00032 *> if any, and using Sturm sequences to count and verify whether each 00033 *> resulting interval has the correct number of singular values (using 00034 *> DSVDCT). Here EPS=TOL*MAX(N/10,1)*MAZHEP, where MACHEP is the 00035 *> machine precision. The routine assumes the singular values are sorted 00036 *> with SVD(1) the largest and SVD(N) smallest. If each interval 00037 *> contains the correct number of singular values, INFO = 0 is returned, 00038 *> otherwise INFO is the index of the first singular value in the first 00039 *> bad interval. 00040 *> \endverbatim 00041 * 00042 * Arguments: 00043 * ========== 00044 * 00045 *> \param[in] N 00046 *> \verbatim 00047 *> N is INTEGER 00048 *> The dimension of the bidiagonal matrix B. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] S 00052 *> \verbatim 00053 *> S is DOUBLE PRECISION array, dimension (N) 00054 *> The diagonal entries of the bidiagonal matrix B. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] E 00058 *> \verbatim 00059 *> E is DOUBLE PRECISION array, dimension (N-1) 00060 *> The superdiagonal entries of the bidiagonal matrix B. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] SVD 00064 *> \verbatim 00065 *> SVD is DOUBLE PRECISION array, dimension (N) 00066 *> The computed singular values to be checked. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] TOL 00070 *> \verbatim 00071 *> TOL is DOUBLE PRECISION 00072 *> Error tolerance for checking, a multiplier of the 00073 *> machine precision. 00074 *> \endverbatim 00075 *> 00076 *> \param[out] INFO 00077 *> \verbatim 00078 *> INFO is INTEGER 00079 *> =0 if the singular values are all correct (to within 00080 *> 1 +- TOL*MAZHEPS) 00081 *> >0 if the interval containing the INFO-th singular value 00082 *> contains the incorrect number of singular values. 00083 *> \endverbatim 00084 * 00085 * Authors: 00086 * ======== 00087 * 00088 *> \author Univ. of Tennessee 00089 *> \author Univ. of California Berkeley 00090 *> \author Univ. of Colorado Denver 00091 *> \author NAG Ltd. 00092 * 00093 *> \date November 2011 00094 * 00095 *> \ingroup double_eig 00096 * 00097 * ===================================================================== 00098 SUBROUTINE DSVDCH( N, S, E, SVD, TOL, INFO ) 00099 * 00100 * -- LAPACK test routine (version 3.4.0) -- 00101 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00102 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00103 * November 2011 00104 * 00105 * .. Scalar Arguments .. 00106 INTEGER INFO, N 00107 DOUBLE PRECISION TOL 00108 * .. 00109 * .. Array Arguments .. 00110 DOUBLE PRECISION E( * ), S( * ), SVD( * ) 00111 * .. 00112 * 00113 * ===================================================================== 00114 * 00115 * .. Parameters .. 00116 DOUBLE PRECISION ONE 00117 PARAMETER ( ONE = 1.0D0 ) 00118 DOUBLE PRECISION ZERO 00119 PARAMETER ( ZERO = 0.0D0 ) 00120 * .. 00121 * .. Local Scalars .. 00122 INTEGER BPNT, COUNT, NUML, NUMU, TPNT 00123 DOUBLE PRECISION EPS, LOWER, OVFL, TUPPR, UNFL, UNFLEP, UPPER 00124 * .. 00125 * .. External Functions .. 00126 DOUBLE PRECISION DLAMCH 00127 EXTERNAL DLAMCH 00128 * .. 00129 * .. External Subroutines .. 00130 EXTERNAL DSVDCT 00131 * .. 00132 * .. Intrinsic Functions .. 00133 INTRINSIC MAX, SQRT 00134 * .. 00135 * .. Executable Statements .. 00136 * 00137 * Get machine constants 00138 * 00139 INFO = 0 00140 IF( N.LE.0 ) 00141 $ RETURN 00142 UNFL = DLAMCH( 'Safe minimum' ) 00143 OVFL = DLAMCH( 'Overflow' ) 00144 EPS = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00145 * 00146 * UNFLEP is chosen so that when an eigenvalue is multiplied by the 00147 * scale factor sqrt(OVFL)*sqrt(sqrt(UNFL))/MX in DSVDCT, it exceeds 00148 * sqrt(UNFL), which is the lower limit for DSVDCT. 00149 * 00150 UNFLEP = ( SQRT( SQRT( UNFL ) ) / SQRT( OVFL ) )*SVD( 1 ) + 00151 $ UNFL / EPS 00152 * 00153 * The value of EPS works best when TOL .GE. 10. 00154 * 00155 EPS = TOL*MAX( N / 10, 1 )*EPS 00156 * 00157 * TPNT points to singular value at right endpoint of interval 00158 * BPNT points to singular value at left endpoint of interval 00159 * 00160 TPNT = 1 00161 BPNT = 1 00162 * 00163 * Begin loop over all intervals 00164 * 00165 10 CONTINUE 00166 UPPER = ( ONE+EPS )*SVD( TPNT ) + UNFLEP 00167 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP 00168 IF( LOWER.LE.UNFLEP ) 00169 $ LOWER = -UPPER 00170 * 00171 * Begin loop merging overlapping intervals 00172 * 00173 20 CONTINUE 00174 IF( BPNT.EQ.N ) 00175 $ GO TO 30 00176 TUPPR = ( ONE+EPS )*SVD( BPNT+1 ) + UNFLEP 00177 IF( TUPPR.LT.LOWER ) 00178 $ GO TO 30 00179 * 00180 * Merge 00181 * 00182 BPNT = BPNT + 1 00183 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP 00184 IF( LOWER.LE.UNFLEP ) 00185 $ LOWER = -UPPER 00186 GO TO 20 00187 30 CONTINUE 00188 * 00189 * Count singular values in interval [ LOWER, UPPER ] 00190 * 00191 CALL DSVDCT( N, S, E, LOWER, NUML ) 00192 CALL DSVDCT( N, S, E, UPPER, NUMU ) 00193 COUNT = NUMU - NUML 00194 IF( LOWER.LT.ZERO ) 00195 $ COUNT = COUNT / 2 00196 IF( COUNT.NE.BPNT-TPNT+1 ) THEN 00197 * 00198 * Wrong number of singular values in interval 00199 * 00200 INFO = TPNT 00201 GO TO 40 00202 END IF 00203 TPNT = BPNT + 1 00204 BPNT = TPNT 00205 IF( TPNT.LE.N ) 00206 $ GO TO 10 00207 40 CONTINUE 00208 RETURN 00209 * 00210 * End of DSVDCH 00211 * 00212 END