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