Browse Files
       1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141
     142
     143
     144
     145
     146
     147
     148
     149
     150
     151
     152
     153
     154
     155
     156
      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