LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dsvdch.f
Go to the documentation of this file.
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
 All Files Functions