dsvdch.f

Go to the documentation of this file.

      SUBROUTINE DSVDCH( N, S, E, SVD, TOL, INFO )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            INFO, N
      DOUBLE PRECISION   TOL
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   E( * ), S( * ), SVD( * )
*     ..
*
*  Purpose
*  =======
*
*  DSVDCH checks to see if SVD(1) ,..., SVD(N) are accurate singular
*  values of the bidiagonal matrix B with diagonal entries
*  S(1) ,..., S(N) and superdiagonal entries E(1) ,..., E(N-1)).
*  It does this by expanding each SVD(I) into an interval
*  [SVD(I) * (1-EPS) , SVD(I) * (1+EPS)], merging overlapping intervals
*  if any, and using Sturm sequences to count and verify whether each
*  resulting interval has the correct number of singular values (using
*  DSVDCT). Here EPS=TOL*MAX(N/10,1)*MAZHEP, where MACHEP is the
*  machine precision. The routine assumes the singular values are sorted
*  with SVD(1) the largest and SVD(N) smallest.  If each interval
*  contains the correct number of singular values, INFO = 0 is returned,
*  otherwise INFO is the index of the first singular value in the first
*  bad interval.
*
*  Arguments
*  ==========
*
*  N       (input) INTEGER
*          The dimension of the bidiagonal matrix B.
*
*  S       (input) DOUBLE PRECISION array, dimension (N)
*          The diagonal entries of the bidiagonal matrix B.
*
*  E       (input) DOUBLE PRECISION array, dimension (N-1)
*          The superdiagonal entries of the bidiagonal matrix B.
*
*  SVD     (input) DOUBLE PRECISION array, dimension (N)
*          The computed singular values to be checked.
*
*  TOL     (input) DOUBLE PRECISION
*          Error tolerance for checking, a multiplier of the
*          machine precision.
*
*  INFO    (output) INTEGER
*          =0 if the singular values are all correct (to within
*             1 +- TOL*MAZHEPS)
*          >0 if the interval containing the INFO-th singular value
*             contains the incorrect number of singular values.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D0 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D0 )
*     ..
*     .. Local Scalars ..
      INTEGER            BPNT, COUNT, NUML, NUMU, TPNT
      DOUBLE PRECISION   EPS, LOWER, OVFL, TUPPR, UNFL, UNFLEP, UPPER
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           DSVDCT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, SQRT
*     ..
*     .. Executable Statements ..
*
*     Get machine constants
*
      INFO = 0
      IF( N.LE.0 )
     $   RETURN
      UNFL = DLAMCH( 'Safe minimum' )
      OVFL = DLAMCH( 'Overflow' )
      EPS = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
*
*     UNFLEP is chosen so that when an eigenvalue is multiplied by the
*     scale factor sqrt(OVFL)*sqrt(sqrt(UNFL))/MX in DSVDCT, it exceeds
*     sqrt(UNFL), which is the lower limit for DSVDCT.
*
      UNFLEP = ( SQRT( SQRT( UNFL ) ) / SQRT( OVFL ) )*SVD( 1 ) +
     $         UNFL / EPS
*
*     The value of EPS works best when TOL .GE. 10.
*
      EPS = TOL*MAX( N / 10, 1 )*EPS
*
*     TPNT points to singular value at right endpoint of interval
*     BPNT points to singular value at left  endpoint of interval
*
      TPNT = 1
      BPNT = 1
*
*     Begin loop over all intervals
*
   10 CONTINUE
      UPPER = ( ONE+EPS )*SVD( TPNT ) + UNFLEP
      LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP
      IF( LOWER.LE.UNFLEP )
     $   LOWER = -UPPER
*
*     Begin loop merging overlapping intervals
*
   20 CONTINUE
      IF( BPNT.EQ.N )
     $   GO TO 30
      TUPPR = ( ONE+EPS )*SVD( BPNT+1 ) + UNFLEP
      IF( TUPPR.LT.LOWER )
     $   GO TO 30
*
*     Merge
*
      BPNT = BPNT + 1
      LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP
      IF( LOWER.LE.UNFLEP )
     $   LOWER = -UPPER
      GO TO 20
   30 CONTINUE
*
*     Count singular values in interval [ LOWER, UPPER ]
*
      CALL DSVDCT( N, S, E, LOWER, NUML )
      CALL DSVDCT( N, S, E, UPPER, NUMU )
      COUNT = NUMU - NUML
      IF( LOWER.LT.ZERO )
     $   COUNT = COUNT / 2
      IF( COUNT.NE.BPNT-TPNT+1 ) THEN
*
*        Wrong number of singular values in interval
*
         INFO = TPNT
         GO TO 40
      END IF
      TPNT = BPNT + 1
      BPNT = TPNT
      IF( TPNT.LE.N )
     $   GO TO 10
   40 CONTINUE
      RETURN
*
*     End of DSVDCH
*
      END