LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zlarrv.f
Go to the documentation of this file.
1 *> \brief \b ZLARRV 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 ZLARRV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarrv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarrv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarrv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLARRV( 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 * COMPLEX*16 Z( LDZ, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> ZLARRV 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 COMPLEX*16 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 ZLARRV.
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 complex16OTHERauxiliary
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 zlarrv( 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  COMPLEX*16 z( ldz, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  INTEGER maxitr
307  parameter( maxitr = 10 )
308  COMPLEX*16 czero
309  parameter( czero = ( 0.0d0, 0.0d0 ) )
310  DOUBLE PRECISION zero, one, two, three, four, half
311  parameter( zero = 0.0d0, one = 1.0d0,
312  $ two = 2.0d0, three = 3.0d0,
313  $ four = 4.0d0, half = 0.5d0)
314 * ..
315 * .. Local Scalars ..
316  LOGICAL eskip, needbs, stp2ii, tryrqc, usedbs, usedrq
317  INTEGER done, i, ibegin, idone, iend, ii, iindc1,
318  $ iindc2, iindr, iindwk, iinfo, im, in, indeig,
319  $ indld, indlld, indwrk, isupmn, isupmx, iter,
320  $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
321  $ ndepth, negcnt, newcls, newfst, newftt, newlst,
322  $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
323  $ oldncl, p, parity, q, wbegin, wend, windex,
324  $ windmn, windpl, zfrom, zto, zusedl, zusedu,
325  $ zusedw
326  INTEGER indin1, indin2
327  DOUBLE PRECISION bstres, bstw, eps, fudge, gap, gaptol, gl, gu,
328  $ lambda, left, lgap, mingma, nrminv, resid,
329  $ rgap, right, rqcorr, rqtol, savgap, sgndef,
330  $ sigma, spdiam, ssigma, tau, tmp, tol, ztz
331 * ..
332 * .. External Functions ..
333  DOUBLE PRECISION dlamch
334  EXTERNAL dlamch
335 * ..
336 * .. External Subroutines ..
337  EXTERNAL dcopy, dlarrb, dlarrf, zdscal, zlar1v,
338  $ zlaset
339 * ..
340 * .. Intrinsic Functions ..
341  INTRINSIC abs, dble, max, min
342  INTRINSIC dcmplx
343 * ..
344 * .. Executable Statements ..
345 * ..
346 
347 * The first N entries of WORK are reserved for the eigenvalues
348  indld = n+1
349  indlld= 2*n+1
350  indin1 = 3*n + 1
351  indin2 = 4*n + 1
352  indwrk = 5*n + 1
353  minwsize = 12 * n
354 
355  DO 5 i= 1,minwsize
356  work( i ) = zero
357  5 continue
358 
359 * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
360 * factorization used to compute the FP vector
361  iindr = 0
362 * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
363 * layer and the one above.
364  iindc1 = n
365  iindc2 = 2*n
366  iindwk = 3*n + 1
367 
368  miniwsize = 7 * n
369  DO 10 i= 1,miniwsize
370  iwork( i ) = 0
371  10 continue
372 
373  zusedl = 1
374  IF(dol.GT.1) THEN
375 * Set lower bound for use of Z
376  zusedl = dol-1
377  ENDIF
378  zusedu = m
379  IF(dou.LT.m) THEN
380 * Set lower bound for use of Z
381  zusedu = dou+1
382  ENDIF
383 * The width of the part of Z that is used
384  zusedw = zusedu - zusedl + 1
385 
386 
387  CALL zlaset( 'Full', n, zusedw, czero, czero,
388  $ z(1,zusedl), ldz )
389 
390  eps = dlamch( 'Precision' )
391  rqtol = two * eps
392 *
393 * Set expert flags for standard code.
394  tryrqc = .true.
395 
396  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
397  ELSE
398 * Only selected eigenpairs are computed. Since the other evalues
399 * are not refined by RQ iteration, bisection has to compute to full
400 * accuracy.
401  rtol1 = four * eps
402  rtol2 = four * eps
403  ENDIF
404 
405 * The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
406 * desired eigenvalues. The support of the nonzero eigenvector
407 * entries is contained in the interval IBEGIN:IEND.
408 * Remark that if k eigenpairs are desired, then the eigenvectors
409 * are stored in k contiguous columns of Z.
410 
411 * DONE is the number of eigenvectors already computed
412  done = 0
413  ibegin = 1
414  wbegin = 1
415  DO 170 jblk = 1, iblock( m )
416  iend = isplit( jblk )
417  sigma = l( iend )
418 * Find the eigenvectors of the submatrix indexed IBEGIN
419 * through IEND.
420  wend = wbegin - 1
421  15 continue
422  IF( wend.LT.m ) THEN
423  IF( iblock( wend+1 ).EQ.jblk ) THEN
424  wend = wend + 1
425  go to 15
426  END IF
427  END IF
428  IF( wend.LT.wbegin ) THEN
429  ibegin = iend + 1
430  go to 170
431  elseif( (wend.LT.dol).OR.(wbegin.GT.dou) ) THEN
432  ibegin = iend + 1
433  wbegin = wend + 1
434  go to 170
435  END IF
436 
437 * Find local spectral diameter of the block
438  gl = gers( 2*ibegin-1 )
439  gu = gers( 2*ibegin )
440  DO 20 i = ibegin+1 , iend
441  gl = min( gers( 2*i-1 ), gl )
442  gu = max( gers( 2*i ), gu )
443  20 continue
444  spdiam = gu - gl
445 
446 * OLDIEN is the last index of the previous block
447  oldien = ibegin - 1
448 * Calculate the size of the current block
449  in = iend - ibegin + 1
450 * The number of eigenvalues in the current block
451  im = wend - wbegin + 1
452 
453 * This is for a 1x1 block
454  IF( ibegin.EQ.iend ) THEN
455  done = done+1
456  z( ibegin, wbegin ) = dcmplx( one, zero )
457  isuppz( 2*wbegin-1 ) = ibegin
458  isuppz( 2*wbegin ) = ibegin
459  w( wbegin ) = w( wbegin ) + sigma
460  work( wbegin ) = w( wbegin )
461  ibegin = iend + 1
462  wbegin = wbegin + 1
463  go to 170
464  END IF
465 
466 * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
467 * Note that these can be approximations, in this case, the corresp.
468 * entries of WERR give the size of the uncertainty interval.
469 * The eigenvalue approximations will be refined when necessary as
470 * high relative accuracy is required for the computation of the
471 * corresponding eigenvectors.
472  CALL dcopy( im, w( wbegin ), 1,
473  $ work( wbegin ), 1 )
474 
475 * We store in W the eigenvalue approximations w.r.t. the original
476 * matrix T.
477  DO 30 i=1,im
478  w(wbegin+i-1) = w(wbegin+i-1)+sigma
479  30 continue
480 
481 
482 * NDEPTH is the current depth of the representation tree
483  ndepth = 0
484 * PARITY is either 1 or 0
485  parity = 1
486 * NCLUS is the number of clusters for the next level of the
487 * representation tree, we start with NCLUS = 1 for the root
488  nclus = 1
489  iwork( iindc1+1 ) = 1
490  iwork( iindc1+2 ) = im
491 
492 * IDONE is the number of eigenvectors already computed in the current
493 * block
494  idone = 0
495 * loop while( IDONE.LT.IM )
496 * generate the representation tree for the current block and
497 * compute the eigenvectors
498  40 continue
499  IF( idone.LT.im ) THEN
500 * This is a crude protection against infinitely deep trees
501  IF( ndepth.GT.m ) THEN
502  info = -2
503  return
504  ENDIF
505 * breadth first processing of the current level of the representation
506 * tree: OLDNCL = number of clusters on current level
507  oldncl = nclus
508 * reset NCLUS to count the number of child clusters
509  nclus = 0
510 *
511  parity = 1 - parity
512  IF( parity.EQ.0 ) THEN
513  oldcls = iindc1
514  newcls = iindc2
515  ELSE
516  oldcls = iindc2
517  newcls = iindc1
518  END IF
519 * Process the clusters on the current level
520  DO 150 i = 1, oldncl
521  j = oldcls + 2*i
522 * OLDFST, OLDLST = first, last index of current cluster.
523 * cluster indices start with 1 and are relative
524 * to WBEGIN when accessing W, WGAP, WERR, Z
525  oldfst = iwork( j-1 )
526  oldlst = iwork( j )
527  IF( ndepth.GT.0 ) THEN
528 * Retrieve relatively robust representation (RRR) of cluster
529 * that has been computed at the previous level
530 * The RRR is stored in Z and overwritten once the eigenvectors
531 * have been computed or when the cluster is refined
532 
533  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
534 * Get representation from location of the leftmost evalue
535 * of the cluster
536  j = wbegin + oldfst - 1
537  ELSE
538  IF(wbegin+oldfst-1.LT.dol) THEN
539 * Get representation from the left end of Z array
540  j = dol - 1
541  elseif(wbegin+oldfst-1.GT.dou) THEN
542 * Get representation from the right end of Z array
543  j = dou
544  ELSE
545  j = wbegin + oldfst - 1
546  ENDIF
547  ENDIF
548  DO 45 k = 1, in - 1
549  d( ibegin+k-1 ) = dble( z( ibegin+k-1,
550  $ j ) )
551  l( ibegin+k-1 ) = dble( z( ibegin+k-1,
552  $ j+1 ) )
553  45 continue
554  d( iend ) = dble( z( iend, j ) )
555  sigma = dble( z( iend, j+1 ) )
556 
557 * Set the corresponding entries in Z to zero
558  CALL zlaset( 'Full', in, 2, czero, czero,
559  $ z( ibegin, j), ldz )
560  END IF
561 
562 * Compute DL and DLL of current RRR
563  DO 50 j = ibegin, iend-1
564  tmp = d( j )*l( j )
565  work( indld-1+j ) = tmp
566  work( indlld-1+j ) = tmp*l( j )
567  50 continue
568 
569  IF( ndepth.GT.0 ) THEN
570 * P and Q are index of the first and last eigenvalue to compute
571 * within the current block
572  p = indexw( wbegin-1+oldfst )
573  q = indexw( wbegin-1+oldlst )
574 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
575 * through the Q-OFFSET elements of these arrays are to be used.
576 * OFFSET = P-OLDFST
577  offset = indexw( wbegin ) - 1
578 * perform limited bisection (if necessary) to get approximate
579 * eigenvalues to the precision needed.
580  CALL dlarrb( in, d( ibegin ),
581  $ work(indlld+ibegin-1),
582  $ p, q, rtol1, rtol2, offset,
583  $ work(wbegin),wgap(wbegin),werr(wbegin),
584  $ work( indwrk ), iwork( iindwk ),
585  $ pivmin, spdiam, in, iinfo )
586  IF( iinfo.NE.0 ) THEN
587  info = -1
588  return
589  ENDIF
590 * We also recompute the extremal gaps. W holds all eigenvalues
591 * of the unshifted matrix and must be used for computation
592 * of WGAP, the entries of WORK might stem from RRRs with
593 * different shifts. The gaps from WBEGIN-1+OLDFST to
594 * WBEGIN-1+OLDLST are correctly computed in DLARRB.
595 * However, we only allow the gaps to become greater since
596 * this is what should happen when we decrease WERR
597  IF( oldfst.GT.1) THEN
598  wgap( wbegin+oldfst-2 ) =
599  $ max(wgap(wbegin+oldfst-2),
600  $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
601  $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
602  ENDIF
603  IF( wbegin + oldlst -1 .LT. wend ) THEN
604  wgap( wbegin+oldlst-1 ) =
605  $ max(wgap(wbegin+oldlst-1),
606  $ w(wbegin+oldlst)-werr(wbegin+oldlst)
607  $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
608  ENDIF
609 * Each time the eigenvalues in WORK get refined, we store
610 * the newly found approximation with all shifts applied in W
611  DO 53 j=oldfst,oldlst
612  w(wbegin+j-1) = work(wbegin+j-1)+sigma
613  53 continue
614  END IF
615 
616 * Process the current node.
617  newfst = oldfst
618  DO 140 j = oldfst, oldlst
619  IF( j.EQ.oldlst ) THEN
620 * we are at the right end of the cluster, this is also the
621 * boundary of the child cluster
622  newlst = j
623  ELSE IF ( wgap( wbegin + j -1).GE.
624  $ minrgp* abs( work(wbegin + j -1) ) ) THEN
625 * the right relative gap is big enough, the child cluster
626 * (NEWFST,..,NEWLST) is well separated from the following
627  newlst = j
628  ELSE
629 * inside a child cluster, the relative gap is not
630 * big enough.
631  goto 140
632  END IF
633 
634 * Compute size of child cluster found
635  newsiz = newlst - newfst + 1
636 
637 * NEWFTT is the place in Z where the new RRR or the computed
638 * eigenvector is to be stored
639  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
640 * Store representation at location of the leftmost evalue
641 * of the cluster
642  newftt = wbegin + newfst - 1
643  ELSE
644  IF(wbegin+newfst-1.LT.dol) THEN
645 * Store representation at the left end of Z array
646  newftt = dol - 1
647  elseif(wbegin+newfst-1.GT.dou) THEN
648 * Store representation at the right end of Z array
649  newftt = dou
650  ELSE
651  newftt = wbegin + newfst - 1
652  ENDIF
653  ENDIF
654 
655  IF( newsiz.GT.1) THEN
656 *
657 * Current child is not a singleton but a cluster.
658 * Compute and store new representation of child.
659 *
660 *
661 * Compute left and right cluster gap.
662 *
663 * LGAP and RGAP are not computed from WORK because
664 * the eigenvalue approximations may stem from RRRs
665 * different shifts. However, W hold all eigenvalues
666 * of the unshifted matrix. Still, the entries in WGAP
667 * have to be computed from WORK since the entries
668 * in W might be of the same order so that gaps are not
669 * exhibited correctly for very close eigenvalues.
670  IF( newfst.EQ.1 ) THEN
671  lgap = max( zero,
672  $ w(wbegin)-werr(wbegin) - vl )
673  ELSE
674  lgap = wgap( wbegin+newfst-2 )
675  ENDIF
676  rgap = wgap( wbegin+newlst-1 )
677 *
678 * Compute left- and rightmost eigenvalue of child
679 * to high precision in order to shift as close
680 * as possible and obtain as large relative gaps
681 * as possible
682 *
683  DO 55 k =1,2
684  IF(k.EQ.1) THEN
685  p = indexw( wbegin-1+newfst )
686  ELSE
687  p = indexw( wbegin-1+newlst )
688  ENDIF
689  offset = indexw( wbegin ) - 1
690  CALL dlarrb( in, d(ibegin),
691  $ work( indlld+ibegin-1 ),p,p,
692  $ rqtol, rqtol, offset,
693  $ work(wbegin),wgap(wbegin),
694  $ werr(wbegin),work( indwrk ),
695  $ iwork( iindwk ), pivmin, spdiam,
696  $ in, iinfo )
697  55 continue
698 *
699  IF((wbegin+newlst-1.LT.dol).OR.
700  $ (wbegin+newfst-1.GT.dou)) THEN
701 * if the cluster contains no desired eigenvalues
702 * skip the computation of that branch of the rep. tree
703 *
704 * We could skip before the refinement of the extremal
705 * eigenvalues of the child, but then the representation
706 * tree could be different from the one when nothing is
707 * skipped. For this reason we skip at this place.
708  idone = idone + newlst - newfst + 1
709  goto 139
710  ENDIF
711 *
712 * Compute RRR of child cluster.
713 * Note that the new RRR is stored in Z
714 *
715 * DLARRF needs LWORK = 2*N
716  CALL dlarrf( in, d( ibegin ), l( ibegin ),
717  $ work(indld+ibegin-1),
718  $ newfst, newlst, work(wbegin),
719  $ wgap(wbegin), werr(wbegin),
720  $ spdiam, lgap, rgap, pivmin, tau,
721  $ work( indin1 ), work( indin2 ),
722  $ work( indwrk ), iinfo )
723 * In the complex case, DLARRF cannot write
724 * the new RRR directly into Z and needs an intermediate
725 * workspace
726  DO 56 k = 1, in-1
727  z( ibegin+k-1, newftt ) =
728  $ dcmplx( work( indin1+k-1 ), zero )
729  z( ibegin+k-1, newftt+1 ) =
730  $ dcmplx( work( indin2+k-1 ), zero )
731  56 continue
732  z( iend, newftt ) =
733  $ dcmplx( work( indin1+in-1 ), zero )
734  IF( iinfo.EQ.0 ) THEN
735 * a new RRR for the cluster was found by DLARRF
736 * update shift and store it
737  ssigma = sigma + tau
738  z( iend, newftt+1 ) = dcmplx( ssigma, zero )
739 * WORK() are the midpoints and WERR() the semi-width
740 * Note that the entries in W are unchanged.
741  DO 116 k = newfst, newlst
742  fudge =
743  $ three*eps*abs(work(wbegin+k-1))
744  work( wbegin + k - 1 ) =
745  $ work( wbegin + k - 1) - tau
746  fudge = fudge +
747  $ four*eps*abs(work(wbegin+k-1))
748 * Fudge errors
749  werr( wbegin + k - 1 ) =
750  $ werr( wbegin + k - 1 ) + fudge
751 * Gaps are not fudged. Provided that WERR is small
752 * when eigenvalues are close, a zero gap indicates
753 * that a new representation is needed for resolving
754 * the cluster. A fudge could lead to a wrong decision
755 * of judging eigenvalues 'separated' which in
756 * reality are not. This could have a negative impact
757 * on the orthogonality of the computed eigenvectors.
758  116 continue
759 
760  nclus = nclus + 1
761  k = newcls + 2*nclus
762  iwork( k-1 ) = newfst
763  iwork( k ) = newlst
764  ELSE
765  info = -2
766  return
767  ENDIF
768  ELSE
769 *
770 * Compute eigenvector of singleton
771 *
772  iter = 0
773 *
774  tol = four * log(dble(in)) * eps
775 *
776  k = newfst
777  windex = wbegin + k - 1
778  windmn = max(windex - 1,1)
779  windpl = min(windex + 1,m)
780  lambda = work( windex )
781  done = done + 1
782 * Check if eigenvector computation is to be skipped
783  IF((windex.LT.dol).OR.
784  $ (windex.GT.dou)) THEN
785  eskip = .true.
786  goto 125
787  ELSE
788  eskip = .false.
789  ENDIF
790  left = work( windex ) - werr( windex )
791  right = work( windex ) + werr( windex )
792  indeig = indexw( windex )
793 * Note that since we compute the eigenpairs for a child,
794 * all eigenvalue approximations are w.r.t the same shift.
795 * In this case, the entries in WORK should be used for
796 * computing the gaps since they exhibit even very small
797 * differences in the eigenvalues, as opposed to the
798 * entries in W which might "look" the same.
799 
800  IF( k .EQ. 1) THEN
801 * In the case RANGE='I' and with not much initial
802 * accuracy in LAMBDA and VL, the formula
803 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
804 * can lead to an overestimation of the left gap and
805 * thus to inadequately early RQI 'convergence'.
806 * Prevent this by forcing a small left gap.
807  lgap = eps*max(abs(left),abs(right))
808  ELSE
809  lgap = wgap(windmn)
810  ENDIF
811  IF( k .EQ. im) THEN
812 * In the case RANGE='I' and with not much initial
813 * accuracy in LAMBDA and VU, the formula
814 * can lead to an overestimation of the right gap and
815 * thus to inadequately early RQI 'convergence'.
816 * Prevent this by forcing a small right gap.
817  rgap = eps*max(abs(left),abs(right))
818  ELSE
819  rgap = wgap(windex)
820  ENDIF
821  gap = min( lgap, rgap )
822  IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
823 * The eigenvector support can become wrong
824 * because significant entries could be cut off due to a
825 * large GAPTOL parameter in LAR1V. Prevent this.
826  gaptol = zero
827  ELSE
828  gaptol = gap * eps
829  ENDIF
830  isupmn = in
831  isupmx = 1
832 * Update WGAP so that it holds the minimum gap
833 * to the left or the right. This is crucial in the
834 * case where bisection is used to ensure that the
835 * eigenvalue is refined up to the required precision.
836 * The correct value is restored afterwards.
837  savgap = wgap(windex)
838  wgap(windex) = gap
839 * We want to use the Rayleigh Quotient Correction
840 * as often as possible since it converges quadratically
841 * when we are close enough to the desired eigenvalue.
842 * However, the Rayleigh Quotient can have the wrong sign
843 * and lead us away from the desired eigenvalue. In this
844 * case, the best we can do is to use bisection.
845  usedbs = .false.
846  usedrq = .false.
847 * Bisection is initially turned off unless it is forced
848  needbs = .NOT.tryrqc
849  120 continue
850 * Check if bisection should be used to refine eigenvalue
851  IF(needbs) THEN
852 * Take the bisection as new iterate
853  usedbs = .true.
854  itmp1 = iwork( iindr+windex )
855  offset = indexw( wbegin ) - 1
856  CALL dlarrb( in, d(ibegin),
857  $ work(indlld+ibegin-1),indeig,indeig,
858  $ zero, two*eps, offset,
859  $ work(wbegin),wgap(wbegin),
860  $ werr(wbegin),work( indwrk ),
861  $ iwork( iindwk ), pivmin, spdiam,
862  $ itmp1, iinfo )
863  IF( iinfo.NE.0 ) THEN
864  info = -3
865  return
866  ENDIF
867  lambda = work( windex )
868 * Reset twist index from inaccurate LAMBDA to
869 * force computation of true MINGMA
870  iwork( iindr+windex ) = 0
871  ENDIF
872 * Given LAMBDA, compute the eigenvector.
873  CALL zlar1v( in, 1, in, lambda, d( ibegin ),
874  $ l( ibegin ), work(indld+ibegin-1),
875  $ work(indlld+ibegin-1),
876  $ pivmin, gaptol, z( ibegin, windex ),
877  $ .NOT.usedbs, negcnt, ztz, mingma,
878  $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
879  $ nrminv, resid, rqcorr, work( indwrk ) )
880  IF(iter .EQ. 0) THEN
881  bstres = resid
882  bstw = lambda
883  elseif(resid.LT.bstres) THEN
884  bstres = resid
885  bstw = lambda
886  ENDIF
887  isupmn = min(isupmn,isuppz( 2*windex-1 ))
888  isupmx = max(isupmx,isuppz( 2*windex ))
889  iter = iter + 1
890 
891 * sin alpha <= |resid|/gap
892 * Note that both the residual and the gap are
893 * proportional to the matrix, so ||T|| doesn't play
894 * a role in the quotient
895 
896 *
897 * Convergence test for Rayleigh-Quotient iteration
898 * (omitted when Bisection has been used)
899 *
900  IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
901  $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
902  $ THEN
903 * We need to check that the RQCORR update doesn't
904 * move the eigenvalue away from the desired one and
905 * towards a neighbor. -> protection with bisection
906  IF(indeig.LE.negcnt) THEN
907 * The wanted eigenvalue lies to the left
908  sgndef = -one
909  ELSE
910 * The wanted eigenvalue lies to the right
911  sgndef = one
912  ENDIF
913 * We only use the RQCORR if it improves the
914 * the iterate reasonably.
915  IF( ( rqcorr*sgndef.GE.zero )
916  $ .AND.( lambda + rqcorr.LE. right)
917  $ .AND.( lambda + rqcorr.GE. left)
918  $ ) THEN
919  usedrq = .true.
920 * Store new midpoint of bisection interval in WORK
921  IF(sgndef.EQ.one) THEN
922 * The current LAMBDA is on the left of the true
923 * eigenvalue
924  left = lambda
925 * We prefer to assume that the error estimate
926 * is correct. We could make the interval not
927 * as a bracket but to be modified if the RQCORR
928 * chooses to. In this case, the RIGHT side should
929 * be modified as follows:
930 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
931  ELSE
932 * The current LAMBDA is on the right of the true
933 * eigenvalue
934  right = lambda
935 * See comment about assuming the error estimate is
936 * correct above.
937 * LEFT = MIN(LEFT, LAMBDA + RQCORR)
938  ENDIF
939  work( windex ) =
940  $ half * (right + left)
941 * Take RQCORR since it has the correct sign and
942 * improves the iterate reasonably
943  lambda = lambda + rqcorr
944 * Update width of error interval
945  werr( windex ) =
946  $ half * (right-left)
947  ELSE
948  needbs = .true.
949  ENDIF
950  IF(right-left.LT.rqtol*abs(lambda)) THEN
951 * The eigenvalue is computed to bisection accuracy
952 * compute eigenvector and stop
953  usedbs = .true.
954  goto 120
955  elseif( iter.LT.maxitr ) THEN
956  goto 120
957  elseif( iter.EQ.maxitr ) THEN
958  needbs = .true.
959  goto 120
960  ELSE
961  info = 5
962  return
963  END IF
964  ELSE
965  stp2ii = .false.
966  IF(usedrq .AND. usedbs .AND.
967  $ bstres.LE.resid) THEN
968  lambda = bstw
969  stp2ii = .true.
970  ENDIF
971  IF (stp2ii) THEN
972 * improve error angle by second step
973  CALL zlar1v( in, 1, in, lambda,
974  $ d( ibegin ), l( ibegin ),
975  $ work(indld+ibegin-1),
976  $ work(indlld+ibegin-1),
977  $ pivmin, gaptol, z( ibegin, windex ),
978  $ .NOT.usedbs, negcnt, ztz, mingma,
979  $ iwork( iindr+windex ),
980  $ isuppz( 2*windex-1 ),
981  $ nrminv, resid, rqcorr, work( indwrk ) )
982  ENDIF
983  work( windex ) = lambda
984  END IF
985 *
986 * Compute FP-vector support w.r.t. whole matrix
987 *
988  isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
989  isuppz( 2*windex ) = isuppz( 2*windex )+oldien
990  zfrom = isuppz( 2*windex-1 )
991  zto = isuppz( 2*windex )
992  isupmn = isupmn + oldien
993  isupmx = isupmx + oldien
994 * Ensure vector is ok if support in the RQI has changed
995  IF(isupmn.LT.zfrom) THEN
996  DO 122 ii = isupmn,zfrom-1
997  z( ii, windex ) = zero
998  122 continue
999  ENDIF
1000  IF(isupmx.GT.zto) THEN
1001  DO 123 ii = zto+1,isupmx
1002  z( ii, windex ) = zero
1003  123 continue
1004  ENDIF
1005  CALL zdscal( zto-zfrom+1, nrminv,
1006  $ z( zfrom, windex ), 1 )
1007  125 continue
1008 * Update W
1009  w( windex ) = lambda+sigma
1010 * Recompute the gaps on the left and right
1011 * But only allow them to become larger and not
1012 * smaller (which can only happen through "bad"
1013 * cancellation and doesn't reflect the theory
1014 * where the initial gaps are underestimated due
1015 * to WERR being too crude.)
1016  IF(.NOT.eskip) THEN
1017  IF( k.GT.1) THEN
1018  wgap( windmn ) = max( wgap(windmn),
1019  $ w(windex)-werr(windex)
1020  $ - w(windmn)-werr(windmn) )
1021  ENDIF
1022  IF( windex.LT.wend ) THEN
1023  wgap( windex ) = max( savgap,
1024  $ w( windpl )-werr( windpl )
1025  $ - w( windex )-werr( windex) )
1026  ENDIF
1027  ENDIF
1028  idone = idone + 1
1029  ENDIF
1030 * here ends the code for the current child
1031 *
1032  139 continue
1033 * Proceed to any remaining child nodes
1034  newfst = j + 1
1035  140 continue
1036  150 continue
1037  ndepth = ndepth + 1
1038  go to 40
1039  END IF
1040  ibegin = iend + 1
1041  wbegin = wend + 1
1042  170 continue
1043 *
1044 
1045  return
1046 *
1047 * End of ZLARRV
1048 *
1049  END