LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Files Functions Typedefs Macros
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 
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:107
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:191
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:195
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:52
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:280
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:229
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:53