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
     157
     158
     159
     160
     161
     162
     163
     164
     165
     166
     167
     168
     169
     170
     171
     172
     173
     174
     175
     176
     177
      SUBROUTINE DSTECHNABEIGTOLWORKINFO )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            INFON
      DOUBLE PRECISION   TOL
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A* ), B* ), EIG* ), WORK* )
*     ..
*
*  Purpose
*  =======
*
*     Let T be the tridiagonal matrix with diagonal entries A(1) ,...,
*     A(N) and offdiagonal entries B(1) ,..., B(N-1)).  DSTECH checks to
*     see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T.
*     It does this by expanding each EIG(I) into an interval
*     [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if
*     any, and using Sturm sequences to count and verify whether each
*     resulting interval has the correct number of eigenvalues (using
*     DSTECT).  Here EPS = TOL*MAZHEPS*MAXEIG, where MACHEPS is the
*     machine precision and MAXEIG is the absolute value of the largest
*     eigenvalue. If each interval contains the correct number of
*     eigenvalues, INFO = 0 is returned, otherwise INFO is the index of
*     the first eigenvalue in the first bad interval.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The dimension of the tridiagonal matrix T.
*
*  A       (input) DOUBLE PRECISION array, dimension (N)
*          The diagonal entries of the tridiagonal matrix T.
*
*  B       (input) DOUBLE PRECISION array, dimension (N-1)
*          The offdiagonal entries of the tridiagonal matrix T.
*
*  EIG     (input) DOUBLE PRECISION array, dimension (N)
*          The purported eigenvalues to be checked.
*
*  TOL     (input) DOUBLE PRECISION
*          Error tolerance for checking, a multiple of the
*          machine precision.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
*
*  INFO    (output) INTEGER
*          0  if the eigenvalues are all correct (to within
*             1 +- TOL*MAZHEPS*MAXEIG)
*          >0 if the interval containing the INFO-th eigenvalue
*             contains the incorrect number of eigenvalues.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D0 )
*     ..
*     .. Local Scalars ..
      INTEGER            BPNTCOUNTIISUBJNUMLNUMUTPNT
      DOUBLE PRECISION   EMINEPSLOWERMXTUPPRUNFLEPUPPER
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           DSTECT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABSMAX
*     ..
*     .. Executable Statements ..
*
*     Check input parameters
*
      INFO = 0
      IFN.EQ.0 )
     $   RETURN
      IFN.LT.0 ) THEN
         INFO = -1
         RETURN
      END IF
      IFTOL.LT.ZERO ) THEN
         INFO = -5
         RETURN
      END IF
*
*     Get machine constants
*
      EPS = DLAMCH'Epsilon' )*DLAMCH'Base' )
      UNFLEP = DLAMCH'Safe minimum' ) / EPS
      EPS = TOL*EPS
*
*     Compute maximum absolute eigenvalue, error tolerance
*
      MX = ABSEIG1 ) )
      DO 10 I = 2N
         MX = MAXMXABSEIGI ) ) )
   10 CONTINUE
      EPS = MAXEPS*MXUNFLEP )
*
*     Sort eigenvalues from EIG into WORK
*
      DO 20 I = 1N
         WORKI ) = EIGI )
   20 CONTINUE
      DO 40 I = 1N - 1
         ISUB = 1
         EMIN = WORK1 )
         DO 30 J = 2N + 1 - I
            IFWORKJ ).LT.EMIN ) THEN
               ISUB = J
               EMIN = WORKJ )
            END IF
   30    CONTINUE
         IFISUB.NE.N+1-I ) THEN
            WORKISUB ) = WORKN+1-I )
            WORKN+1-I ) = EMIN
         END IF
   40 CONTINUE
*
*     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
*
   50 CONTINUE
      UPPER = WORKTPNT ) + EPS
      LOWER = WORKBPNT ) - EPS
*
*     Begin loop merging overlapping intervals
*
   60 CONTINUE
      IFBPNT.EQ.N )
     $   GO TO 70
      TUPPR = WORKBPNT+1 ) + EPS
      IFTUPPR.LT.LOWER )
     $   GO TO 70
*
*     Merge
*
      BPNT = BPNT + 1
      LOWER = WORKBPNT ) - EPS
      GO TO 60
   70 CONTINUE
*
*     Count singular values in interval [ LOWER, UPPER ]
*
      CALL DSTECTNABLOWERNUML )
      CALL DSTECTNABUPPERNUMU )
      COUNT = NUMU - NUML
      IFCOUNT.NE.BPNT-TPNT+1 ) THEN
*
*        Wrong number of singular values in interval
*
         INFO = TPNT
         GO TO 80
      END IF
      TPNT = BPNT + 1
      BPNT = TPNT
      IFTPNT.LE.N )
     $   GO TO 50
   80 CONTINUE
      RETURN
*
*     End of DSTECH
*
      END