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 NEBLW2(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      DOUBLE PRECISION   V
51*     ..
52*     .. Array Arguments ..
53      DOUBLE PRECISION   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) DOUBLE PRECISION
70*          Compute the number of eigenvalues less than or equal to V.
71*
72*  D       (input) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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      DOUBLE PRECISION   ZERO, ONE
95      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
96*     ..
97*     .. Local Scalars ..
98      INTEGER            J, MOUT
99      DOUBLE PRECISION   PIVMIN, SAFEMN, TMP1, ULP
100*     ..
101*     .. External Functions ..
102      DOUBLE PRECISION   DLAMCH
103      EXTERNAL           DLAMCH
104*     ..
105*     .. External Subroutines ..
106      EXTERNAL           XERBLA
107*     ..
108*     .. Intrinsic Functions ..
109      INTRINSIC          ABS, MAX
110#include "blas_lapack.data"
111*     ..
112*     .. Executable Statements ..
113*
114      INFO = 0
115*
116*     Check for Errors
117*
118      IF( N.LT.0 )
119     $   INFO = -1
120*
121      IF( INFO.NE.0 ) THEN
122         CALL XERBLA( 'NEWBLW2', -INFO )
123         RETURN
124      END IF
125*
126*     Quick return if possible
127*
128      NEBLW2 = 0
129      IF( N.EQ.0 )
130     $   RETURN
131*
132*     Get machine constants
133*
134c      SAFEMN = DLAMCH( 'S' )
135c      ULP = DLAMCH( 'P' )
136       SAFEMN = DLAMCHS
137       ULP = DLAMCHP
138*
139*     Special Case when N=1
140*
141      IF( N.EQ.1 ) THEN
142         IF( D( 1 ) .LE. V ) THEN
143            NEBLW2 = 1
144         ELSE
145            NEBLW2 = 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      NEBLW2 = MOUT
184*
185*     End of NEWBLW2
186*
187      RETURN
188      END
189