LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlarrv.f
Go to the documentation of this file.
1 *> \brief \b DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLARRV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
22 * ISPLIT, M, DOL, DOU, MINRGP,
23 * RTOL1, RTOL2, W, WERR, WGAP,
24 * IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
25 * WORK, IWORK, INFO )
26 *
27 * .. Scalar Arguments ..
28 * INTEGER DOL, DOU, INFO, LDZ, M, N
29 * DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
33 * $ ISUPPZ( * ), IWORK( * )
34 * DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
35 * $ WGAP( * ), WORK( * )
36 * DOUBLE PRECISION Z( LDZ, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> DLARRV computes the eigenvectors of the tridiagonal matrix
46 *> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
47 *> The input eigenvalues should have been computed by DLARRE.
48 *> \endverbatim
49 *
50 * Arguments:
51 * ==========
52 *
53 *> \param[in] N
54 *> \verbatim
55 *> N is INTEGER
56 *> The order of the matrix. N >= 0.
57 *> \endverbatim
58 *>
59 *> \param[in] VL
60 *> \verbatim
61 *> VL is DOUBLE PRECISION
62 *> \endverbatim
63 *>
64 *> \param[in] VU
65 *> \verbatim
66 *> VU is DOUBLE PRECISION
67 *> Lower and upper bounds of the interval that contains the desired
68 *> eigenvalues. VL < VU. Needed to compute gaps on the left or right
69 *> end of the extremal eigenvalues in the desired RANGE.
70 *> \endverbatim
71 *>
72 *> \param[in,out] D
73 *> \verbatim
74 *> D is DOUBLE PRECISION array, dimension (N)
75 *> On entry, the N diagonal elements of the diagonal matrix D.
76 *> On exit, D may be overwritten.
77 *> \endverbatim
78 *>
79 *> \param[in,out] L
80 *> \verbatim
81 *> L is DOUBLE PRECISION array, dimension (N)
82 *> On entry, the (N-1) subdiagonal elements of the unit
83 *> bidiagonal matrix L are in elements 1 to N-1 of L
84 *> (if the matrix is not splitted.) At the end of each block
85 *> is stored the corresponding shift as given by DLARRE.
86 *> On exit, L is overwritten.
87 *> \endverbatim
88 *>
89 *> \param[in] PIVMIN
90 *> \verbatim
91 *> PIVMIN is DOUBLE PRECISION
92 *> The minimum pivot allowed in the Sturm sequence.
93 *> \endverbatim
94 *>
95 *> \param[in] ISPLIT
96 *> \verbatim
97 *> ISPLIT is INTEGER array, dimension (N)
98 *> The splitting points, at which T breaks up into blocks.
99 *> The first block consists of rows/columns 1 to
100 *> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
101 *> through ISPLIT( 2 ), etc.
102 *> \endverbatim
103 *>
104 *> \param[in] M
105 *> \verbatim
106 *> M is INTEGER
107 *> The total number of input eigenvalues. 0 <= M <= N.
108 *> \endverbatim
109 *>
110 *> \param[in] DOL
111 *> \verbatim
112 *> DOL is INTEGER
113 *> \endverbatim
114 *>
115 *> \param[in] DOU
116 *> \verbatim
117 *> DOU is INTEGER
118 *> If the user wants to compute only selected eigenvectors from all
119 *> the eigenvalues supplied, he can specify an index range DOL:DOU.
120 *> Or else the setting DOL=1, DOU=M should be applied.
121 *> Note that DOL and DOU refer to the order in which the eigenvalues
122 *> are stored in W.
123 *> If the user wants to compute only selected eigenpairs, then
124 *> the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
125 *> computed eigenvectors. All other columns of Z are set to zero.
126 *> \endverbatim
127 *>
128 *> \param[in] MINRGP
129 *> \verbatim
130 *> MINRGP is DOUBLE PRECISION
131 *> \endverbatim
132 *>
133 *> \param[in] RTOL1
134 *> \verbatim
135 *> RTOL1 is DOUBLE PRECISION
136 *> \endverbatim
137 *>
138 *> \param[in] RTOL2
139 *> \verbatim
140 *> RTOL2 is DOUBLE PRECISION
141 *> Parameters for bisection.
142 *> An interval [LEFT,RIGHT] has converged if
143 *> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
144 *> \endverbatim
145 *>
146 *> \param[in,out] W
147 *> \verbatim
148 *> W is DOUBLE PRECISION array, dimension (N)
149 *> The first M elements of W contain the APPROXIMATE eigenvalues for
150 *> which eigenvectors are to be computed. The eigenvalues
151 *> should be grouped by split-off block and ordered from
152 *> smallest to largest within the block ( The output array
153 *> W from DLARRE is expected here ). Furthermore, they are with
154 *> respect to the shift of the corresponding root representation
155 *> for their block. On exit, W holds the eigenvalues of the
156 *> UNshifted matrix.
157 *> \endverbatim
158 *>
159 *> \param[in,out] WERR
160 *> \verbatim
161 *> WERR is DOUBLE PRECISION array, dimension (N)
162 *> The first M elements contain the semiwidth of the uncertainty
163 *> interval of the corresponding eigenvalue in W
164 *> \endverbatim
165 *>
166 *> \param[in,out] WGAP
167 *> \verbatim
168 *> WGAP is DOUBLE PRECISION array, dimension (N)
169 *> The separation from the right neighbor eigenvalue in W.
170 *> \endverbatim
171 *>
172 *> \param[in] IBLOCK
173 *> \verbatim
174 *> IBLOCK is INTEGER array, dimension (N)
175 *> The indices of the blocks (submatrices) associated with the
176 *> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
177 *> W(i) belongs to the first block from the top, =2 if W(i)
178 *> belongs to the second block, etc.
179 *> \endverbatim
180 *>
181 *> \param[in] INDEXW
182 *> \verbatim
183 *> INDEXW is INTEGER array, dimension (N)
184 *> The indices of the eigenvalues within each block (submatrix);
185 *> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
186 *> i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
187 *> \endverbatim
188 *>
189 *> \param[in] GERS
190 *> \verbatim
191 *> GERS is DOUBLE PRECISION array, dimension (2*N)
192 *> The N Gerschgorin intervals (the i-th Gerschgorin interval
193 *> is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
194 *> be computed from the original UNshifted matrix.
195 *> \endverbatim
196 *>
197 *> \param[out] Z
198 *> \verbatim
199 *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
200 *> If INFO = 0, the first M columns of Z contain the
201 *> orthonormal eigenvectors of the matrix T
202 *> corresponding to the input eigenvalues, with the i-th
203 *> column of Z holding the eigenvector associated with W(i).
204 *> Note: the user must ensure that at least max(1,M) columns are
205 *> supplied in the array Z.
206 *> \endverbatim
207 *>
208 *> \param[in] LDZ
209 *> \verbatim
210 *> LDZ is INTEGER
211 *> The leading dimension of the array Z. LDZ >= 1, and if
212 *> JOBZ = 'V', LDZ >= max(1,N).
213 *> \endverbatim
214 *>
215 *> \param[out] ISUPPZ
216 *> \verbatim
217 *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
218 *> The support of the eigenvectors in Z, i.e., the indices
219 *> indicating the nonzero elements in Z. The I-th eigenvector
220 *> is nonzero only in elements ISUPPZ( 2*I-1 ) through
221 *> ISUPPZ( 2*I ).
222 *> \endverbatim
223 *>
224 *> \param[out] WORK
225 *> \verbatim
226 *> WORK is DOUBLE PRECISION array, dimension (12*N)
227 *> \endverbatim
228 *>
229 *> \param[out] IWORK
230 *> \verbatim
231 *> IWORK is INTEGER array, dimension (7*N)
232 *> \endverbatim
233 *>
234 *> \param[out] INFO
235 *> \verbatim
236 *> INFO is INTEGER
237 *> = 0: successful exit
238 *>
239 *> > 0: A problem occured in DLARRV.
240 *> < 0: One of the called subroutines signaled an internal problem.
241 *> Needs inspection of the corresponding parameter IINFO
242 *> for further information.
243 *>
244 *> =-1: Problem in DLARRB when refining a child's eigenvalues.
245 *> =-2: Problem in DLARRF when computing the RRR of a child.
246 *> When a child is inside a tight cluster, it can be difficult
247 *> to find an RRR. A partial remedy from the user's point of
248 *> view is to make the parameter MINRGP smaller and recompile.
249 *> However, as the orthogonality of the computed vectors is
250 *> proportional to 1/MINRGP, the user should be aware that
251 *> he might be trading in precision when he decreases MINRGP.
252 *> =-3: Problem in DLARRB when refining a single eigenvalue
253 *> after the Rayleigh correction was rejected.
254 *> = 5: The Rayleigh Quotient Iteration failed to converge to
255 *> full accuracy in MAXITR steps.
256 *> \endverbatim
257 *
258 * Authors:
259 * ========
260 *
261 *> \author Univ. of Tennessee
262 *> \author Univ. of California Berkeley
263 *> \author Univ. of Colorado Denver
264 *> \author NAG Ltd.
265 *
266 *> \date September 2012
267 *
268 *> \ingroup doubleOTHERauxiliary
269 *
270 *> \par Contributors:
271 * ==================
272 *>
273 *> Beresford Parlett, University of California, Berkeley, USA \n
274 *> Jim Demmel, University of California, Berkeley, USA \n
275 *> Inderjit Dhillon, University of Texas, Austin, USA \n
276 *> Osni Marques, LBNL/NERSC, USA \n
277 *> Christof Voemel, University of California, Berkeley, USA
278 *
279 * =====================================================================
280  SUBROUTINE dlarrv( N, VL, VU, D, L, PIVMIN,
281  $ isplit, m, dol, dou, minrgp,
282  $ rtol1, rtol2, w, werr, wgap,
283  $ iblock, indexw, gers, z, ldz, isuppz,
284  $ work, iwork, info )
285 *
286 * -- LAPACK auxiliary routine (version 3.4.2) --
287 * -- LAPACK is a software package provided by Univ. of Tennessee, --
288 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
289 * September 2012
290 *
291 * .. Scalar Arguments ..
292  INTEGER dol, dou, info, ldz, m, n
293  DOUBLE PRECISION minrgp, pivmin, rtol1, rtol2, vl, vu
294 * ..
295 * .. Array Arguments ..
296  INTEGER iblock( * ), indexw( * ), isplit( * ),
297  $ isuppz( * ), iwork( * )
298  DOUBLE PRECISION d( * ), gers( * ), l( * ), w( * ), werr( * ),
299  $ wgap( * ), work( * )
300  DOUBLE PRECISION z( ldz, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  INTEGER maxitr
307  parameter( maxitr = 10 )
308  DOUBLE PRECISION zero, one, two, three, four, half
309  parameter( zero = 0.0d0, one = 1.0d0,
310  $ two = 2.0d0, three = 3.0d0,
311  $ four = 4.0d0, half = 0.5d0)
312 * ..
313 * .. Local Scalars ..
314  LOGICAL eskip, needbs, stp2ii, tryrqc, usedbs, usedrq
315  INTEGER done, i, ibegin, idone, iend, ii, iindc1,
316  $ iindc2, iindr, iindwk, iinfo, im, in, indeig,
317  $ indld, indlld, indwrk, isupmn, isupmx, iter,
318  $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
319  $ ndepth, negcnt, newcls, newfst, newftt, newlst,
320  $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
321  $ oldncl, p, parity, q, wbegin, wend, windex,
322  $ windmn, windpl, zfrom, zto, zusedl, zusedu,
323  $ zusedw
324  DOUBLE PRECISION bstres, bstw, eps, fudge, gap, gaptol, gl, gu,
325  $ lambda, left, lgap, mingma, nrminv, resid,
326  $ rgap, right, rqcorr, rqtol, savgap, sgndef,
327  $ sigma, spdiam, ssigma, tau, tmp, tol, ztz
328 * ..
329 * .. External Functions ..
330  DOUBLE PRECISION dlamch
331  EXTERNAL dlamch
332 * ..
333 * .. External Subroutines ..
334  EXTERNAL dcopy, dlar1v, dlarrb, dlarrf, dlaset,
335  $ dscal
336 * ..
337 * .. Intrinsic Functions ..
338  INTRINSIC abs, dble, max, min
339 * ..
340 * .. Executable Statements ..
341 * ..
342 
343 * The first N entries of WORK are reserved for the eigenvalues
344  indld = n+1
345  indlld= 2*n+1
346  indwrk= 3*n+1
347  minwsize = 12 * n
348 
349  DO 5 i= 1,minwsize
350  work( i ) = zero
351  5 continue
352 
353 * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
354 * factorization used to compute the FP vector
355  iindr = 0
356 * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
357 * layer and the one above.
358  iindc1 = n
359  iindc2 = 2*n
360  iindwk = 3*n + 1
361 
362  miniwsize = 7 * n
363  DO 10 i= 1,miniwsize
364  iwork( i ) = 0
365  10 continue
366 
367  zusedl = 1
368  IF(dol.GT.1) THEN
369 * Set lower bound for use of Z
370  zusedl = dol-1
371  ENDIF
372  zusedu = m
373  IF(dou.LT.m) THEN
374 * Set lower bound for use of Z
375  zusedu = dou+1
376  ENDIF
377 * The width of the part of Z that is used
378  zusedw = zusedu - zusedl + 1
379 
380 
381  CALL dlaset( 'Full', n, zusedw, zero, zero,
382  $ z(1,zusedl), ldz )
383 
384  eps = dlamch( 'Precision' )
385  rqtol = two * eps
386 *
387 * Set expert flags for standard code.
388  tryrqc = .true.
389 
390  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
391  ELSE
392 * Only selected eigenpairs are computed. Since the other evalues
393 * are not refined by RQ iteration, bisection has to compute to full
394 * accuracy.
395  rtol1 = four * eps
396  rtol2 = four * eps
397  ENDIF
398 
399 * The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
400 * desired eigenvalues. The support of the nonzero eigenvector
401 * entries is contained in the interval IBEGIN:IEND.
402 * Remark that if k eigenpairs are desired, then the eigenvectors
403 * are stored in k contiguous columns of Z.
404 
405 * DONE is the number of eigenvectors already computed
406  done = 0
407  ibegin = 1
408  wbegin = 1
409  DO 170 jblk = 1, iblock( m )
410  iend = isplit( jblk )
411  sigma = l( iend )
412 * Find the eigenvectors of the submatrix indexed IBEGIN
413 * through IEND.
414  wend = wbegin - 1
415  15 continue
416  IF( wend.LT.m ) THEN
417  IF( iblock( wend+1 ).EQ.jblk ) THEN
418  wend = wend + 1
419  go to 15
420  END IF
421  END IF
422  IF( wend.LT.wbegin ) THEN
423  ibegin = iend + 1
424  go to 170
425  elseif( (wend.LT.dol).OR.(wbegin.GT.dou) ) THEN
426  ibegin = iend + 1
427  wbegin = wend + 1
428  go to 170
429  END IF
430 
431 * Find local spectral diameter of the block
432  gl = gers( 2*ibegin-1 )
433  gu = gers( 2*ibegin )
434  DO 20 i = ibegin+1 , iend
435  gl = min( gers( 2*i-1 ), gl )
436  gu = max( gers( 2*i ), gu )
437  20 continue
438  spdiam = gu - gl
439 
440 * OLDIEN is the last index of the previous block
441  oldien = ibegin - 1
442 * Calculate the size of the current block
443  in = iend - ibegin + 1
444 * The number of eigenvalues in the current block
445  im = wend - wbegin + 1
446 
447 * This is for a 1x1 block
448  IF( ibegin.EQ.iend ) THEN
449  done = done+1
450  z( ibegin, wbegin ) = one
451  isuppz( 2*wbegin-1 ) = ibegin
452  isuppz( 2*wbegin ) = ibegin
453  w( wbegin ) = w( wbegin ) + sigma
454  work( wbegin ) = w( wbegin )
455  ibegin = iend + 1
456  wbegin = wbegin + 1
457  go to 170
458  END IF
459 
460 * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
461 * Note that these can be approximations, in this case, the corresp.
462 * entries of WERR give the size of the uncertainty interval.
463 * The eigenvalue approximations will be refined when necessary as
464 * high relative accuracy is required for the computation of the
465 * corresponding eigenvectors.
466  CALL dcopy( im, w( wbegin ), 1,
467  $ work( wbegin ), 1 )
468 
469 * We store in W the eigenvalue approximations w.r.t. the original
470 * matrix T.
471  DO 30 i=1,im
472  w(wbegin+i-1) = w(wbegin+i-1)+sigma
473  30 continue
474 
475 
476 * NDEPTH is the current depth of the representation tree
477  ndepth = 0
478 * PARITY is either 1 or 0
479  parity = 1
480 * NCLUS is the number of clusters for the next level of the
481 * representation tree, we start with NCLUS = 1 for the root
482  nclus = 1
483  iwork( iindc1+1 ) = 1
484  iwork( iindc1+2 ) = im
485 
486 * IDONE is the number of eigenvectors already computed in the current
487 * block
488  idone = 0
489 * loop while( IDONE.LT.IM )
490 * generate the representation tree for the current block and
491 * compute the eigenvectors
492  40 continue
493  IF( idone.LT.im ) THEN
494 * This is a crude protection against infinitely deep trees
495  IF( ndepth.GT.m ) THEN
496  info = -2
497  return
498  ENDIF
499 * breadth first processing of the current level of the representation
500 * tree: OLDNCL = number of clusters on current level
501  oldncl = nclus
502 * reset NCLUS to count the number of child clusters
503  nclus = 0
504 *
505  parity = 1 - parity
506  IF( parity.EQ.0 ) THEN
507  oldcls = iindc1
508  newcls = iindc2
509  ELSE
510  oldcls = iindc2
511  newcls = iindc1
512  END IF
513 * Process the clusters on the current level
514  DO 150 i = 1, oldncl
515  j = oldcls + 2*i
516 * OLDFST, OLDLST = first, last index of current cluster.
517 * cluster indices start with 1 and are relative
518 * to WBEGIN when accessing W, WGAP, WERR, Z
519  oldfst = iwork( j-1 )
520  oldlst = iwork( j )
521  IF( ndepth.GT.0 ) THEN
522 * Retrieve relatively robust representation (RRR) of cluster
523 * that has been computed at the previous level
524 * The RRR is stored in Z and overwritten once the eigenvectors
525 * have been computed or when the cluster is refined
526 
527  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
528 * Get representation from location of the leftmost evalue
529 * of the cluster
530  j = wbegin + oldfst - 1
531  ELSE
532  IF(wbegin+oldfst-1.LT.dol) THEN
533 * Get representation from the left end of Z array
534  j = dol - 1
535  elseif(wbegin+oldfst-1.GT.dou) THEN
536 * Get representation from the right end of Z array
537  j = dou
538  ELSE
539  j = wbegin + oldfst - 1
540  ENDIF
541  ENDIF
542  CALL dcopy( in, z( ibegin, j ), 1, d( ibegin ), 1 )
543  CALL dcopy( in-1, z( ibegin, j+1 ), 1, l( ibegin ),
544  $ 1 )
545  sigma = z( iend, j+1 )
546 
547 * Set the corresponding entries in Z to zero
548  CALL dlaset( 'Full', in, 2, zero, zero,
549  $ z( ibegin, j), ldz )
550  END IF
551 
552 * Compute DL and DLL of current RRR
553  DO 50 j = ibegin, iend-1
554  tmp = d( j )*l( j )
555  work( indld-1+j ) = tmp
556  work( indlld-1+j ) = tmp*l( j )
557  50 continue
558 
559  IF( ndepth.GT.0 ) THEN
560 * P and Q are index of the first and last eigenvalue to compute
561 * within the current block
562  p = indexw( wbegin-1+oldfst )
563  q = indexw( wbegin-1+oldlst )
564 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
565 * through the Q-OFFSET elements of these arrays are to be used.
566 * OFFSET = P-OLDFST
567  offset = indexw( wbegin ) - 1
568 * perform limited bisection (if necessary) to get approximate
569 * eigenvalues to the precision needed.
570  CALL dlarrb( in, d( ibegin ),
571  $ work(indlld+ibegin-1),
572  $ p, q, rtol1, rtol2, offset,
573  $ work(wbegin),wgap(wbegin),werr(wbegin),
574  $ work( indwrk ), iwork( iindwk ),
575  $ pivmin, spdiam, in, iinfo )
576  IF( iinfo.NE.0 ) THEN
577  info = -1
578  return
579  ENDIF
580 * We also recompute the extremal gaps. W holds all eigenvalues
581 * of the unshifted matrix and must be used for computation
582 * of WGAP, the entries of WORK might stem from RRRs with
583 * different shifts. The gaps from WBEGIN-1+OLDFST to
584 * WBEGIN-1+OLDLST are correctly computed in DLARRB.
585 * However, we only allow the gaps to become greater since
586 * this is what should happen when we decrease WERR
587  IF( oldfst.GT.1) THEN
588  wgap( wbegin+oldfst-2 ) =
589  $ max(wgap(wbegin+oldfst-2),
590  $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
591  $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
592  ENDIF
593  IF( wbegin + oldlst -1 .LT. wend ) THEN
594  wgap( wbegin+oldlst-1 ) =
595  $ max(wgap(wbegin+oldlst-1),
596  $ w(wbegin+oldlst)-werr(wbegin+oldlst)
597  $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
598  ENDIF
599 * Each time the eigenvalues in WORK get refined, we store
600 * the newly found approximation with all shifts applied in W
601  DO 53 j=oldfst,oldlst
602  w(wbegin+j-1) = work(wbegin+j-1)+sigma
603  53 continue
604  END IF
605 
606 * Process the current node.
607  newfst = oldfst
608  DO 140 j = oldfst, oldlst
609  IF( j.EQ.oldlst ) THEN
610 * we are at the right end of the cluster, this is also the
611 * boundary of the child cluster
612  newlst = j
613  ELSE IF ( wgap( wbegin + j -1).GE.
614  $ minrgp* abs( work(wbegin + j -1) ) ) THEN
615 * the right relative gap is big enough, the child cluster
616 * (NEWFST,..,NEWLST) is well separated from the following
617  newlst = j
618  ELSE
619 * inside a child cluster, the relative gap is not
620 * big enough.
621  goto 140
622  END IF
623 
624 * Compute size of child cluster found
625  newsiz = newlst - newfst + 1
626 
627 * NEWFTT is the place in Z where the new RRR or the computed
628 * eigenvector is to be stored
629  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
630 * Store representation at location of the leftmost evalue
631 * of the cluster
632  newftt = wbegin + newfst - 1
633  ELSE
634  IF(wbegin+newfst-1.LT.dol) THEN
635 * Store representation at the left end of Z array
636  newftt = dol - 1
637  elseif(wbegin+newfst-1.GT.dou) THEN
638 * Store representation at the right end of Z array
639  newftt = dou
640  ELSE
641  newftt = wbegin + newfst - 1
642  ENDIF
643  ENDIF
644 
645  IF( newsiz.GT.1) THEN
646 *
647 * Current child is not a singleton but a cluster.
648 * Compute and store new representation of child.
649 *
650 *
651 * Compute left and right cluster gap.
652 *
653 * LGAP and RGAP are not computed from WORK because
654 * the eigenvalue approximations may stem from RRRs
655 * different shifts. However, W hold all eigenvalues
656 * of the unshifted matrix. Still, the entries in WGAP
657 * have to be computed from WORK since the entries
658 * in W might be of the same order so that gaps are not
659 * exhibited correctly for very close eigenvalues.
660  IF( newfst.EQ.1 ) THEN
661  lgap = max( zero,
662  $ w(wbegin)-werr(wbegin) - vl )
663  ELSE
664  lgap = wgap( wbegin+newfst-2 )
665  ENDIF
666  rgap = wgap( wbegin+newlst-1 )
667 *
668 * Compute left- and rightmost eigenvalue of child
669 * to high precision in order to shift as close
670 * as possible and obtain as large relative gaps
671 * as possible
672 *
673  DO 55 k =1,2
674  IF(k.EQ.1) THEN
675  p = indexw( wbegin-1+newfst )
676  ELSE
677  p = indexw( wbegin-1+newlst )
678  ENDIF
679  offset = indexw( wbegin ) - 1
680  CALL dlarrb( in, d(ibegin),
681  $ work( indlld+ibegin-1 ),p,p,
682  $ rqtol, rqtol, offset,
683  $ work(wbegin),wgap(wbegin),
684  $ werr(wbegin),work( indwrk ),
685  $ iwork( iindwk ), pivmin, spdiam,
686  $ in, iinfo )
687  55 continue
688 *
689  IF((wbegin+newlst-1.LT.dol).OR.
690  $ (wbegin+newfst-1.GT.dou)) THEN
691 * if the cluster contains no desired eigenvalues
692 * skip the computation of that branch of the rep. tree
693 *
694 * We could skip before the refinement of the extremal
695 * eigenvalues of the child, but then the representation
696 * tree could be different from the one when nothing is
697 * skipped. For this reason we skip at this place.
698  idone = idone + newlst - newfst + 1
699  goto 139
700  ENDIF
701 *
702 * Compute RRR of child cluster.
703 * Note that the new RRR is stored in Z
704 *
705 * DLARRF needs LWORK = 2*N
706  CALL dlarrf( in, d( ibegin ), l( ibegin ),
707  $ work(indld+ibegin-1),
708  $ newfst, newlst, work(wbegin),
709  $ wgap(wbegin), werr(wbegin),
710  $ spdiam, lgap, rgap, pivmin, tau,
711  $ z(ibegin, newftt),z(ibegin, newftt+1),
712  $ work( indwrk ), iinfo )
713  IF( iinfo.EQ.0 ) THEN
714 * a new RRR for the cluster was found by DLARRF
715 * update shift and store it
716  ssigma = sigma + tau
717  z( iend, newftt+1 ) = ssigma
718 * WORK() are the midpoints and WERR() the semi-width
719 * Note that the entries in W are unchanged.
720  DO 116 k = newfst, newlst
721  fudge =
722  $ three*eps*abs(work(wbegin+k-1))
723  work( wbegin + k - 1 ) =
724  $ work( wbegin + k - 1) - tau
725  fudge = fudge +
726  $ four*eps*abs(work(wbegin+k-1))
727 * Fudge errors
728  werr( wbegin + k - 1 ) =
729  $ werr( wbegin + k - 1 ) + fudge
730 * Gaps are not fudged. Provided that WERR is small
731 * when eigenvalues are close, a zero gap indicates
732 * that a new representation is needed for resolving
733 * the cluster. A fudge could lead to a wrong decision
734 * of judging eigenvalues 'separated' which in
735 * reality are not. This could have a negative impact
736 * on the orthogonality of the computed eigenvectors.
737  116 continue
738 
739  nclus = nclus + 1
740  k = newcls + 2*nclus
741  iwork( k-1 ) = newfst
742  iwork( k ) = newlst
743  ELSE
744  info = -2
745  return
746  ENDIF
747  ELSE
748 *
749 * Compute eigenvector of singleton
750 *
751  iter = 0
752 *
753  tol = four * log(dble(in)) * eps
754 *
755  k = newfst
756  windex = wbegin + k - 1
757  windmn = max(windex - 1,1)
758  windpl = min(windex + 1,m)
759  lambda = work( windex )
760  done = done + 1
761 * Check if eigenvector computation is to be skipped
762  IF((windex.LT.dol).OR.
763  $ (windex.GT.dou)) THEN
764  eskip = .true.
765  goto 125
766  ELSE
767  eskip = .false.
768  ENDIF
769  left = work( windex ) - werr( windex )
770  right = work( windex ) + werr( windex )
771  indeig = indexw( windex )
772 * Note that since we compute the eigenpairs for a child,
773 * all eigenvalue approximations are w.r.t the same shift.
774 * In this case, the entries in WORK should be used for
775 * computing the gaps since they exhibit even very small
776 * differences in the eigenvalues, as opposed to the
777 * entries in W which might "look" the same.
778 
779  IF( k .EQ. 1) THEN
780 * In the case RANGE='I' and with not much initial
781 * accuracy in LAMBDA and VL, the formula
782 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
783 * can lead to an overestimation of the left gap and
784 * thus to inadequately early RQI 'convergence'.
785 * Prevent this by forcing a small left gap.
786  lgap = eps*max(abs(left),abs(right))
787  ELSE
788  lgap = wgap(windmn)
789  ENDIF
790  IF( k .EQ. im) THEN
791 * In the case RANGE='I' and with not much initial
792 * accuracy in LAMBDA and VU, the formula
793 * can lead to an overestimation of the right gap and
794 * thus to inadequately early RQI 'convergence'.
795 * Prevent this by forcing a small right gap.
796  rgap = eps*max(abs(left),abs(right))
797  ELSE
798  rgap = wgap(windex)
799  ENDIF
800  gap = min( lgap, rgap )
801  IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
802 * The eigenvector support can become wrong
803 * because significant entries could be cut off due to a
804 * large GAPTOL parameter in LAR1V. Prevent this.
805  gaptol = zero
806  ELSE
807  gaptol = gap * eps
808  ENDIF
809  isupmn = in
810  isupmx = 1
811 * Update WGAP so that it holds the minimum gap
812 * to the left or the right. This is crucial in the
813 * case where bisection is used to ensure that the
814 * eigenvalue is refined up to the required precision.
815 * The correct value is restored afterwards.
816  savgap = wgap(windex)
817  wgap(windex) = gap
818 * We want to use the Rayleigh Quotient Correction
819 * as often as possible since it converges quadratically
820 * when we are close enough to the desired eigenvalue.
821 * However, the Rayleigh Quotient can have the wrong sign
822 * and lead us away from the desired eigenvalue. In this
823 * case, the best we can do is to use bisection.
824  usedbs = .false.
825  usedrq = .false.
826 * Bisection is initially turned off unless it is forced
827  needbs = .NOT.tryrqc
828  120 continue
829 * Check if bisection should be used to refine eigenvalue
830  IF(needbs) THEN
831 * Take the bisection as new iterate
832  usedbs = .true.
833  itmp1 = iwork( iindr+windex )
834  offset = indexw( wbegin ) - 1
835  CALL dlarrb( in, d(ibegin),
836  $ work(indlld+ibegin-1),indeig,indeig,
837  $ zero, two*eps, offset,
838  $ work(wbegin),wgap(wbegin),
839  $ werr(wbegin),work( indwrk ),
840  $ iwork( iindwk ), pivmin, spdiam,
841  $ itmp1, iinfo )
842  IF( iinfo.NE.0 ) THEN
843  info = -3
844  return
845  ENDIF
846  lambda = work( windex )
847 * Reset twist index from inaccurate LAMBDA to
848 * force computation of true MINGMA
849  iwork( iindr+windex ) = 0
850  ENDIF
851 * Given LAMBDA, compute the eigenvector.
852  CALL dlar1v( in, 1, in, lambda, d( ibegin ),
853  $ l( ibegin ), work(indld+ibegin-1),
854  $ work(indlld+ibegin-1),
855  $ pivmin, gaptol, z( ibegin, windex ),
856  $ .NOT.usedbs, negcnt, ztz, mingma,
857  $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
858  $ nrminv, resid, rqcorr, work( indwrk ) )
859  IF(iter .EQ. 0) THEN
860  bstres = resid
861  bstw = lambda
862  elseif(resid.LT.bstres) THEN
863  bstres = resid
864  bstw = lambda
865  ENDIF
866  isupmn = min(isupmn,isuppz( 2*windex-1 ))
867  isupmx = max(isupmx,isuppz( 2*windex ))
868  iter = iter + 1
869 
870 * sin alpha <= |resid|/gap
871 * Note that both the residual and the gap are
872 * proportional to the matrix, so ||T|| doesn't play
873 * a role in the quotient
874 
875 *
876 * Convergence test for Rayleigh-Quotient iteration
877 * (omitted when Bisection has been used)
878 *
879  IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
880  $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
881  $ THEN
882 * We need to check that the RQCORR update doesn't
883 * move the eigenvalue away from the desired one and
884 * towards a neighbor. -> protection with bisection
885  IF(indeig.LE.negcnt) THEN
886 * The wanted eigenvalue lies to the left
887  sgndef = -one
888  ELSE
889 * The wanted eigenvalue lies to the right
890  sgndef = one
891  ENDIF
892 * We only use the RQCORR if it improves the
893 * the iterate reasonably.
894  IF( ( rqcorr*sgndef.GE.zero )
895  $ .AND.( lambda + rqcorr.LE. right)
896  $ .AND.( lambda + rqcorr.GE. left)
897  $ ) THEN
898  usedrq = .true.
899 * Store new midpoint of bisection interval in WORK
900  IF(sgndef.EQ.one) THEN
901 * The current LAMBDA is on the left of the true
902 * eigenvalue
903  left = lambda
904 * We prefer to assume that the error estimate
905 * is correct. We could make the interval not
906 * as a bracket but to be modified if the RQCORR
907 * chooses to. In this case, the RIGHT side should
908 * be modified as follows:
909 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
910  ELSE
911 * The current LAMBDA is on the right of the true
912 * eigenvalue
913  right = lambda
914 * See comment about assuming the error estimate is
915 * correct above.
916 * LEFT = MIN(LEFT, LAMBDA + RQCORR)
917  ENDIF
918  work( windex ) =
919  $ half * (right + left)
920 * Take RQCORR since it has the correct sign and
921 * improves the iterate reasonably
922  lambda = lambda + rqcorr
923 * Update width of error interval
924  werr( windex ) =
925  $ half * (right-left)
926  ELSE
927  needbs = .true.
928  ENDIF
929  IF(right-left.LT.rqtol*abs(lambda)) THEN
930 * The eigenvalue is computed to bisection accuracy
931 * compute eigenvector and stop
932  usedbs = .true.
933  goto 120
934  elseif( iter.LT.maxitr ) THEN
935  goto 120
936  elseif( iter.EQ.maxitr ) THEN
937  needbs = .true.
938  goto 120
939  ELSE
940  info = 5
941  return
942  END IF
943  ELSE
944  stp2ii = .false.
945  IF(usedrq .AND. usedbs .AND.
946  $ bstres.LE.resid) THEN
947  lambda = bstw
948  stp2ii = .true.
949  ENDIF
950  IF (stp2ii) THEN
951 * improve error angle by second step
952  CALL dlar1v( in, 1, in, lambda,
953  $ d( ibegin ), l( ibegin ),
954  $ work(indld+ibegin-1),
955  $ work(indlld+ibegin-1),
956  $ pivmin, gaptol, z( ibegin, windex ),
957  $ .NOT.usedbs, negcnt, ztz, mingma,
958  $ iwork( iindr+windex ),
959  $ isuppz( 2*windex-1 ),
960  $ nrminv, resid, rqcorr, work( indwrk ) )
961  ENDIF
962  work( windex ) = lambda
963  END IF
964 *
965 * Compute FP-vector support w.r.t. whole matrix
966 *
967  isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
968  isuppz( 2*windex ) = isuppz( 2*windex )+oldien
969  zfrom = isuppz( 2*windex-1 )
970  zto = isuppz( 2*windex )
971  isupmn = isupmn + oldien
972  isupmx = isupmx + oldien
973 * Ensure vector is ok if support in the RQI has changed
974  IF(isupmn.LT.zfrom) THEN
975  DO 122 ii = isupmn,zfrom-1
976  z( ii, windex ) = zero
977  122 continue
978  ENDIF
979  IF(isupmx.GT.zto) THEN
980  DO 123 ii = zto+1,isupmx
981  z( ii, windex ) = zero
982  123 continue
983  ENDIF
984  CALL dscal( zto-zfrom+1, nrminv,
985  $ z( zfrom, windex ), 1 )
986  125 continue
987 * Update W
988  w( windex ) = lambda+sigma
989 * Recompute the gaps on the left and right
990 * But only allow them to become larger and not
991 * smaller (which can only happen through "bad"
992 * cancellation and doesn't reflect the theory
993 * where the initial gaps are underestimated due
994 * to WERR being too crude.)
995  IF(.NOT.eskip) THEN
996  IF( k.GT.1) THEN
997  wgap( windmn ) = max( wgap(windmn),
998  $ w(windex)-werr(windex)
999  $ - w(windmn)-werr(windmn) )
1000  ENDIF
1001  IF( windex.LT.wend ) THEN
1002  wgap( windex ) = max( savgap,
1003  $ w( windpl )-werr( windpl )
1004  $ - w( windex )-werr( windex) )
1005  ENDIF
1006  ENDIF
1007  idone = idone + 1
1008  ENDIF
1009 * here ends the code for the current child
1010 *
1011  139 continue
1012 * Proceed to any remaining child nodes
1013  newfst = j + 1
1014  140 continue
1015  150 continue
1016  ndepth = ndepth + 1
1017  go to 40
1018  END IF
1019  ibegin = iend + 1
1020  wbegin = wend + 1
1021  170 continue
1022 *
1023 
1024  return
1025 *
1026 * End of DLARRV
1027 *
1028  END