LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
clarrv.f
Go to the documentation of this file.
1 *> \brief \b CLARRV 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 CLARRV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarrv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarrv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarrv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLARRV( 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 * REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
33 * $ ISUPPZ( * ), IWORK( * )
34 * REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
35 * $ WGAP( * ), WORK( * )
36 * COMPLEX Z( LDZ, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> CLARRV 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 SLARRE.
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 REAL
62 *> \endverbatim
63 *>
64 *> \param[in] VU
65 *> \verbatim
66 *> VU is REAL
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 REAL 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 REAL 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 SLARRE.
86 *> On exit, L is overwritten.
87 *> \endverbatim
88 *>
89 *> \param[in] PIVMIN
90 *> \verbatim
91 *> PIVMIN is REAL
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 REAL
131 *> \endverbatim
132 *>
133 *> \param[in] RTOL1
134 *> \verbatim
135 *> RTOL1 is REAL
136 *> \endverbatim
137 *>
138 *> \param[in] RTOL2
139 *> \verbatim
140 *> RTOL2 is REAL
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 REAL 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 SLARRE 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 REAL 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 REAL 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 REAL 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 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 REAL 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 CLARRV.
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 SLARRB when refining a child's eigenvalues.
245 *> =-2: Problem in SLARRF 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 SLARRB 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 complexOTHERauxiliary
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 clarrv( 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  REAL minrgp, pivmin, rtol1, rtol2, vl, vu
294 * ..
295 * .. Array Arguments ..
296  INTEGER iblock( * ), indexw( * ), isplit( * ),
297  $ isuppz( * ), iwork( * )
298  REAL d( * ), gers( * ), l( * ), w( * ), werr( * ),
299  $ wgap( * ), work( * )
300  COMPLEX z( ldz, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  INTEGER maxitr
307  parameter( maxitr = 10 )
308  COMPLEX czero
309  parameter( czero = ( 0.0e0, 0.0e0 ) )
310  REAL zero, one, two, three, four, half
311  parameter( zero = 0.0e0, one = 1.0e0,
312  $ two = 2.0e0, three = 3.0e0,
313  $ four = 4.0e0, half = 0.5e0)
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  REAL 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  REAL slamch
334  EXTERNAL slamch
335 * ..
336 * .. External Subroutines ..
337  EXTERNAL clar1v, claset, csscal, scopy, slarrb,
338  $ slarrf
339 * ..
340 * .. Intrinsic Functions ..
341  INTRINSIC abs, REAL, max, min
342  INTRINSIC cmplx
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 claset( 'Full', n, zusedw, czero, czero,
388  $ z(1,zusedl), ldz )
389 
390  eps = slamch( '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 ) = cmplx( 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 scopy( 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 ) = REAL( Z( IBEGIN+K-1, $ J ) )
550  l( ibegin+k-1 ) = REAL( Z( IBEGIN+K-1, $ J+1 ) )
551  45 continue
552  d( iend ) = REAL( Z( IEND, J ) )
553  sigma = REAL( Z( IEND, J+1 ) )
554 
555 * Set the corresponding entries in Z to zero
556  CALL claset( 'Full', in, 2, czero, czero,
557  $ z( ibegin, j), ldz )
558  END IF
559 
560 * Compute DL and DLL of current RRR
561  DO 50 j = ibegin, iend-1
562  tmp = d( j )*l( j )
563  work( indld-1+j ) = tmp
564  work( indlld-1+j ) = tmp*l( j )
565  50 continue
566 
567  IF( ndepth.GT.0 ) THEN
568 * P and Q are index of the first and last eigenvalue to compute
569 * within the current block
570  p = indexw( wbegin-1+oldfst )
571  q = indexw( wbegin-1+oldlst )
572 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
573 * through the Q-OFFSET elements of these arrays are to be used.
574 * OFFSET = P-OLDFST
575  offset = indexw( wbegin ) - 1
576 * perform limited bisection (if necessary) to get approximate
577 * eigenvalues to the precision needed.
578  CALL slarrb( in, d( ibegin ),
579  $ work(indlld+ibegin-1),
580  $ p, q, rtol1, rtol2, offset,
581  $ work(wbegin),wgap(wbegin),werr(wbegin),
582  $ work( indwrk ), iwork( iindwk ),
583  $ pivmin, spdiam, in, iinfo )
584  IF( iinfo.NE.0 ) THEN
585  info = -1
586  return
587  ENDIF
588 * We also recompute the extremal gaps. W holds all eigenvalues
589 * of the unshifted matrix and must be used for computation
590 * of WGAP, the entries of WORK might stem from RRRs with
591 * different shifts. The gaps from WBEGIN-1+OLDFST to
592 * WBEGIN-1+OLDLST are correctly computed in SLARRB.
593 * However, we only allow the gaps to become greater since
594 * this is what should happen when we decrease WERR
595  IF( oldfst.GT.1) THEN
596  wgap( wbegin+oldfst-2 ) =
597  $ max(wgap(wbegin+oldfst-2),
598  $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
599  $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
600  ENDIF
601  IF( wbegin + oldlst -1 .LT. wend ) THEN
602  wgap( wbegin+oldlst-1 ) =
603  $ max(wgap(wbegin+oldlst-1),
604  $ w(wbegin+oldlst)-werr(wbegin+oldlst)
605  $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
606  ENDIF
607 * Each time the eigenvalues in WORK get refined, we store
608 * the newly found approximation with all shifts applied in W
609  DO 53 j=oldfst,oldlst
610  w(wbegin+j-1) = work(wbegin+j-1)+sigma
611  53 continue
612  END IF
613 
614 * Process the current node.
615  newfst = oldfst
616  DO 140 j = oldfst, oldlst
617  IF( j.EQ.oldlst ) THEN
618 * we are at the right end of the cluster, this is also the
619 * boundary of the child cluster
620  newlst = j
621  ELSE IF ( wgap( wbegin + j -1).GE.
622  $ minrgp* abs( work(wbegin + j -1) ) ) THEN
623 * the right relative gap is big enough, the child cluster
624 * (NEWFST,..,NEWLST) is well separated from the following
625  newlst = j
626  ELSE
627 * inside a child cluster, the relative gap is not
628 * big enough.
629  goto 140
630  END IF
631 
632 * Compute size of child cluster found
633  newsiz = newlst - newfst + 1
634 
635 * NEWFTT is the place in Z where the new RRR or the computed
636 * eigenvector is to be stored
637  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
638 * Store representation at location of the leftmost evalue
639 * of the cluster
640  newftt = wbegin + newfst - 1
641  ELSE
642  IF(wbegin+newfst-1.LT.dol) THEN
643 * Store representation at the left end of Z array
644  newftt = dol - 1
645  elseif(wbegin+newfst-1.GT.dou) THEN
646 * Store representation at the right end of Z array
647  newftt = dou
648  ELSE
649  newftt = wbegin + newfst - 1
650  ENDIF
651  ENDIF
652 
653  IF( newsiz.GT.1) THEN
654 *
655 * Current child is not a singleton but a cluster.
656 * Compute and store new representation of child.
657 *
658 *
659 * Compute left and right cluster gap.
660 *
661 * LGAP and RGAP are not computed from WORK because
662 * the eigenvalue approximations may stem from RRRs
663 * different shifts. However, W hold all eigenvalues
664 * of the unshifted matrix. Still, the entries in WGAP
665 * have to be computed from WORK since the entries
666 * in W might be of the same order so that gaps are not
667 * exhibited correctly for very close eigenvalues.
668  IF( newfst.EQ.1 ) THEN
669  lgap = max( zero,
670  $ w(wbegin)-werr(wbegin) - vl )
671  ELSE
672  lgap = wgap( wbegin+newfst-2 )
673  ENDIF
674  rgap = wgap( wbegin+newlst-1 )
675 *
676 * Compute left- and rightmost eigenvalue of child
677 * to high precision in order to shift as close
678 * as possible and obtain as large relative gaps
679 * as possible
680 *
681  DO 55 k =1,2
682  IF(k.EQ.1) THEN
683  p = indexw( wbegin-1+newfst )
684  ELSE
685  p = indexw( wbegin-1+newlst )
686  ENDIF
687  offset = indexw( wbegin ) - 1
688  CALL slarrb( in, d(ibegin),
689  $ work( indlld+ibegin-1 ),p,p,
690  $ rqtol, rqtol, offset,
691  $ work(wbegin),wgap(wbegin),
692  $ werr(wbegin),work( indwrk ),
693  $ iwork( iindwk ), pivmin, spdiam,
694  $ in, iinfo )
695  55 continue
696 *
697  IF((wbegin+newlst-1.LT.dol).OR.
698  $ (wbegin+newfst-1.GT.dou)) THEN
699 * if the cluster contains no desired eigenvalues
700 * skip the computation of that branch of the rep. tree
701 *
702 * We could skip before the refinement of the extremal
703 * eigenvalues of the child, but then the representation
704 * tree could be different from the one when nothing is
705 * skipped. For this reason we skip at this place.
706  idone = idone + newlst - newfst + 1
707  goto 139
708  ENDIF
709 *
710 * Compute RRR of child cluster.
711 * Note that the new RRR is stored in Z
712 *
713 * SLARRF needs LWORK = 2*N
714  CALL slarrf( in, d( ibegin ), l( ibegin ),
715  $ work(indld+ibegin-1),
716  $ newfst, newlst, work(wbegin),
717  $ wgap(wbegin), werr(wbegin),
718  $ spdiam, lgap, rgap, pivmin, tau,
719  $ work( indin1 ), work( indin2 ),
720  $ work( indwrk ), iinfo )
721 * In the complex case, SLARRF cannot write
722 * the new RRR directly into Z and needs an intermediate
723 * workspace
724  DO 56 k = 1, in-1
725  z( ibegin+k-1, newftt ) =
726  $ cmplx( work( indin1+k-1 ), zero )
727  z( ibegin+k-1, newftt+1 ) =
728  $ cmplx( work( indin2+k-1 ), zero )
729  56 continue
730  z( iend, newftt ) =
731  $ cmplx( work( indin1+in-1 ), zero )
732  IF( iinfo.EQ.0 ) THEN
733 * a new RRR for the cluster was found by SLARRF
734 * update shift and store it
735  ssigma = sigma + tau
736  z( iend, newftt+1 ) = cmplx( ssigma, zero )
737 * WORK() are the midpoints and WERR() the semi-width
738 * Note that the entries in W are unchanged.
739  DO 116 k = newfst, newlst
740  fudge =
741  $ three*eps*abs(work(wbegin+k-1))
742  work( wbegin + k - 1 ) =
743  $ work( wbegin + k - 1) - tau
744  fudge = fudge +
745  $ four*eps*abs(work(wbegin+k-1))
746 * Fudge errors
747  werr( wbegin + k - 1 ) =
748  $ werr( wbegin + k - 1 ) + fudge
749 * Gaps are not fudged. Provided that WERR is small
750 * when eigenvalues are close, a zero gap indicates
751 * that a new representation is needed for resolving
752 * the cluster. A fudge could lead to a wrong decision
753 * of judging eigenvalues 'separated' which in
754 * reality are not. This could have a negative impact
755 * on the orthogonality of the computed eigenvectors.
756  116 continue
757 
758  nclus = nclus + 1
759  k = newcls + 2*nclus
760  iwork( k-1 ) = newfst
761  iwork( k ) = newlst
762  ELSE
763  info = -2
764  return
765  ENDIF
766  ELSE
767 *
768 * Compute eigenvector of singleton
769 *
770  iter = 0
771 *
772  tol = four * log(REAL(in)) * eps
773 *
774  k = newfst
775  windex = wbegin + k - 1
776  windmn = max(windex - 1,1)
777  windpl = min(windex + 1,m)
778  lambda = work( windex )
779  done = done + 1
780 * Check if eigenvector computation is to be skipped
781  IF((windex.LT.dol).OR.
782  $ (windex.GT.dou)) THEN
783  eskip = .true.
784  goto 125
785  ELSE
786  eskip = .false.
787  ENDIF
788  left = work( windex ) - werr( windex )
789  right = work( windex ) + werr( windex )
790  indeig = indexw( windex )
791 * Note that since we compute the eigenpairs for a child,
792 * all eigenvalue approximations are w.r.t the same shift.
793 * In this case, the entries in WORK should be used for
794 * computing the gaps since they exhibit even very small
795 * differences in the eigenvalues, as opposed to the
796 * entries in W which might "look" the same.
797 
798  IF( k .EQ. 1) THEN
799 * In the case RANGE='I' and with not much initial
800 * accuracy in LAMBDA and VL, the formula
801 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
802 * can lead to an overestimation of the left gap and
803 * thus to inadequately early RQI 'convergence'.
804 * Prevent this by forcing a small left gap.
805  lgap = eps*max(abs(left),abs(right))
806  ELSE
807  lgap = wgap(windmn)
808  ENDIF
809  IF( k .EQ. im) THEN
810 * In the case RANGE='I' and with not much initial
811 * accuracy in LAMBDA and VU, the formula
812 * can lead to an overestimation of the right gap and
813 * thus to inadequately early RQI 'convergence'.
814 * Prevent this by forcing a small right gap.
815  rgap = eps*max(abs(left),abs(right))
816  ELSE
817  rgap = wgap(windex)
818  ENDIF
819  gap = min( lgap, rgap )
820  IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
821 * The eigenvector support can become wrong
822 * because significant entries could be cut off due to a
823 * large GAPTOL parameter in LAR1V. Prevent this.
824  gaptol = zero
825  ELSE
826  gaptol = gap * eps
827  ENDIF
828  isupmn = in
829  isupmx = 1
830 * Update WGAP so that it holds the minimum gap
831 * to the left or the right. This is crucial in the
832 * case where bisection is used to ensure that the
833 * eigenvalue is refined up to the required precision.
834 * The correct value is restored afterwards.
835  savgap = wgap(windex)
836  wgap(windex) = gap
837 * We want to use the Rayleigh Quotient Correction
838 * as often as possible since it converges quadratically
839 * when we are close enough to the desired eigenvalue.
840 * However, the Rayleigh Quotient can have the wrong sign
841 * and lead us away from the desired eigenvalue. In this
842 * case, the best we can do is to use bisection.
843  usedbs = .false.
844  usedrq = .false.
845 * Bisection is initially turned off unless it is forced
846  needbs = .NOT.tryrqc
847  120 continue
848 * Check if bisection should be used to refine eigenvalue
849  IF(needbs) THEN
850 * Take the bisection as new iterate
851  usedbs = .true.
852  itmp1 = iwork( iindr+windex )
853  offset = indexw( wbegin ) - 1
854  CALL slarrb( in, d(ibegin),
855  $ work(indlld+ibegin-1),indeig,indeig,
856  $ zero, two*eps, offset,
857  $ work(wbegin),wgap(wbegin),
858  $ werr(wbegin),work( indwrk ),
859  $ iwork( iindwk ), pivmin, spdiam,
860  $ itmp1, iinfo )
861  IF( iinfo.NE.0 ) THEN
862  info = -3
863  return
864  ENDIF
865  lambda = work( windex )
866 * Reset twist index from inaccurate LAMBDA to
867 * force computation of true MINGMA
868  iwork( iindr+windex ) = 0
869  ENDIF
870 * Given LAMBDA, compute the eigenvector.
871  CALL clar1v( in, 1, in, lambda, d( ibegin ),
872  $ l( ibegin ), work(indld+ibegin-1),
873  $ work(indlld+ibegin-1),
874  $ pivmin, gaptol, z( ibegin, windex ),
875  $ .NOT.usedbs, negcnt, ztz, mingma,
876  $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
877  $ nrminv, resid, rqcorr, work( indwrk ) )
878  IF(iter .EQ. 0) THEN
879  bstres = resid
880  bstw = lambda
881  elseif(resid.LT.bstres) THEN
882  bstres = resid
883  bstw = lambda
884  ENDIF
885  isupmn = min(isupmn,isuppz( 2*windex-1 ))
886  isupmx = max(isupmx,isuppz( 2*windex ))
887  iter = iter + 1
888 
889 * sin alpha <= |resid|/gap
890 * Note that both the residual and the gap are
891 * proportional to the matrix, so ||T|| doesn't play
892 * a role in the quotient
893 
894 *
895 * Convergence test for Rayleigh-Quotient iteration
896 * (omitted when Bisection has been used)
897 *
898  IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
899  $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
900  $ THEN
901 * We need to check that the RQCORR update doesn't
902 * move the eigenvalue away from the desired one and
903 * towards a neighbor. -> protection with bisection
904  IF(indeig.LE.negcnt) THEN
905 * The wanted eigenvalue lies to the left
906  sgndef = -one
907  ELSE
908 * The wanted eigenvalue lies to the right
909  sgndef = one
910  ENDIF
911 * We only use the RQCORR if it improves the
912 * the iterate reasonably.
913  IF( ( rqcorr*sgndef.GE.zero )
914  $ .AND.( lambda + rqcorr.LE. right)
915  $ .AND.( lambda + rqcorr.GE. left)
916  $ ) THEN
917  usedrq = .true.
918 * Store new midpoint of bisection interval in WORK
919  IF(sgndef.EQ.one) THEN
920 * The current LAMBDA is on the left of the true
921 * eigenvalue
922  left = lambda
923 * We prefer to assume that the error estimate
924 * is correct. We could make the interval not
925 * as a bracket but to be modified if the RQCORR
926 * chooses to. In this case, the RIGHT side should
927 * be modified as follows:
928 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
929  ELSE
930 * The current LAMBDA is on the right of the true
931 * eigenvalue
932  right = lambda
933 * See comment about assuming the error estimate is
934 * correct above.
935 * LEFT = MIN(LEFT, LAMBDA + RQCORR)
936  ENDIF
937  work( windex ) =
938  $ half * (right + left)
939 * Take RQCORR since it has the correct sign and
940 * improves the iterate reasonably
941  lambda = lambda + rqcorr
942 * Update width of error interval
943  werr( windex ) =
944  $ half * (right-left)
945  ELSE
946  needbs = .true.
947  ENDIF
948  IF(right-left.LT.rqtol*abs(lambda)) THEN
949 * The eigenvalue is computed to bisection accuracy
950 * compute eigenvector and stop
951  usedbs = .true.
952  goto 120
953  elseif( iter.LT.maxitr ) THEN
954  goto 120
955  elseif( iter.EQ.maxitr ) THEN
956  needbs = .true.
957  goto 120
958  ELSE
959  info = 5
960  return
961  END IF
962  ELSE
963  stp2ii = .false.
964  IF(usedrq .AND. usedbs .AND.
965  $ bstres.LE.resid) THEN
966  lambda = bstw
967  stp2ii = .true.
968  ENDIF
969  IF (stp2ii) THEN
970 * improve error angle by second step
971  CALL clar1v( in, 1, in, lambda,
972  $ d( ibegin ), l( ibegin ),
973  $ work(indld+ibegin-1),
974  $ work(indlld+ibegin-1),
975  $ pivmin, gaptol, z( ibegin, windex ),
976  $ .NOT.usedbs, negcnt, ztz, mingma,
977  $ iwork( iindr+windex ),
978  $ isuppz( 2*windex-1 ),
979  $ nrminv, resid, rqcorr, work( indwrk ) )
980  ENDIF
981  work( windex ) = lambda
982  END IF
983 *
984 * Compute FP-vector support w.r.t. whole matrix
985 *
986  isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
987  isuppz( 2*windex ) = isuppz( 2*windex )+oldien
988  zfrom = isuppz( 2*windex-1 )
989  zto = isuppz( 2*windex )
990  isupmn = isupmn + oldien
991  isupmx = isupmx + oldien
992 * Ensure vector is ok if support in the RQI has changed
993  IF(isupmn.LT.zfrom) THEN
994  DO 122 ii = isupmn,zfrom-1
995  z( ii, windex ) = zero
996  122 continue
997  ENDIF
998  IF(isupmx.GT.zto) THEN
999  DO 123 ii = zto+1,isupmx
1000  z( ii, windex ) = zero
1001  123 continue
1002  ENDIF
1003  CALL csscal( zto-zfrom+1, nrminv,
1004  $ z( zfrom, windex ), 1 )
1005  125 continue
1006 * Update W
1007  w( windex ) = lambda+sigma
1008 * Recompute the gaps on the left and right
1009 * But only allow them to become larger and not
1010 * smaller (which can only happen through "bad"
1011 * cancellation and doesn't reflect the theory
1012 * where the initial gaps are underestimated due
1013 * to WERR being too crude.)
1014  IF(.NOT.eskip) THEN
1015  IF( k.GT.1) THEN
1016  wgap( windmn ) = max( wgap(windmn),
1017  $ w(windex)-werr(windex)
1018  $ - w(windmn)-werr(windmn) )
1019  ENDIF
1020  IF( windex.LT.wend ) THEN
1021  wgap( windex ) = max( savgap,
1022  $ w( windpl )-werr( windpl )
1023  $ - w( windex )-werr( windex) )
1024  ENDIF
1025  ENDIF
1026  idone = idone + 1
1027  ENDIF
1028 * here ends the code for the current child
1029 *
1030  139 continue
1031 * Proceed to any remaining child nodes
1032  newfst = j + 1
1033  140 continue
1034  150 continue
1035  ndepth = ndepth + 1
1036  go to 40
1037  END IF
1038  ibegin = iend + 1
1039  wbegin = wend + 1
1040  170 continue
1041 *
1042 
1043  return
1044 *
1045 * End of CLARRV
1046 *
1047  END
1048