LAPACK: TESTING/EIG/dsvdct.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSVDCT( N, S, E, 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   E( * ), S( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N
00019 *  tridiagonal matrix T which are less than or equal to SHIFT.  T is
00020 *  formed by putting zeros on the diagonal and making the off-diagonals
00021 *  equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N).  If SHIFT is
00022 *  positive, NUM is equal to N plus the number of singular values of a
00023 *  bidiagonal matrix B less than or equal to SHIFT.  Here B has diagonal
00024 *  entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1).
00025 *  If SHIFT is negative, NUM is equal to the number of singular values
00026 *  of B greater than or equal to -SHIFT.
00027 *
00028 *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00029 *  Matrix", Report CS41, Computer Science Dept., Stanford University,
00030 *  July 21, 1966
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  N       (input) INTEGER
00036 *          The dimension of the bidiagonal matrix B.
00037 *
00038 *  S       (input) DOUBLE PRECISION array, dimension (N)
00039 *          The diagonal entries of the bidiagonal matrix B.
00040 *
00041 *  E       (input) DOUBLE PRECISION array of dimension (N-1)
00042 *          The superdiagonal entries of the bidiagonal matrix B.
00043 *
00044 *  SHIFT   (input) DOUBLE PRECISION
00045 *          The shift, used as described under Purpose.
00046 *
00047 *  NUM     (output) INTEGER
00048 *          The number of eigenvalues of T less than or equal to SHIFT.
00049 *
00050 *  =====================================================================
00051 *
00052 *     .. Parameters ..
00053       DOUBLE PRECISION   ONE
00054       PARAMETER          ( ONE = 1.0D0 )
00055       DOUBLE PRECISION   ZERO
00056       PARAMETER          ( ZERO = 0.0D0 )
00057 *     ..
00058 *     .. Local Scalars ..
00059       INTEGER            I
00060       DOUBLE PRECISION   M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
00061      $                   TOM, U, UNFL
00062 *     ..
00063 *     .. External Functions ..
00064       DOUBLE PRECISION   DLAMCH
00065       EXTERNAL           DLAMCH
00066 *     ..
00067 *     .. Intrinsic Functions ..
00068       INTRINSIC          ABS, MAX, SQRT
00069 *     ..
00070 *     .. Executable Statements ..
00071 *
00072 *     Get machine constants
00073 *
00074       UNFL = 2*DLAMCH( 'Safe minimum' )
00075       OVFL = ONE / UNFL
00076 *
00077 *     Find largest entry
00078 *
00079       MX = ABS( S( 1 ) )
00080       DO 10 I = 1, N - 1
00081          MX = MAX( MX, ABS( S( I+1 ) ), ABS( E( I ) ) )
00082    10 CONTINUE
00083 *
00084       IF( MX.EQ.ZERO ) THEN
00085          IF( SHIFT.LT.ZERO ) THEN
00086             NUM = 0
00087          ELSE
00088             NUM = 2*N
00089          END IF
00090          RETURN
00091       END IF
00092 *
00093 *     Compute scale factors as in Kahan's report
00094 *
00095       SUN = SQRT( UNFL )
00096       SSUN = SQRT( SUN )
00097       SOV = SQRT( OVFL )
00098       TOM = SSUN*SOV
00099       IF( MX.LE.ONE ) THEN
00100          M1 = ONE / MX
00101          M2 = TOM
00102       ELSE
00103          M1 = ONE
00104          M2 = TOM / MX
00105       END IF
00106 *
00107 *     Begin counting
00108 *
00109       U = ONE
00110       NUM = 0
00111       SSHIFT = ( SHIFT*M1 )*M2
00112       U = -SSHIFT
00113       IF( U.LE.SUN ) THEN
00114          IF( U.LE.ZERO ) THEN
00115             NUM = NUM + 1
00116             IF( U.GT.-SUN )
00117      $         U = -SUN
00118          ELSE
00119             U = SUN
00120          END IF
00121       END IF
00122       TMP = ( S( 1 )*M1 )*M2
00123       U = -TMP*( TMP / U ) - SSHIFT
00124       IF( U.LE.SUN ) THEN
00125          IF( U.LE.ZERO ) THEN
00126             NUM = NUM + 1
00127             IF( U.GT.-SUN )
00128      $         U = -SUN
00129          ELSE
00130             U = SUN
00131          END IF
00132       END IF
00133       DO 20 I = 1, N - 1
00134          TMP = ( E( I )*M1 )*M2
00135          U = -TMP*( TMP / U ) - SSHIFT
00136          IF( U.LE.SUN ) THEN
00137             IF( U.LE.ZERO ) THEN
00138                NUM = NUM + 1
00139                IF( U.GT.-SUN )
00140      $            U = -SUN
00141             ELSE
00142                U = SUN
00143             END IF
00144          END IF
00145          TMP = ( S( I+1 )*M1 )*M2
00146          U = -TMP*( TMP / U ) - SSHIFT
00147          IF( U.LE.SUN ) THEN
00148             IF( U.LE.ZERO ) THEN
00149                NUM = NUM + 1
00150                IF( U.GT.-SUN )
00151      $            U = -SUN
00152             ELSE
00153                U = SUN
00154             END IF
00155          END IF
00156    20 CONTINUE
00157       RETURN
00158 *
00159 *     End of DSVDCT
00160 *
00161       END