1*
2* $Id$
3*
4#include "blas_lapackf.h"
5*======================================================================
6*
7* DISCLAIMER
8*
9* This material was prepared as an account of work sponsored by an
10* agency of the United States Government.  Neither the United States
11* Government nor the United States Department of Energy, nor Battelle,
12* nor any of their employees, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
13* ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
14* COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT,
15* SOFTWARE, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT
16* INFRINGE PRIVATELY OWNED RIGHTS.
17*
18* ACKNOWLEDGMENT
19*
20* This software and its documentation were produced with Government
21* support under Contract Number DE-AC06-76RLO-1830 awarded by the United
22* States Department of Energy.  The Government retains a paid-up
23* non-exclusive, irrevocable worldwide license to reproduce, prepare
24* derivative works, perform publicly and display publicly by or for the
25* Government, including the right to distribute to other Government
26* contractors.
27*
28*======================================================================
29*
30*  -- PEIGS  routine (version 2.1) --
31*     Pacific Northwest Laboratory
32*     July 28, 1995
33*
34*======================================================================
35      INTEGER FUNCTION SNEBLW2(N, V, D, E, WORK, INFO )
36*
37*
38*     Same functionality as neblw1 from the Eispack library
39*     but using LAPACK code.
40*
41*     WARNING:  in this code E(1) is the first subdiagonal
42*     element of T and E(N) is junk.  While in
43*     neblw1 E(1) is junk and E(2) is the first
44*     subdiagonal element of T.
45*
46*     D.M. Elwood 09/26/94
47*
48*     .. Scalar Arguments ..
49      INTEGER            N, INFO
50      REAL   V
51*     ..
52*     .. Array Arguments ..
53      REAL   D( * ), E( * ), WORK( * )
54*     ..
55*
56*     Purpose
57*     =======
58*
59*     NEBLW2 computes the number of eigenvalues of a symmetric
60*     tridiagonal matrix, T, which are less than or equal to V.  Uses
61*     code from LAPACK's DSTEBZ and DLAEBZ.
62*
63*     Arguments
64*     =========
65*
66*     N       (input) INTEGER
67*     The dimension of the tridiagonal matrix T.  N >= 0.
68*
69*     V       (input) REAL
70*     Compute the number of eigenvalues less than or equal to V.
71*
72*     D       (input) REAL array, dimension (N)
73*          The n diagonal elements of the tridiagonal matrix T.  To
74*          avoid overflow, the matrix must be scaled so that its largest
75*          entry is no greater than overflow**(1/2) * underflow**(1/4)
76*          in absolute value, and for greatest accuracy, it should not
77*          be much smaller than that.
78*
79*  E       (input) REAL array, dimension (N-1)
80*          The (n-1) off-diagonal elements of the tridiagonal matrix T.
81*          To avoid overflow, the matrix must be scaled so that its
82*          largest entry is no greater than overflow**(1/2) *
83*          underflow**(1/4) in absolute value, and for greatest
84*          accuracy, it should not be much smaller than that.
85*
86*  WORK    (workspace) REAL array, dimension (N)
87*
88*  INFO    (output) INTEGER
89*          = 0:  successful exit
90*          < 0:  if INFO = -i, the i-th argument had an illegal value
91*  =====================================================================
92*
93*     .. Parameters ..
94      REAL   ZERO, ONE
95      PARAMETER          ( ZERO = 0.0, ONE = 1.0 )
96*     ..
97*     .. Local Scalars ..
98      INTEGER            J, MOUT
99      REAL   PIVMIN, SAFEMN, TMP1, ULP
100*     ..
101*     .. External Functions ..
102*     ..
103*     .. External Subroutines ..
104      EXTERNAL           XERBLA
105*     ..
106*     .. Intrinsic Functions ..
107      INTRINSIC          ABS, MAX
108#include "blas_lapack.data"
109*     ..
110*     .. Executable Statements ..
111*
112      INFO = 0
113*
114*     Check for Errors
115*
116      IF( N.LT.0 )
117     $   INFO = -1
118*
119      IF( INFO.NE.0 ) THEN
120         CALL XERBLA( 'NEWBLW2', -INFO )
121         RETURN
122      END IF
123*
124*     Quick return if possible
125*
126      SNEBLW2 = 0
127      IF( N.EQ.0 )
128     $   RETURN
129*
130*     Get machine constants
131*
132c      SAFEMN = SLAMCH( 'S' )
133c      ULP = SLAMCH( 'P' )
134c
135      SAFEMN = SLAMCHS
136      ULP = SLAMCHP
137
138*
139*     Special Case when N=1
140*
141      IF( N.EQ.1 ) THEN
142         IF( D( 1 ) .LE. V ) THEN
143            SNEBLW2 = 1
144         ELSE
145            SNEBLW2 = 0
146         END IF
147         RETURN
148      END IF
149*
150*     Compute Splitting Points
151*     Set WORK = E * E
152*
153      WORK( N ) = ZERO
154      PIVMIN = ONE
155*
156      DO 10 J = 2, N
157         TMP1 = E( J-1 )**2
158         IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
159            WORK( J-1 ) = ZERO
160         ELSE
161            WORK( J-1 ) = TMP1
162            PIVMIN = MAX( PIVMIN, TMP1 )
163         END IF
164   10 CONTINUE
165      PIVMIN = PIVMIN*SAFEMN
166*
167*     Compute the number of eigenvalues in the initial intervals.
168*
169      TMP1 = D( 1 ) - V
170      IF( ABS( TMP1 ).LT.PIVMIN )
171     $   TMP1 = -PIVMIN
172      MOUT = 0
173      IF( TMP1.LE.ZERO )
174     $   MOUT = 1
175*
176      DO 20 J = 2, N
177         TMP1 = D( J ) - WORK( J-1 ) / TMP1 - V
178         IF( ABS( TMP1 ).LT.PIVMIN )
179     $      TMP1 = -PIVMIN
180         IF( TMP1.LE.ZERO )
181     $      MOUT = MOUT + 1
182   20 CONTINUE
183      SNEBLW2 = MOUT
184*
185*     End of NEWBLW2
186*
187      RETURN
188      END
189