LAPACK 3.3.0

dstect.f

Go to the documentation of this file.
00001       SUBROUTINE DSTECT( N, A, B, SHIFT, NUM )
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            N, NUM
00009       DOUBLE PRECISION   SHIFT
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   A( * ), B( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *     DSTECT counts the number NUM of eigenvalues of a tridiagonal
00019 *     matrix T which are less than or equal to SHIFT. T has
00020 *     diagonal entries A(1), ... , A(N), and offdiagonal entries
00021 *     B(1), ..., B(N-1).
00022 *     See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00023 *     Matrix", Report CS41, Computer Science Dept., Stanford
00024 *     University, July 21, 1966
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  N       (input) INTEGER
00030 *          The dimension of the tridiagonal matrix T.
00031 *
00032 *  A       (input) DOUBLE PRECISION array, dimension (N)
00033 *          The diagonal entries of the tridiagonal matrix T.
00034 *
00035 *  B       (input) DOUBLE PRECISION array, dimension (N-1)
00036 *          The offdiagonal entries of the tridiagonal matrix T.
00037 *
00038 *  SHIFT   (input) DOUBLE PRECISION
00039 *          The shift, used as described under Purpose.
00040 *
00041 *  NUM     (output) INTEGER
00042 *          The number of eigenvalues of T less than or equal
00043 *          to SHIFT.
00044 *
00045 *  =====================================================================
00046 *
00047 *     .. Parameters ..
00048       DOUBLE PRECISION   ZERO, ONE, THREE
00049       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, THREE = 3.0D0 )
00050 *     ..
00051 *     .. Local Scalars ..
00052       INTEGER            I
00053       DOUBLE PRECISION   M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
00054      $                   TOM, U, UNFL
00055 *     ..
00056 *     .. External Functions ..
00057       DOUBLE PRECISION   DLAMCH
00058       EXTERNAL           DLAMCH
00059 *     ..
00060 *     .. Intrinsic Functions ..
00061       INTRINSIC          ABS, MAX, SQRT
00062 *     ..
00063 *     .. Executable Statements ..
00064 *
00065 *     Get machine constants
00066 *
00067       UNFL = DLAMCH( 'Safe minimum' )
00068       OVFL = DLAMCH( 'Overflow' )
00069 *
00070 *     Find largest entry
00071 *
00072       MX = ABS( A( 1 ) )
00073       DO 10 I = 1, N - 1
00074          MX = MAX( MX, ABS( A( I+1 ) ), ABS( B( I ) ) )
00075    10 CONTINUE
00076 *
00077 *     Handle easy cases, including zero matrix
00078 *
00079       IF( SHIFT.GE.THREE*MX ) THEN
00080          NUM = N
00081          RETURN
00082       END IF
00083       IF( SHIFT.LT.-THREE*MX ) THEN
00084          NUM = 0
00085          RETURN
00086       END IF
00087 *
00088 *     Compute scale factors as in Kahan's report
00089 *     At this point, MX .NE. 0 so we can divide by it
00090 *
00091       SUN = SQRT( UNFL )
00092       SSUN = SQRT( SUN )
00093       SOV = SQRT( OVFL )
00094       TOM = SSUN*SOV
00095       IF( MX.LE.ONE ) THEN
00096          M1 = ONE / MX
00097          M2 = TOM
00098       ELSE
00099          M1 = ONE
00100          M2 = TOM / MX
00101       END IF
00102 *
00103 *     Begin counting
00104 *
00105       NUM = 0
00106       SSHIFT = ( SHIFT*M1 )*M2
00107       U = ( A( 1 )*M1 )*M2 - SSHIFT
00108       IF( U.LE.SUN ) THEN
00109          IF( U.LE.ZERO ) THEN
00110             NUM = NUM + 1
00111             IF( U.GT.-SUN )
00112      $         U = -SUN
00113          ELSE
00114             U = SUN
00115          END IF
00116       END IF
00117       DO 20 I = 2, N
00118          TMP = ( B( I-1 )*M1 )*M2
00119          U = ( ( A( I )*M1 )*M2-TMP*( TMP / U ) ) - SSHIFT
00120          IF( U.LE.SUN ) THEN
00121             IF( U.LE.ZERO ) THEN
00122                NUM = NUM + 1
00123                IF( U.GT.-SUN )
00124      $            U = -SUN
00125             ELSE
00126                U = SUN
00127             END IF
00128          END IF
00129    20 CONTINUE
00130       RETURN
00131 *
00132 *     End of DSTECT
00133 *
00134       END
 All Files Functions