LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dstech.f
Go to the documentation of this file.
00001 *> \brief \b DSTECH
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DSTECH( N, A, B, EIG, TOL, WORK, INFO )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       INTEGER            INFO, N
00015 *       DOUBLE PRECISION   TOL
00016 *       ..
00017 *       .. Array Arguments ..
00018 *       DOUBLE PRECISION   A( * ), B( * ), EIG( * ), WORK( * )
00019 *       ..
00020 *  
00021 *
00022 *> \par Purpose:
00023 *  =============
00024 *>
00025 *> \verbatim
00026 *>
00027 *>    Let T be the tridiagonal matrix with diagonal entries A(1) ,...,
00028 *>    A(N) and offdiagonal entries B(1) ,..., B(N-1)).  DSTECH checks to
00029 *>    see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T.
00030 *>    It does this by expanding each EIG(I) into an interval
00031 *>    [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if
00032 *>    any, and using Sturm sequences to count and verify whether each
00033 *>    resulting interval has the correct number of eigenvalues (using
00034 *>    DSTECT).  Here EPS = TOL*MAZHEPS*MAXEIG, where MACHEPS is the
00035 *>    machine precision and MAXEIG is the absolute value of the largest
00036 *>    eigenvalue. If each interval contains the correct number of
00037 *>    eigenvalues, INFO = 0 is returned, otherwise INFO is the index of
00038 *>    the first eigenvalue in the first bad interval.
00039 *> \endverbatim
00040 *
00041 *  Arguments:
00042 *  ==========
00043 *
00044 *> \param[in] N
00045 *> \verbatim
00046 *>          N is INTEGER
00047 *>          The dimension of the tridiagonal matrix T.
00048 *> \endverbatim
00049 *>
00050 *> \param[in] A
00051 *> \verbatim
00052 *>          A is DOUBLE PRECISION array, dimension (N)
00053 *>          The diagonal entries of the tridiagonal matrix T.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] B
00057 *> \verbatim
00058 *>          B is DOUBLE PRECISION array, dimension (N-1)
00059 *>          The offdiagonal entries of the tridiagonal matrix T.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] EIG
00063 *> \verbatim
00064 *>          EIG is DOUBLE PRECISION array, dimension (N)
00065 *>          The purported eigenvalues to be checked.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] TOL
00069 *> \verbatim
00070 *>          TOL is DOUBLE PRECISION
00071 *>          Error tolerance for checking, a multiple of the
00072 *>          machine precision.
00073 *> \endverbatim
00074 *>
00075 *> \param[out] WORK
00076 *> \verbatim
00077 *>          WORK is DOUBLE PRECISION array, dimension (N)
00078 *> \endverbatim
00079 *>
00080 *> \param[out] INFO
00081 *> \verbatim
00082 *>          INFO is INTEGER
00083 *>          0  if the eigenvalues are all correct (to within
00084 *>             1 +- TOL*MAZHEPS*MAXEIG)
00085 *>          >0 if the interval containing the INFO-th eigenvalue
00086 *>             contains the incorrect number of eigenvalues.
00087 *> \endverbatim
00088 *
00089 *  Authors:
00090 *  ========
00091 *
00092 *> \author Univ. of Tennessee 
00093 *> \author Univ. of California Berkeley 
00094 *> \author Univ. of Colorado Denver 
00095 *> \author NAG Ltd. 
00096 *
00097 *> \date November 2011
00098 *
00099 *> \ingroup double_eig
00100 *
00101 *  =====================================================================
00102       SUBROUTINE DSTECH( N, A, B, EIG, TOL, WORK, INFO )
00103 *
00104 *  -- LAPACK test routine (version 3.4.0) --
00105 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00106 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00107 *     November 2011
00108 *
00109 *     .. Scalar Arguments ..
00110       INTEGER            INFO, N
00111       DOUBLE PRECISION   TOL
00112 *     ..
00113 *     .. Array Arguments ..
00114       DOUBLE PRECISION   A( * ), B( * ), EIG( * ), WORK( * )
00115 *     ..
00116 *
00117 *  =====================================================================
00118 *
00119 *     .. Parameters ..
00120       DOUBLE PRECISION   ZERO
00121       PARAMETER          ( ZERO = 0.0D0 )
00122 *     ..
00123 *     .. Local Scalars ..
00124       INTEGER            BPNT, COUNT, I, ISUB, J, NUML, NUMU, TPNT
00125       DOUBLE PRECISION   EMIN, EPS, LOWER, MX, TUPPR, UNFLEP, UPPER
00126 *     ..
00127 *     .. External Functions ..
00128       DOUBLE PRECISION   DLAMCH
00129       EXTERNAL           DLAMCH
00130 *     ..
00131 *     .. External Subroutines ..
00132       EXTERNAL           DSTECT
00133 *     ..
00134 *     .. Intrinsic Functions ..
00135       INTRINSIC          ABS, MAX
00136 *     ..
00137 *     .. Executable Statements ..
00138 *
00139 *     Check input parameters
00140 *
00141       INFO = 0
00142       IF( N.EQ.0 )
00143      $   RETURN
00144       IF( N.LT.0 ) THEN
00145          INFO = -1
00146          RETURN
00147       END IF
00148       IF( TOL.LT.ZERO ) THEN
00149          INFO = -5
00150          RETURN
00151       END IF
00152 *
00153 *     Get machine constants
00154 *
00155       EPS = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00156       UNFLEP = DLAMCH( 'Safe minimum' ) / EPS
00157       EPS = TOL*EPS
00158 *
00159 *     Compute maximum absolute eigenvalue, error tolerance
00160 *
00161       MX = ABS( EIG( 1 ) )
00162       DO 10 I = 2, N
00163          MX = MAX( MX, ABS( EIG( I ) ) )
00164    10 CONTINUE
00165       EPS = MAX( EPS*MX, UNFLEP )
00166 *
00167 *     Sort eigenvalues from EIG into WORK
00168 *
00169       DO 20 I = 1, N
00170          WORK( I ) = EIG( I )
00171    20 CONTINUE
00172       DO 40 I = 1, N - 1
00173          ISUB = 1
00174          EMIN = WORK( 1 )
00175          DO 30 J = 2, N + 1 - I
00176             IF( WORK( J ).LT.EMIN ) THEN
00177                ISUB = J
00178                EMIN = WORK( J )
00179             END IF
00180    30    CONTINUE
00181          IF( ISUB.NE.N+1-I ) THEN
00182             WORK( ISUB ) = WORK( N+1-I )
00183             WORK( N+1-I ) = EMIN
00184          END IF
00185    40 CONTINUE
00186 *
00187 *     TPNT points to singular value at right endpoint of interval
00188 *     BPNT points to singular value at left  endpoint of interval
00189 *
00190       TPNT = 1
00191       BPNT = 1
00192 *
00193 *     Begin loop over all intervals
00194 *
00195    50 CONTINUE
00196       UPPER = WORK( TPNT ) + EPS
00197       LOWER = WORK( BPNT ) - EPS
00198 *
00199 *     Begin loop merging overlapping intervals
00200 *
00201    60 CONTINUE
00202       IF( BPNT.EQ.N )
00203      $   GO TO 70
00204       TUPPR = WORK( BPNT+1 ) + EPS
00205       IF( TUPPR.LT.LOWER )
00206      $   GO TO 70
00207 *
00208 *     Merge
00209 *
00210       BPNT = BPNT + 1
00211       LOWER = WORK( BPNT ) - EPS
00212       GO TO 60
00213    70 CONTINUE
00214 *
00215 *     Count singular values in interval [ LOWER, UPPER ]
00216 *
00217       CALL DSTECT( N, A, B, LOWER, NUML )
00218       CALL DSTECT( N, A, B, UPPER, NUMU )
00219       COUNT = NUMU - NUML
00220       IF( COUNT.NE.BPNT-TPNT+1 ) THEN
00221 *
00222 *        Wrong number of singular values in interval
00223 *
00224          INFO = TPNT
00225          GO TO 80
00226       END IF
00227       TPNT = BPNT + 1
00228       BPNT = TPNT
00229       IF( TPNT.LE.N )
00230      $   GO TO 50
00231    80 CONTINUE
00232       RETURN
00233 *
00234 *     End of DSTECH
00235 *
00236       END
 All Files Functions