LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dlarrv.f
1 *> \brief \b DLARRV 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 DLARRV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
22 * ISPLIT, M, DOL, DOU, MINRGP,
23 * RTOL1, RTOL2, W, WERR, WGAP,
24 * IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
25 * WORK, IWORK, INFO )
26 *
27 * .. Scalar Arguments ..
28 * INTEGER DOL, DOU, INFO, LDZ, M, N
29 * DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
33 * $ ISUPPZ( * ), IWORK( * )
34 * DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
35 * $ WGAP( * ), WORK( * )
36 * DOUBLE PRECISION Z( LDZ, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> DLARRV computes the eigenvectors of the tridiagonal matrix
46 *> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
47 *> The input eigenvalues should have been computed by DLARRE.
48 *> \endverbatim
49 *
50 * Arguments:
51 * ==========
52 *
53 *> \param[in] N
54 *> \verbatim
55 *> N is INTEGER
56 *> The order of the matrix. N >= 0.
57 *> \endverbatim
58 *>
59 *> \param[in] VL
60 *> \verbatim
61 *> VL is DOUBLE PRECISION
62 *> 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 DOUBLE PRECISION
70 *> Upper bound of the interval that contains the desired
71 *> eigenvalues. VL < VU.
72 *> Note: VU is currently not used by this implementation of DLARRV, VU is
73 *> passed to DLARRV because it could be used compute gaps on the right end
74 *> of the extremal eigenvalues. However, with not much initial accuracy in
75 *> LAMBDA and VU, the formula can lead to an overestimation of the right gap
76 *> and thus to inadequately early RQI 'convergence'. This is currently
77 *> prevented this by forcing a small right gap. And so it turns out that VU
78 *> is currently not used by this implementation of DLARRV.
79 *> \endverbatim
80 *>
81 *> \param[in,out] D
82 *> \verbatim
83 *> D is DOUBLE PRECISION array, dimension (N)
84 *> On entry, the N diagonal elements of the diagonal matrix D.
85 *> On exit, D may be overwritten.
86 *> \endverbatim
87 *>
88 *> \param[in,out] L
89 *> \verbatim
90 *> L is DOUBLE PRECISION array, dimension (N)
91 *> On entry, the (N-1) subdiagonal elements of the unit
92 *> bidiagonal matrix L are in elements 1 to N-1 of L
93 *> (if the matrix is not split.) At the end of each block
94 *> is stored the corresponding shift as given by DLARRE.
95 *> On exit, L is overwritten.
96 *> \endverbatim
97 *>
98 *> \param[in] PIVMIN
99 *> \verbatim
100 *> PIVMIN is DOUBLE PRECISION
101 *> The minimum pivot allowed in the Sturm sequence.
102 *> \endverbatim
103 *>
104 *> \param[in] ISPLIT
105 *> \verbatim
106 *> ISPLIT is INTEGER array, dimension (N)
107 *> The splitting points, at which T breaks up into blocks.
108 *> The first block consists of rows/columns 1 to
109 *> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
110 *> through ISPLIT( 2 ), etc.
111 *> \endverbatim
112 *>
113 *> \param[in] M
114 *> \verbatim
115 *> M is INTEGER
116 *> The total number of input eigenvalues. 0 <= M <= N.
117 *> \endverbatim
118 *>
119 *> \param[in] DOL
120 *> \verbatim
121 *> DOL is INTEGER
122 *> \endverbatim
123 *>
124 *> \param[in] DOU
125 *> \verbatim
126 *> DOU is INTEGER
127 *> If the user wants to compute only selected eigenvectors from all
128 *> the eigenvalues supplied, he can specify an index range DOL:DOU.
129 *> Or else the setting DOL=1, DOU=M should be applied.
130 *> Note that DOL and DOU refer to the order in which the eigenvalues
131 *> are stored in W.
132 *> If the user wants to compute only selected eigenpairs, then
133 *> the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
134 *> computed eigenvectors. All other columns of Z are set to zero.
135 *> \endverbatim
136 *>
137 *> \param[in] MINRGP
138 *> \verbatim
139 *> MINRGP is DOUBLE PRECISION
140 *> \endverbatim
141 *>
142 *> \param[in] RTOL1
143 *> \verbatim
144 *> RTOL1 is DOUBLE PRECISION
145 *> \endverbatim
146 *>
147 *> \param[in] RTOL2
148 *> \verbatim
149 *> RTOL2 is DOUBLE PRECISION
150 *> Parameters for bisection.
151 *> An interval [LEFT,RIGHT] has converged if
152 *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
153 *> \endverbatim
154 *>
155 *> \param[in,out] W
156 *> \verbatim
157 *> W is DOUBLE PRECISION array, dimension (N)
158 *> The first M elements of W contain the APPROXIMATE eigenvalues for
159 *> which eigenvectors are to be computed. The eigenvalues
160 *> should be grouped by split-off block and ordered from
161 *> smallest to largest within the block ( The output array
162 *> W from DLARRE is expected here ). Furthermore, they are with
163 *> respect to the shift of the corresponding root representation
164 *> for their block. On exit, W holds the eigenvalues of the
165 *> UNshifted matrix.
166 *> \endverbatim
167 *>
168 *> \param[in,out] WERR
169 *> \verbatim
170 *> WERR is DOUBLE PRECISION array, dimension (N)
171 *> The first M elements contain the semiwidth of the uncertainty
172 *> interval of the corresponding eigenvalue in W
173 *> \endverbatim
174 *>
175 *> \param[in,out] WGAP
176 *> \verbatim
177 *> WGAP is DOUBLE PRECISION array, dimension (N)
178 *> The separation from the right neighbor eigenvalue in W.
179 *> \endverbatim
180 *>
181 *> \param[in] IBLOCK
182 *> \verbatim
183 *> IBLOCK is INTEGER array, dimension (N)
184 *> The indices of the blocks (submatrices) associated with the
185 *> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
186 *> W(i) belongs to the first block from the top, =2 if W(i)
187 *> belongs to the second block, etc.
188 *> \endverbatim
189 *>
190 *> \param[in] INDEXW
191 *> \verbatim
192 *> INDEXW is INTEGER array, dimension (N)
193 *> The indices of the eigenvalues within each block (submatrix);
194 *> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
195 *> i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
196 *> \endverbatim
197 *>
198 *> \param[in] GERS
199 *> \verbatim
200 *> GERS is DOUBLE PRECISION array, dimension (2*N)
201 *> The N Gerschgorin intervals (the i-th Gerschgorin interval
202 *> is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
203 *> be computed from the original UNshifted matrix.
204 *> \endverbatim
205 *>
206 *> \param[out] Z
207 *> \verbatim
208 *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
209 *> If INFO = 0, the first M columns of Z contain the
210 *> orthonormal eigenvectors of the matrix T
211 *> corresponding to the input eigenvalues, with the i-th
212 *> column of Z holding the eigenvector associated with W(i).
213 *> Note: the user must ensure that at least max(1,M) columns are
214 *> supplied in the array Z.
215 *> \endverbatim
216 *>
217 *> \param[in] LDZ
218 *> \verbatim
219 *> LDZ is INTEGER
220 *> The leading dimension of the array Z. LDZ >= 1, and if
221 *> JOBZ = 'V', LDZ >= max(1,N).
222 *> \endverbatim
223 *>
224 *> \param[out] ISUPPZ
225 *> \verbatim
226 *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
227 *> The support of the eigenvectors in Z, i.e., the indices
228 *> indicating the nonzero elements in Z. The I-th eigenvector
229 *> is nonzero only in elements ISUPPZ( 2*I-1 ) through
230 *> ISUPPZ( 2*I ).
231 *> \endverbatim
232 *>
233 *> \param[out] WORK
234 *> \verbatim
235 *> WORK is DOUBLE PRECISION array, dimension (12*N)
236 *> \endverbatim
237 *>
238 *> \param[out] IWORK
239 *> \verbatim
240 *> IWORK is INTEGER array, dimension (7*N)
241 *> \endverbatim
242 *>
243 *> \param[out] INFO
244 *> \verbatim
245 *> INFO is INTEGER
246 *> = 0: successful exit
247 *>
248 *> > 0: A problem occurred in DLARRV.
249 *> < 0: One of the called subroutines signaled an internal problem.
250 *> Needs inspection of the corresponding parameter IINFO
251 *> for further information.
252 *>
253 *> =-1: Problem in DLARRB when refining a child's eigenvalues.
254 *> =-2: Problem in DLARRF when computing the RRR of a child.
255 *> When a child is inside a tight cluster, it can be difficult
256 *> to find an RRR. A partial remedy from the user's point of
257 *> view is to make the parameter MINRGP smaller and recompile.
258 *> However, as the orthogonality of the computed vectors is
259 *> proportional to 1/MINRGP, the user should be aware that
260 *> he might be trading in precision when he decreases MINRGP.
261 *> =-3: Problem in DLARRB when refining a single eigenvalue
262 *> after the Rayleigh correction was rejected.
263 *> = 5: The Rayleigh Quotient Iteration failed to converge to
264 *> full accuracy in MAXITR steps.
265 *> \endverbatim
266 *
267 * Authors:
268 * ========
269 *
270 *> \author Univ. of Tennessee
271 *> \author Univ. of California Berkeley
272 *> \author Univ. of Colorado Denver
273 *> \author NAG Ltd.
274 *
275 *> \ingroup larrv
276 *
277 *> \par Contributors:
278 * ==================
279 *>
280 *> Beresford Parlett, University of California, Berkeley, USA \n
281 *> Jim Demmel, University of California, Berkeley, USA \n
282 *> Inderjit Dhillon, University of Texas, Austin, USA \n
283 *> Osni Marques, LBNL/NERSC, USA \n
284 *> Christof Voemel, University of California, Berkeley, USA
285 *
286 * =====================================================================
287  SUBROUTINE dlarrv( N, VL, VU, D, L, PIVMIN,
288  $ ISPLIT, M, DOL, DOU, MINRGP,
289  $ RTOL1, RTOL2, W, WERR, WGAP,
290  $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
291  $ WORK, IWORK, INFO )
292 *
293 * -- LAPACK auxiliary routine --
294 * -- LAPACK is a software package provided by Univ. of Tennessee, --
295 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
296 *
297 * .. Scalar Arguments ..
298  INTEGER DOL, DOU, INFO, LDZ, M, N
299  DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
300 * ..
301 * .. Array Arguments ..
302  INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
303  $ isuppz( * ), iwork( * )
304  DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
305  $ wgap( * ), work( * )
306  DOUBLE PRECISION Z( ldz, * )
307 * ..
308 *
309 * =====================================================================
310 *
311 * .. Parameters ..
312  INTEGER MAXITR
313  parameter( maxitr = 10 )
314  DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, HALF
315  parameter( zero = 0.0d0, one = 1.0d0,
316  $ two = 2.0d0, three = 3.0d0,
317  $ four = 4.0d0, half = 0.5d0)
318 * ..
319 * .. Local Scalars ..
320  LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
321  INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
322  $ iindc2, iindr, iindwk, iinfo, im, in, indeig,
323  $ indld, indlld, indwrk, isupmn, isupmx, iter,
324  $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
325  $ ndepth, negcnt, newcls, newfst, newftt, newlst,
326  $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
327  $ oldncl, p, parity, q, wbegin, wend, windex,
328  $ windmn, windpl, zfrom, zto, zusedl, zusedu,
329  $ zusedw
330  DOUBLE PRECISION 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  DOUBLE PRECISION DLAMCH
337  EXTERNAL dlamch
338 * ..
339 * .. External Subroutines ..
340  EXTERNAL dcopy, dlar1v, dlarrb, dlarrf, dlaset,
341  $ dscal
342 * ..
343 * .. Intrinsic Functions ..
344  INTRINSIC abs, dble, max, min
345 * ..
346 * .. Executable Statements ..
347 * ..
348 
349  info = 0
350 *
351 * Quick return if possible
352 *
353  IF( (n.LE.0).OR.(m.LE.0) ) THEN
354  RETURN
355  END IF
356 *
357 * The first N entries of WORK are reserved for the eigenvalues
358  indld = n+1
359  indlld= 2*n+1
360  indwrk= 3*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 dlaset( 'Full', n, zusedw, zero, zero,
396  $ z(1,zusedl), ldz )
397 
398  eps = dlamch( '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 ) = one
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 dcopy( 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  CALL dcopy( in, z( ibegin, j ), 1, d( ibegin ), 1 )
557  CALL dcopy( in-1, z( ibegin, j+1 ), 1, l( ibegin ),
558  $ 1 )
559  sigma = z( iend, j+1 )
560 
561 * Set the corresponding entries in Z to zero
562  CALL dlaset( 'Full', in, 2, zero, zero,
563  $ z( ibegin, j), ldz )
564  END IF
565 
566 * Compute DL and DLL of current RRR
567  DO 50 j = ibegin, iend-1
568  tmp = d( j )*l( j )
569  work( indld-1+j ) = tmp
570  work( indlld-1+j ) = tmp*l( j )
571  50 CONTINUE
572 
573  IF( ndepth.GT.0 ) THEN
574 * P and Q are index of the first and last eigenvalue to compute
575 * within the current block
576  p = indexw( wbegin-1+oldfst )
577  q = indexw( wbegin-1+oldlst )
578 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
579 * through the Q-OFFSET elements of these arrays are to be used.
580 * OFFSET = P-OLDFST
581  offset = indexw( wbegin ) - 1
582 * perform limited bisection (if necessary) to get approximate
583 * eigenvalues to the precision needed.
584  CALL dlarrb( in, d( ibegin ),
585  $ work(indlld+ibegin-1),
586  $ p, q, rtol1, rtol2, offset,
587  $ work(wbegin),wgap(wbegin),werr(wbegin),
588  $ work( indwrk ), iwork( iindwk ),
589  $ pivmin, spdiam, in, iinfo )
590  IF( iinfo.NE.0 ) THEN
591  info = -1
592  RETURN
593  ENDIF
594 * We also recompute the extremal gaps. W holds all eigenvalues
595 * of the unshifted matrix and must be used for computation
596 * of WGAP, the entries of WORK might stem from RRRs with
597 * different shifts. The gaps from WBEGIN-1+OLDFST to
598 * WBEGIN-1+OLDLST are correctly computed in DLARRB.
599 * However, we only allow the gaps to become greater since
600 * this is what should happen when we decrease WERR
601  IF( oldfst.GT.1) THEN
602  wgap( wbegin+oldfst-2 ) =
603  $ max(wgap(wbegin+oldfst-2),
604  $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
605  $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
606  ENDIF
607  IF( wbegin + oldlst -1 .LT. wend ) THEN
608  wgap( wbegin+oldlst-1 ) =
609  $ max(wgap(wbegin+oldlst-1),
610  $ w(wbegin+oldlst)-werr(wbegin+oldlst)
611  $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
612  ENDIF
613 * Each time the eigenvalues in WORK get refined, we store
614 * the newly found approximation with all shifts applied in W
615  DO 53 j=oldfst,oldlst
616  w(wbegin+j-1) = work(wbegin+j-1)+sigma
617  53 CONTINUE
618  END IF
619 
620 * Process the current node.
621  newfst = oldfst
622  DO 140 j = oldfst, oldlst
623  IF( j.EQ.oldlst ) THEN
624 * we are at the right end of the cluster, this is also the
625 * boundary of the child cluster
626  newlst = j
627  ELSE IF ( wgap( wbegin + j -1).GE.
628  $ minrgp* abs( work(wbegin + j -1) ) ) THEN
629 * the right relative gap is big enough, the child cluster
630 * (NEWFST,..,NEWLST) is well separated from the following
631  newlst = j
632  ELSE
633 * inside a child cluster, the relative gap is not
634 * big enough.
635  GOTO 140
636  END IF
637 
638 * Compute size of child cluster found
639  newsiz = newlst - newfst + 1
640 
641 * NEWFTT is the place in Z where the new RRR or the computed
642 * eigenvector is to be stored
643  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
644 * Store representation at location of the leftmost evalue
645 * of the cluster
646  newftt = wbegin + newfst - 1
647  ELSE
648  IF(wbegin+newfst-1.LT.dol) THEN
649 * Store representation at the left end of Z array
650  newftt = dol - 1
651  ELSEIF(wbegin+newfst-1.GT.dou) THEN
652 * Store representation at the right end of Z array
653  newftt = dou
654  ELSE
655  newftt = wbegin + newfst - 1
656  ENDIF
657  ENDIF
658 
659  IF( newsiz.GT.1) THEN
660 *
661 * Current child is not a singleton but a cluster.
662 * Compute and store new representation of child.
663 *
664 *
665 * Compute left and right cluster gap.
666 *
667 * LGAP and RGAP are not computed from WORK because
668 * the eigenvalue approximations may stem from RRRs
669 * different shifts. However, W hold all eigenvalues
670 * of the unshifted matrix. Still, the entries in WGAP
671 * have to be computed from WORK since the entries
672 * in W might be of the same order so that gaps are not
673 * exhibited correctly for very close eigenvalues.
674  IF( newfst.EQ.1 ) THEN
675  lgap = max( zero,
676  $ w(wbegin)-werr(wbegin) - vl )
677  ELSE
678  lgap = wgap( wbegin+newfst-2 )
679  ENDIF
680  rgap = wgap( wbegin+newlst-1 )
681 *
682 * Compute left- and rightmost eigenvalue of child
683 * to high precision in order to shift as close
684 * as possible and obtain as large relative gaps
685 * as possible
686 *
687  DO 55 k =1,2
688  IF(k.EQ.1) THEN
689  p = indexw( wbegin-1+newfst )
690  ELSE
691  p = indexw( wbegin-1+newlst )
692  ENDIF
693  offset = indexw( wbegin ) - 1
694  CALL dlarrb( in, d(ibegin),
695  $ work( indlld+ibegin-1 ),p,p,
696  $ rqtol, rqtol, offset,
697  $ work(wbegin),wgap(wbegin),
698  $ werr(wbegin),work( indwrk ),
699  $ iwork( iindwk ), pivmin, spdiam,
700  $ in, iinfo )
701  55 CONTINUE
702 *
703  IF((wbegin+newlst-1.LT.dol).OR.
704  $ (wbegin+newfst-1.GT.dou)) THEN
705 * if the cluster contains no desired eigenvalues
706 * skip the computation of that branch of the rep. tree
707 *
708 * We could skip before the refinement of the extremal
709 * eigenvalues of the child, but then the representation
710 * tree could be different from the one when nothing is
711 * skipped. For this reason we skip at this place.
712  idone = idone + newlst - newfst + 1
713  GOTO 139
714  ENDIF
715 *
716 * Compute RRR of child cluster.
717 * Note that the new RRR is stored in Z
718 *
719 * DLARRF needs LWORK = 2*N
720  CALL dlarrf( in, d( ibegin ), l( ibegin ),
721  $ work(indld+ibegin-1),
722  $ newfst, newlst, work(wbegin),
723  $ wgap(wbegin), werr(wbegin),
724  $ spdiam, lgap, rgap, pivmin, tau,
725  $ z(ibegin, newftt),z(ibegin, newftt+1),
726  $ work( indwrk ), iinfo )
727  IF( iinfo.EQ.0 ) THEN
728 * a new RRR for the cluster was found by DLARRF
729 * update shift and store it
730  ssigma = sigma + tau
731  z( iend, newftt+1 ) = ssigma
732 * WORK() are the midpoints and WERR() the semi-width
733 * Note that the entries in W are unchanged.
734  DO 116 k = newfst, newlst
735  fudge =
736  $ three*eps*abs(work(wbegin+k-1))
737  work( wbegin + k - 1 ) =
738  $ work( wbegin + k - 1) - tau
739  fudge = fudge +
740  $ four*eps*abs(work(wbegin+k-1))
741 * Fudge errors
742  werr( wbegin + k - 1 ) =
743  $ werr( wbegin + k - 1 ) + fudge
744 * Gaps are not fudged. Provided that WERR is small
745 * when eigenvalues are close, a zero gap indicates
746 * that a new representation is needed for resolving
747 * the cluster. A fudge could lead to a wrong decision
748 * of judging eigenvalues 'separated' which in
749 * reality are not. This could have a negative impact
750 * on the orthogonality of the computed eigenvectors.
751  116 CONTINUE
752 
753  nclus = nclus + 1
754  k = newcls + 2*nclus
755  iwork( k-1 ) = newfst
756  iwork( k ) = newlst
757  ELSE
758  info = -2
759  RETURN
760  ENDIF
761  ELSE
762 *
763 * Compute eigenvector of singleton
764 *
765  iter = 0
766 *
767  tol = four * log(dble(in)) * eps
768 *
769  k = newfst
770  windex = wbegin + k - 1
771  windmn = max(windex - 1,1)
772  windpl = min(windex + 1,m)
773  lambda = work( windex )
774  done = done + 1
775 * Check if eigenvector computation is to be skipped
776  IF((windex.LT.dol).OR.
777  $ (windex.GT.dou)) THEN
778  eskip = .true.
779  GOTO 125
780  ELSE
781  eskip = .false.
782  ENDIF
783  left = work( windex ) - werr( windex )
784  right = work( windex ) + werr( windex )
785  indeig = indexw( windex )
786 * Note that since we compute the eigenpairs for a child,
787 * all eigenvalue approximations are w.r.t the same shift.
788 * In this case, the entries in WORK should be used for
789 * computing the gaps since they exhibit even very small
790 * differences in the eigenvalues, as opposed to the
791 * entries in W which might "look" the same.
792 
793  IF( k .EQ. 1) THEN
794 * In the case RANGE='I' and with not much initial
795 * accuracy in LAMBDA and VL, the formula
796 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
797 * can lead to an overestimation of the left gap and
798 * thus to inadequately early RQI 'convergence'.
799 * Prevent this by forcing a small left gap.
800  lgap = eps*max(abs(left),abs(right))
801  ELSE
802  lgap = wgap(windmn)
803  ENDIF
804  IF( k .EQ. im) THEN
805 * In the case RANGE='I' and with not much initial
806 * accuracy in LAMBDA and VU, the formula
807 * can lead to an overestimation of the right gap and
808 * thus to inadequately early RQI 'convergence'.
809 * Prevent this by forcing a small right gap.
810  rgap = eps*max(abs(left),abs(right))
811  ELSE
812  rgap = wgap(windex)
813  ENDIF
814  gap = min( lgap, rgap )
815  IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
816 * The eigenvector support can become wrong
817 * because significant entries could be cut off due to a
818 * large GAPTOL parameter in LAR1V. Prevent this.
819  gaptol = zero
820  ELSE
821  gaptol = gap * eps
822  ENDIF
823  isupmn = in
824  isupmx = 1
825 * Update WGAP so that it holds the minimum gap
826 * to the left or the right. This is crucial in the
827 * case where bisection is used to ensure that the
828 * eigenvalue is refined up to the required precision.
829 * The correct value is restored afterwards.
830  savgap = wgap(windex)
831  wgap(windex) = gap
832 * We want to use the Rayleigh Quotient Correction
833 * as often as possible since it converges quadratically
834 * when we are close enough to the desired eigenvalue.
835 * However, the Rayleigh Quotient can have the wrong sign
836 * and lead us away from the desired eigenvalue. In this
837 * case, the best we can do is to use bisection.
838  usedbs = .false.
839  usedrq = .false.
840 * Bisection is initially turned off unless it is forced
841  needbs = .NOT.tryrqc
842  120 CONTINUE
843 * Check if bisection should be used to refine eigenvalue
844  IF(needbs) THEN
845 * Take the bisection as new iterate
846  usedbs = .true.
847  itmp1 = iwork( iindr+windex )
848  offset = indexw( wbegin ) - 1
849  CALL dlarrb( in, d(ibegin),
850  $ work(indlld+ibegin-1),indeig,indeig,
851  $ zero, two*eps, offset,
852  $ work(wbegin),wgap(wbegin),
853  $ werr(wbegin),work( indwrk ),
854  $ iwork( iindwk ), pivmin, spdiam,
855  $ itmp1, iinfo )
856  IF( iinfo.NE.0 ) THEN
857  info = -3
858  RETURN
859  ENDIF
860  lambda = work( windex )
861 * Reset twist index from inaccurate LAMBDA to
862 * force computation of true MINGMA
863  iwork( iindr+windex ) = 0
864  ENDIF
865 * Given LAMBDA, compute the eigenvector.
866  CALL dlar1v( in, 1, in, lambda, d( ibegin ),
867  $ l( ibegin ), work(indld+ibegin-1),
868  $ work(indlld+ibegin-1),
869  $ pivmin, gaptol, z( ibegin, windex ),
870  $ .NOT.usedbs, negcnt, ztz, mingma,
871  $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
872  $ nrminv, resid, rqcorr, work( indwrk ) )
873  IF(iter .EQ. 0) THEN
874  bstres = resid
875  bstw = lambda
876  ELSEIF(resid.LT.bstres) THEN
877  bstres = resid
878  bstw = lambda
879  ENDIF
880  isupmn = min(isupmn,isuppz( 2*windex-1 ))
881  isupmx = max(isupmx,isuppz( 2*windex ))
882  iter = iter + 1
883 
884 * sin alpha <= |resid|/gap
885 * Note that both the residual and the gap are
886 * proportional to the matrix, so ||T|| doesn't play
887 * a role in the quotient
888 
889 *
890 * Convergence test for Rayleigh-Quotient iteration
891 * (omitted when Bisection has been used)
892 *
893  IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
894  $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
895  $ THEN
896 * We need to check that the RQCORR update doesn't
897 * move the eigenvalue away from the desired one and
898 * towards a neighbor. -> protection with bisection
899  IF(indeig.LE.negcnt) THEN
900 * The wanted eigenvalue lies to the left
901  sgndef = -one
902  ELSE
903 * The wanted eigenvalue lies to the right
904  sgndef = one
905  ENDIF
906 * We only use the RQCORR if it improves the
907 * the iterate reasonably.
908  IF( ( rqcorr*sgndef.GE.zero )
909  $ .AND.( lambda + rqcorr.LE. right)
910  $ .AND.( lambda + rqcorr.GE. left)
911  $ ) THEN
912  usedrq = .true.
913 * Store new midpoint of bisection interval in WORK
914  IF(sgndef.EQ.one) THEN
915 * The current LAMBDA is on the left of the true
916 * eigenvalue
917  left = lambda
918 * We prefer to assume that the error estimate
919 * is correct. We could make the interval not
920 * as a bracket but to be modified if the RQCORR
921 * chooses to. In this case, the RIGHT side should
922 * be modified as follows:
923 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
924  ELSE
925 * The current LAMBDA is on the right of the true
926 * eigenvalue
927  right = lambda
928 * See comment about assuming the error estimate is
929 * correct above.
930 * LEFT = MIN(LEFT, LAMBDA + RQCORR)
931  ENDIF
932  work( windex ) =
933  $ half * (right + left)
934 * Take RQCORR since it has the correct sign and
935 * improves the iterate reasonably
936  lambda = lambda + rqcorr
937 * Update width of error interval
938  werr( windex ) =
939  $ half * (right-left)
940  ELSE
941  needbs = .true.
942  ENDIF
943  IF(right-left.LT.rqtol*abs(lambda)) THEN
944 * The eigenvalue is computed to bisection accuracy
945 * compute eigenvector and stop
946  usedbs = .true.
947  GOTO 120
948  ELSEIF( iter.LT.maxitr ) THEN
949  GOTO 120
950  ELSEIF( iter.EQ.maxitr ) THEN
951  needbs = .true.
952  GOTO 120
953  ELSE
954  info = 5
955  RETURN
956  END IF
957  ELSE
958  stp2ii = .false.
959  IF(usedrq .AND. usedbs .AND.
960  $ bstres.LE.resid) THEN
961  lambda = bstw
962  stp2ii = .true.
963  ENDIF
964  IF (stp2ii) THEN
965 * improve error angle by second step
966  CALL dlar1v( in, 1, in, lambda,
967  $ d( ibegin ), l( ibegin ),
968  $ work(indld+ibegin-1),
969  $ work(indlld+ibegin-1),
970  $ pivmin, gaptol, z( ibegin, windex ),
971  $ .NOT.usedbs, negcnt, ztz, mingma,
972  $ iwork( iindr+windex ),
973  $ isuppz( 2*windex-1 ),
974  $ nrminv, resid, rqcorr, work( indwrk ) )
975  ENDIF
976  work( windex ) = lambda
977  END IF
978 *
979 * Compute FP-vector support w.r.t. whole matrix
980 *
981  isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
982  isuppz( 2*windex ) = isuppz( 2*windex )+oldien
983  zfrom = isuppz( 2*windex-1 )
984  zto = isuppz( 2*windex )
985  isupmn = isupmn + oldien
986  isupmx = isupmx + oldien
987 * Ensure vector is ok if support in the RQI has changed
988  IF(isupmn.LT.zfrom) THEN
989  DO 122 ii = isupmn,zfrom-1
990  z( ii, windex ) = zero
991  122 CONTINUE
992  ENDIF
993  IF(isupmx.GT.zto) THEN
994  DO 123 ii = zto+1,isupmx
995  z( ii, windex ) = zero
996  123 CONTINUE
997  ENDIF
998  CALL dscal( zto-zfrom+1, nrminv,
999  $ z( zfrom, windex ), 1 )
1000  125 CONTINUE
1001 * Update W
1002  w( windex ) = lambda+sigma
1003 * Recompute the gaps on the left and right
1004 * But only allow them to become larger and not
1005 * smaller (which can only happen through "bad"
1006 * cancellation and doesn't reflect the theory
1007 * where the initial gaps are underestimated due
1008 * to WERR being too crude.)
1009  IF(.NOT.eskip) THEN
1010  IF( k.GT.1) THEN
1011  wgap( windmn ) = max( wgap(windmn),
1012  $ w(windex)-werr(windex)
1013  $ - w(windmn)-werr(windmn) )
1014  ENDIF
1015  IF( windex.LT.wend ) THEN
1016  wgap( windex ) = max( savgap,
1017  $ w( windpl )-werr( windpl )
1018  $ - w( windex )-werr( windex) )
1019  ENDIF
1020  ENDIF
1021  idone = idone + 1
1022  ENDIF
1023 * here ends the code for the current child
1024 *
1025  139 CONTINUE
1026 * Proceed to any remaining child nodes
1027  newfst = j + 1
1028  140 CONTINUE
1029  150 CONTINUE
1030  ndepth = ndepth + 1
1031  GO TO 40
1032  END IF
1033  ibegin = iend + 1
1034  wbegin = wend + 1
1035  170 CONTINUE
1036 *
1037 
1038  RETURN
1039 *
1040 * End of DLARRV
1041 *
1042  END
dlar1v
subroutine dlar1v(N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)
DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the...
Definition: dlar1v.f:230
dlarrf
subroutine dlarrf(N, D, L, LD, CLSTRT, CLEND, W, WGAP, WERR, SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, DPLUS, LPLUS, WORK, INFO)
DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is rela...
Definition: dlarrf.f:193
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
dlarrb
subroutine dlarrb(N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
DLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition: dlarrb.f:196
dlarrv
subroutine dlarrv(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)
DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition: dlarrv.f:292
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:110
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82