LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ssvdct.f
Go to the documentation of this file.
00001 *> \brief \b SSVDCT
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 SSVDCT( N, S, E, SHIFT, NUM )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       INTEGER            N, NUM
00015 *       REAL               SHIFT
00016 *       ..
00017 *       .. Array Arguments ..
00018 *       REAL               E( * ), S( * )
00019 *       ..
00020 *  
00021 *
00022 *> \par Purpose:
00023 *  =============
00024 *>
00025 *> \verbatim
00026 *>
00027 *> SSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N
00028 *> tridiagonal matrix T which are less than or equal to SHIFT.  T is
00029 *> formed by putting zeros on the diagonal and making the off-diagonals
00030 *> equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N).  If SHIFT is
00031 *> positive, NUM is equal to N plus the number of singular values of a
00032 *> bidiagonal matrix B less than or equal to SHIFT.  Here B has diagonal
00033 *> entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1).
00034 *> If SHIFT is negative, NUM is equal to the number of singular values
00035 *> of B greater than or equal to -SHIFT.
00036 *>
00037 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00038 *> Matrix", Report CS41, Computer Science Dept., Stanford University,
00039 *> July 21, 1966
00040 *> \endverbatim
00041 *
00042 *  Arguments:
00043 *  ==========
00044 *
00045 *> \param[in] N
00046 *> \verbatim
00047 *>          N is INTEGER
00048 *>          The dimension of the bidiagonal matrix B.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] S
00052 *> \verbatim
00053 *>          S is REAL array, dimension (N)
00054 *>          The diagonal entries of the bidiagonal matrix B.
00055 *> \endverbatim
00056 *>
00057 *> \param[in] E
00058 *> \verbatim
00059 *>          E is REAL array of dimension (N-1)
00060 *>          The superdiagonal entries of the bidiagonal matrix B.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] SHIFT
00064 *> \verbatim
00065 *>          SHIFT is REAL
00066 *>          The shift, used as described under Purpose.
00067 *> \endverbatim
00068 *>
00069 *> \param[out] NUM
00070 *> \verbatim
00071 *>          NUM is INTEGER
00072 *>          The number of eigenvalues of T less than or equal to SHIFT.
00073 *> \endverbatim
00074 *
00075 *  Authors:
00076 *  ========
00077 *
00078 *> \author Univ. of Tennessee 
00079 *> \author Univ. of California Berkeley 
00080 *> \author Univ. of Colorado Denver 
00081 *> \author NAG Ltd. 
00082 *
00083 *> \date November 2011
00084 *
00085 *> \ingroup single_eig
00086 *
00087 *  =====================================================================
00088       SUBROUTINE SSVDCT( N, S, E, SHIFT, NUM )
00089 *
00090 *  -- LAPACK test routine (version 3.4.0) --
00091 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00092 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00093 *     November 2011
00094 *
00095 *     .. Scalar Arguments ..
00096       INTEGER            N, NUM
00097       REAL               SHIFT
00098 *     ..
00099 *     .. Array Arguments ..
00100       REAL               E( * ), S( * )
00101 *     ..
00102 *
00103 *  =====================================================================
00104 *
00105 *     .. Parameters ..
00106       REAL               ONE
00107       PARAMETER          ( ONE = 1.0E0 )
00108       REAL               ZERO
00109       PARAMETER          ( ZERO = 0.0E0 )
00110 *     ..
00111 *     .. Local Scalars ..
00112       INTEGER            I
00113       REAL               M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
00114      $                   TOM, U, UNFL
00115 *     ..
00116 *     .. External Functions ..
00117       REAL               SLAMCH
00118       EXTERNAL           SLAMCH
00119 *     ..
00120 *     .. Intrinsic Functions ..
00121       INTRINSIC          ABS, MAX, SQRT
00122 *     ..
00123 *     .. Executable Statements ..
00124 *
00125 *     Get machine constants
00126 *
00127       UNFL = 2*SLAMCH( 'Safe minimum' )
00128       OVFL = ONE / UNFL
00129 *
00130 *     Find largest entry
00131 *
00132       MX = ABS( S( 1 ) )
00133       DO 10 I = 1, N - 1
00134          MX = MAX( MX, ABS( S( I+1 ) ), ABS( E( I ) ) )
00135    10 CONTINUE
00136 *
00137       IF( MX.EQ.ZERO ) THEN
00138          IF( SHIFT.LT.ZERO ) THEN
00139             NUM = 0
00140          ELSE
00141             NUM = 2*N
00142          END IF
00143          RETURN
00144       END IF
00145 *
00146 *     Compute scale factors as in Kahan's report
00147 *
00148       SUN = SQRT( UNFL )
00149       SSUN = SQRT( SUN )
00150       SOV = SQRT( OVFL )
00151       TOM = SSUN*SOV
00152       IF( MX.LE.ONE ) THEN
00153          M1 = ONE / MX
00154          M2 = TOM
00155       ELSE
00156          M1 = ONE
00157          M2 = TOM / MX
00158       END IF
00159 *
00160 *     Begin counting
00161 *
00162       U = ONE
00163       NUM = 0
00164       SSHIFT = ( SHIFT*M1 )*M2
00165       U = -SSHIFT
00166       IF( U.LE.SUN ) THEN
00167          IF( U.LE.ZERO ) THEN
00168             NUM = NUM + 1
00169             IF( U.GT.-SUN )
00170      $         U = -SUN
00171          ELSE
00172             U = SUN
00173          END IF
00174       END IF
00175       TMP = ( S( 1 )*M1 )*M2
00176       U = -TMP*( TMP / U ) - SSHIFT
00177       IF( U.LE.SUN ) THEN
00178          IF( U.LE.ZERO ) THEN
00179             NUM = NUM + 1
00180             IF( U.GT.-SUN )
00181      $         U = -SUN
00182          ELSE
00183             U = SUN
00184          END IF
00185       END IF
00186       DO 20 I = 1, N - 1
00187          TMP = ( E( I )*M1 )*M2
00188          U = -TMP*( TMP / U ) - SSHIFT
00189          IF( U.LE.SUN ) THEN
00190             IF( U.LE.ZERO ) THEN
00191                NUM = NUM + 1
00192                IF( U.GT.-SUN )
00193      $            U = -SUN
00194             ELSE
00195                U = SUN
00196             END IF
00197          END IF
00198          TMP = ( S( I+1 )*M1 )*M2
00199          U = -TMP*( TMP / U ) - SSHIFT
00200          IF( U.LE.SUN ) THEN
00201             IF( U.LE.ZERO ) THEN
00202                NUM = NUM + 1
00203                IF( U.GT.-SUN )
00204      $            U = -SUN
00205             ELSE
00206                U = SUN
00207             END IF
00208          END IF
00209    20 CONTINUE
00210       RETURN
00211 *
00212 *     End of SSVDCT
00213 *
00214       END
 All Files Functions