LAPACK: TESTING/EIG/dstech.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSTECH( N, A, B, EIG, TOL, WORK, 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   A( * ), B( * ), EIG( * ), WORK( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *     Let T be the tridiagonal matrix with diagonal entries A(1) ,...,
00019 *     A(N) and offdiagonal entries B(1) ,..., B(N-1)).  DSTECH checks to
00020 *     see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T.
00021 *     It does this by expanding each EIG(I) into an interval
00022 *     [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if
00023 *     any, and using Sturm sequences to count and verify whether each
00024 *     resulting interval has the correct number of eigenvalues (using
00025 *     DSTECT).  Here EPS = TOL*MAZHEPS*MAXEIG, where MACHEPS is the
00026 *     machine precision and MAXEIG is the absolute value of the largest
00027 *     eigenvalue. If each interval contains the correct number of
00028 *     eigenvalues, INFO = 0 is returned, otherwise INFO is the index of
00029 *     the first eigenvalue in the first bad interval.
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  N       (input) INTEGER
00035 *          The dimension of the tridiagonal matrix T.
00036 *
00037 *  A       (input) DOUBLE PRECISION array, dimension (N)
00038 *          The diagonal entries of the tridiagonal matrix T.
00039 *
00040 *  B       (input) DOUBLE PRECISION array, dimension (N-1)
00041 *          The offdiagonal entries of the tridiagonal matrix T.
00042 *
00043 *  EIG     (input) DOUBLE PRECISION array, dimension (N)
00044 *          The purported eigenvalues to be checked.
00045 *
00046 *  TOL     (input) DOUBLE PRECISION
00047 *          Error tolerance for checking, a multiple of the
00048 *          machine precision.
00049 *
00050 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00051 *
00052 *  INFO    (output) INTEGER
00053 *          0  if the eigenvalues are all correct (to within
00054 *             1 +- TOL*MAZHEPS*MAXEIG)
00055 *          >0 if the interval containing the INFO-th eigenvalue
00056 *             contains the incorrect number of eigenvalues.
00057 *
00058 *  =====================================================================
00059 *
00060 *     .. Parameters ..
00061       DOUBLE PRECISION   ZERO
00062       PARAMETER          ( ZERO = 0.0D0 )
00063 *     ..
00064 *     .. Local Scalars ..
00065       INTEGER            BPNT, COUNT, I, ISUB, J, NUML, NUMU, TPNT
00066       DOUBLE PRECISION   EMIN, EPS, LOWER, MX, TUPPR, UNFLEP, UPPER
00067 *     ..
00068 *     .. External Functions ..
00069       DOUBLE PRECISION   DLAMCH
00070       EXTERNAL           DLAMCH
00071 *     ..
00072 *     .. External Subroutines ..
00073       EXTERNAL           DSTECT
00074 *     ..
00075 *     .. Intrinsic Functions ..
00076       INTRINSIC          ABS, MAX
00077 *     ..
00078 *     .. Executable Statements ..
00079 *
00080 *     Check input parameters
00081 *
00082       INFO = 0
00083       IF( N.EQ.0 )
00084      $   RETURN
00085       IF( N.LT.0 ) THEN
00086          INFO = -1
00087          RETURN
00088       END IF
00089       IF( TOL.LT.ZERO ) THEN
00090          INFO = -5
00091          RETURN
00092       END IF
00093 *
00094 *     Get machine constants
00095 *
00096       EPS = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00097       UNFLEP = DLAMCH( 'Safe minimum' ) / EPS
00098       EPS = TOL*EPS
00099 *
00100 *     Compute maximum absolute eigenvalue, error tolerance
00101 *
00102       MX = ABS( EIG( 1 ) )
00103       DO 10 I = 2, N
00104          MX = MAX( MX, ABS( EIG( I ) ) )
00105    10 CONTINUE
00106       EPS = MAX( EPS*MX, UNFLEP )
00107 *
00108 *     Sort eigenvalues from EIG into WORK
00109 *
00110       DO 20 I = 1, N
00111          WORK( I ) = EIG( I )
00112    20 CONTINUE
00113       DO 40 I = 1, N - 1
00114          ISUB = 1
00115          EMIN = WORK( 1 )
00116          DO 30 J = 2, N + 1 - I
00117             IF( WORK( J ).LT.EMIN ) THEN
00118                ISUB = J
00119                EMIN = WORK( J )
00120             END IF
00121    30    CONTINUE
00122          IF( ISUB.NE.N+1-I ) THEN
00123             WORK( ISUB ) = WORK( N+1-I )
00124             WORK( N+1-I ) = EMIN
00125          END IF
00126    40 CONTINUE
00127 *
00128 *     TPNT points to singular value at right endpoint of interval
00129 *     BPNT points to singular value at left  endpoint of interval
00130 *
00131       TPNT = 1
00132       BPNT = 1
00133 *
00134 *     Begin loop over all intervals
00135 *
00136    50 CONTINUE
00137       UPPER = WORK( TPNT ) + EPS
00138       LOWER = WORK( BPNT ) - EPS
00139 *
00140 *     Begin loop merging overlapping intervals
00141 *
00142    60 CONTINUE
00143       IF( BPNT.EQ.N )
00144      $   GO TO 70
00145       TUPPR = WORK( BPNT+1 ) + EPS
00146       IF( TUPPR.LT.LOWER )
00147      $   GO TO 70
00148 *
00149 *     Merge
00150 *
00151       BPNT = BPNT + 1
00152       LOWER = WORK( BPNT ) - EPS
00153       GO TO 60
00154    70 CONTINUE
00155 *
00156 *     Count singular values in interval [ LOWER, UPPER ]
00157 *
00158       CALL DSTECT( N, A, B, LOWER, NUML )
00159       CALL DSTECT( N, A, B, UPPER, NUMU )
00160       COUNT = NUMU - NUML
00161       IF( COUNT.NE.BPNT-TPNT+1 ) THEN
00162 *
00163 *        Wrong number of singular values in interval
00164 *
00165          INFO = TPNT
00166          GO TO 80
00167       END IF
00168       TPNT = BPNT + 1
00169       BPNT = TPNT
00170       IF( TPNT.LE.N )
00171      $   GO TO 50
00172    80 CONTINUE
00173       RETURN
00174 *
00175 *     End of DSTECH
00176 *
00177       END