LAPACK 3.3.1
Linear Algebra PACKage

dsvdch.f

Go to the documentation of this file.
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
 All Files Functions