LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dstech.f
Go to the documentation of this file.
1 *> \brief \b DSTECH
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE DSTECH( N, A, B, EIG, TOL, WORK, INFO )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER INFO, N
15 * DOUBLE PRECISION TOL
16 * ..
17 * .. Array Arguments ..
18 * DOUBLE PRECISION A( * ), B( * ), EIG( * ), WORK( * )
19 * ..
20 *
21 *
22 *> \par Purpose:
23 * =============
24 *>
25 *> \verbatim
26 *>
27 *> Let T be the tridiagonal matrix with diagonal entries A(1) ,...,
28 *> A(N) and offdiagonal entries B(1) ,..., B(N-1)). DSTECH checks to
29 *> see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T.
30 *> It does this by expanding each EIG(I) into an interval
31 *> [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if
32 *> any, and using Sturm sequences to count and verify whether each
33 *> resulting interval has the correct number of eigenvalues (using
34 *> DSTECT). Here EPS = TOL*MAZHEPS*MAXEIG, where MACHEPS is the
35 *> machine precision and MAXEIG is the absolute value of the largest
36 *> eigenvalue. If each interval contains the correct number of
37 *> eigenvalues, INFO = 0 is returned, otherwise INFO is the index of
38 *> the first eigenvalue in the first bad interval.
39 *> \endverbatim
40 *
41 * Arguments:
42 * ==========
43 *
44 *> \param[in] N
45 *> \verbatim
46 *> N is INTEGER
47 *> The dimension of the tridiagonal matrix T.
48 *> \endverbatim
49 *>
50 *> \param[in] A
51 *> \verbatim
52 *> A is DOUBLE PRECISION array, dimension (N)
53 *> The diagonal entries of the tridiagonal matrix T.
54 *> \endverbatim
55 *>
56 *> \param[in] B
57 *> \verbatim
58 *> B is DOUBLE PRECISION array, dimension (N-1)
59 *> The offdiagonal entries of the tridiagonal matrix T.
60 *> \endverbatim
61 *>
62 *> \param[in] EIG
63 *> \verbatim
64 *> EIG is DOUBLE PRECISION array, dimension (N)
65 *> The purported eigenvalues to be checked.
66 *> \endverbatim
67 *>
68 *> \param[in] TOL
69 *> \verbatim
70 *> TOL is DOUBLE PRECISION
71 *> Error tolerance for checking, a multiple of the
72 *> machine precision.
73 *> \endverbatim
74 *>
75 *> \param[out] WORK
76 *> \verbatim
77 *> WORK is DOUBLE PRECISION array, dimension (N)
78 *> \endverbatim
79 *>
80 *> \param[out] INFO
81 *> \verbatim
82 *> INFO is INTEGER
83 *> 0 if the eigenvalues are all correct (to within
84 *> 1 +- TOL*MAZHEPS*MAXEIG)
85 *> >0 if the interval containing the INFO-th eigenvalue
86 *> contains the incorrect number of eigenvalues.
87 *> \endverbatim
88 *
89 * Authors:
90 * ========
91 *
92 *> \author Univ. of Tennessee
93 *> \author Univ. of California Berkeley
94 *> \author Univ. of Colorado Denver
95 *> \author NAG Ltd.
96 *
97 *> \date November 2011
98 *
99 *> \ingroup double_eig
100 *
101 * =====================================================================
102  SUBROUTINE dstech( N, A, B, EIG, TOL, WORK, INFO )
103 *
104 * -- LAPACK test routine (version 3.4.0) --
105 * -- LAPACK is a software package provided by Univ. of Tennessee, --
106 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
107 * November 2011
108 *
109 * .. Scalar Arguments ..
110  INTEGER info, n
111  DOUBLE PRECISION tol
112 * ..
113 * .. Array Arguments ..
114  DOUBLE PRECISION a( * ), b( * ), eig( * ), work( * )
115 * ..
116 *
117 * =====================================================================
118 *
119 * .. Parameters ..
120  DOUBLE PRECISION zero
121  parameter( zero = 0.0d0 )
122 * ..
123 * .. Local Scalars ..
124  INTEGER bpnt, count, i, isub, j, numl, numu, tpnt
125  DOUBLE PRECISION emin, eps, lower, mx, tuppr, unflep, upper
126 * ..
127 * .. External Functions ..
128  DOUBLE PRECISION dlamch
129  EXTERNAL dlamch
130 * ..
131 * .. External Subroutines ..
132  EXTERNAL dstect
133 * ..
134 * .. Intrinsic Functions ..
135  INTRINSIC abs, max
136 * ..
137 * .. Executable Statements ..
138 *
139 * Check input parameters
140 *
141  info = 0
142  IF( n.EQ.0 )
143  $ RETURN
144  IF( n.LT.0 ) THEN
145  info = -1
146  RETURN
147  END IF
148  IF( tol.LT.zero ) THEN
149  info = -5
150  RETURN
151  END IF
152 *
153 * Get machine constants
154 *
155  eps = dlamch( 'Epsilon' )*dlamch( 'Base' )
156  unflep = dlamch( 'Safe minimum' ) / eps
157  eps = tol*eps
158 *
159 * Compute maximum absolute eigenvalue, error tolerance
160 *
161  mx = abs( eig( 1 ) )
162  DO 10 i = 2, n
163  mx = max( mx, abs( eig( i ) ) )
164  10 CONTINUE
165  eps = max( eps*mx, unflep )
166 *
167 * Sort eigenvalues from EIG into WORK
168 *
169  DO 20 i = 1, n
170  work( i ) = eig( i )
171  20 CONTINUE
172  DO 40 i = 1, n - 1
173  isub = 1
174  emin = work( 1 )
175  DO 30 j = 2, n + 1 - i
176  IF( work( j ).LT.emin ) THEN
177  isub = j
178  emin = work( j )
179  END IF
180  30 CONTINUE
181  IF( isub.NE.n+1-i ) THEN
182  work( isub ) = work( n+1-i )
183  work( n+1-i ) = emin
184  END IF
185  40 CONTINUE
186 *
187 * TPNT points to singular value at right endpoint of interval
188 * BPNT points to singular value at left endpoint of interval
189 *
190  tpnt = 1
191  bpnt = 1
192 *
193 * Begin loop over all intervals
194 *
195  50 CONTINUE
196  upper = work( tpnt ) + eps
197  lower = work( bpnt ) - eps
198 *
199 * Begin loop merging overlapping intervals
200 *
201  60 CONTINUE
202  IF( bpnt.EQ.n )
203  $ go to 70
204  tuppr = work( bpnt+1 ) + eps
205  IF( tuppr.LT.lower )
206  $ go to 70
207 *
208 * Merge
209 *
210  bpnt = bpnt + 1
211  lower = work( bpnt ) - eps
212  go to 60
213  70 CONTINUE
214 *
215 * Count singular values in interval [ LOWER, UPPER ]
216 *
217  CALL dstect( n, a, b, lower, numl )
218  CALL dstect( n, a, b, upper, numu )
219  count = numu - numl
220  IF( count.NE.bpnt-tpnt+1 ) THEN
221 *
222 * Wrong number of singular values in interval
223 *
224  info = tpnt
225  go to 80
226  END IF
227  tpnt = bpnt + 1
228  bpnt = tpnt
229  IF( tpnt.LE.n )
230  $ go to 50
231  80 CONTINUE
232  RETURN
233 *
234 * End of DSTECH
235 *
236  END