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