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