LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zlarrv.f
Go to the documentation of this file.
1 *> \brief \b ZLARRV 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 ZLARRV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarrv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarrv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarrv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLARRV( 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 * COMPLEX*16 Z( LDZ, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> ZLARRV 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 COMPLEX*16 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 ZLARRV.
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 complex16OTHERauxiliary
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 zlarrv( 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  COMPLEX*16 Z( ldz, * )
304 * ..
305 *
306 * =====================================================================
307 *
308 * .. Parameters ..
309  INTEGER MAXITR
310  parameter ( maxitr = 10 )
311  COMPLEX*16 CZERO
312  parameter ( czero = ( 0.0d0, 0.0d0 ) )
313  DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, HALF
314  parameter ( zero = 0.0d0, one = 1.0d0,
315  $ two = 2.0d0, three = 3.0d0,
316  $ four = 4.0d0, half = 0.5d0)
317 * ..
318 * .. Local Scalars ..
319  LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
320  INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
321  $ iindc2, iindr, iindwk, iinfo, im, in, indeig,
322  $ indld, indlld, indwrk, isupmn, isupmx, iter,
323  $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
324  $ ndepth, negcnt, newcls, newfst, newftt, newlst,
325  $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
326  $ oldncl, p, parity, q, wbegin, wend, windex,
327  $ windmn, windpl, zfrom, zto, zusedl, zusedu,
328  $ zusedw
329  INTEGER INDIN1, INDIN2
330  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, dlarrb, dlarrf, zdscal, zlar1v,
341  $ zlaset
342 * ..
343 * .. Intrinsic Functions ..
344  INTRINSIC abs, dble, max, min
345  INTRINSIC dcmplx
346 * ..
347 * .. Executable Statements ..
348 * ..
349 
350  info = 0
351 * The first N entries of WORK are reserved for the eigenvalues
352  indld = n+1
353  indlld= 2*n+1
354  indin1 = 3*n + 1
355  indin2 = 4*n + 1
356  indwrk = 5*n + 1
357  minwsize = 12 * n
358 
359  DO 5 i= 1,minwsize
360  work( i ) = zero
361  5 CONTINUE
362 
363 * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
364 * factorization used to compute the FP vector
365  iindr = 0
366 * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
367 * layer and the one above.
368  iindc1 = n
369  iindc2 = 2*n
370  iindwk = 3*n + 1
371 
372  miniwsize = 7 * n
373  DO 10 i= 1,miniwsize
374  iwork( i ) = 0
375  10 CONTINUE
376 
377  zusedl = 1
378  IF(dol.GT.1) THEN
379 * Set lower bound for use of Z
380  zusedl = dol-1
381  ENDIF
382  zusedu = m
383  IF(dou.LT.m) THEN
384 * Set lower bound for use of Z
385  zusedu = dou+1
386  ENDIF
387 * The width of the part of Z that is used
388  zusedw = zusedu - zusedl + 1
389 
390 
391  CALL zlaset( 'Full', n, zusedw, czero, czero,
392  $ z(1,zusedl), ldz )
393 
394  eps = dlamch( 'Precision' )
395  rqtol = two * eps
396 *
397 * Set expert flags for standard code.
398  tryrqc = .true.
399 
400  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
401  ELSE
402 * Only selected eigenpairs are computed. Since the other evalues
403 * are not refined by RQ iteration, bisection has to compute to full
404 * accuracy.
405  rtol1 = four * eps
406  rtol2 = four * eps
407  ENDIF
408 
409 * The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
410 * desired eigenvalues. The support of the nonzero eigenvector
411 * entries is contained in the interval IBEGIN:IEND.
412 * Remark that if k eigenpairs are desired, then the eigenvectors
413 * are stored in k contiguous columns of Z.
414 
415 * DONE is the number of eigenvectors already computed
416  done = 0
417  ibegin = 1
418  wbegin = 1
419  DO 170 jblk = 1, iblock( m )
420  iend = isplit( jblk )
421  sigma = l( iend )
422 * Find the eigenvectors of the submatrix indexed IBEGIN
423 * through IEND.
424  wend = wbegin - 1
425  15 CONTINUE
426  IF( wend.LT.m ) THEN
427  IF( iblock( wend+1 ).EQ.jblk ) THEN
428  wend = wend + 1
429  GO TO 15
430  END IF
431  END IF
432  IF( wend.LT.wbegin ) THEN
433  ibegin = iend + 1
434  GO TO 170
435  ELSEIF( (wend.LT.dol).OR.(wbegin.GT.dou) ) THEN
436  ibegin = iend + 1
437  wbegin = wend + 1
438  GO TO 170
439  END IF
440 
441 * Find local spectral diameter of the block
442  gl = gers( 2*ibegin-1 )
443  gu = gers( 2*ibegin )
444  DO 20 i = ibegin+1 , iend
445  gl = min( gers( 2*i-1 ), gl )
446  gu = max( gers( 2*i ), gu )
447  20 CONTINUE
448  spdiam = gu - gl
449 
450 * OLDIEN is the last index of the previous block
451  oldien = ibegin - 1
452 * Calculate the size of the current block
453  in = iend - ibegin + 1
454 * The number of eigenvalues in the current block
455  im = wend - wbegin + 1
456 
457 * This is for a 1x1 block
458  IF( ibegin.EQ.iend ) THEN
459  done = done+1
460  z( ibegin, wbegin ) = dcmplx( one, zero )
461  isuppz( 2*wbegin-1 ) = ibegin
462  isuppz( 2*wbegin ) = ibegin
463  w( wbegin ) = w( wbegin ) + sigma
464  work( wbegin ) = w( wbegin )
465  ibegin = iend + 1
466  wbegin = wbegin + 1
467  GO TO 170
468  END IF
469 
470 * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
471 * Note that these can be approximations, in this case, the corresp.
472 * entries of WERR give the size of the uncertainty interval.
473 * The eigenvalue approximations will be refined when necessary as
474 * high relative accuracy is required for the computation of the
475 * corresponding eigenvectors.
476  CALL dcopy( im, w( wbegin ), 1,
477  $ work( wbegin ), 1 )
478 
479 * We store in W the eigenvalue approximations w.r.t. the original
480 * matrix T.
481  DO 30 i=1,im
482  w(wbegin+i-1) = w(wbegin+i-1)+sigma
483  30 CONTINUE
484 
485 
486 * NDEPTH is the current depth of the representation tree
487  ndepth = 0
488 * PARITY is either 1 or 0
489  parity = 1
490 * NCLUS is the number of clusters for the next level of the
491 * representation tree, we start with NCLUS = 1 for the root
492  nclus = 1
493  iwork( iindc1+1 ) = 1
494  iwork( iindc1+2 ) = im
495 
496 * IDONE is the number of eigenvectors already computed in the current
497 * block
498  idone = 0
499 * loop while( IDONE.LT.IM )
500 * generate the representation tree for the current block and
501 * compute the eigenvectors
502  40 CONTINUE
503  IF( idone.LT.im ) THEN
504 * This is a crude protection against infinitely deep trees
505  IF( ndepth.GT.m ) THEN
506  info = -2
507  RETURN
508  ENDIF
509 * breadth first processing of the current level of the representation
510 * tree: OLDNCL = number of clusters on current level
511  oldncl = nclus
512 * reset NCLUS to count the number of child clusters
513  nclus = 0
514 *
515  parity = 1 - parity
516  IF( parity.EQ.0 ) THEN
517  oldcls = iindc1
518  newcls = iindc2
519  ELSE
520  oldcls = iindc2
521  newcls = iindc1
522  END IF
523 * Process the clusters on the current level
524  DO 150 i = 1, oldncl
525  j = oldcls + 2*i
526 * OLDFST, OLDLST = first, last index of current cluster.
527 * cluster indices start with 1 and are relative
528 * to WBEGIN when accessing W, WGAP, WERR, Z
529  oldfst = iwork( j-1 )
530  oldlst = iwork( j )
531  IF( ndepth.GT.0 ) THEN
532 * Retrieve relatively robust representation (RRR) of cluster
533 * that has been computed at the previous level
534 * The RRR is stored in Z and overwritten once the eigenvectors
535 * have been computed or when the cluster is refined
536 
537  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
538 * Get representation from location of the leftmost evalue
539 * of the cluster
540  j = wbegin + oldfst - 1
541  ELSE
542  IF(wbegin+oldfst-1.LT.dol) THEN
543 * Get representation from the left end of Z array
544  j = dol - 1
545  ELSEIF(wbegin+oldfst-1.GT.dou) THEN
546 * Get representation from the right end of Z array
547  j = dou
548  ELSE
549  j = wbegin + oldfst - 1
550  ENDIF
551  ENDIF
552  DO 45 k = 1, in - 1
553  d( ibegin+k-1 ) = dble( z( ibegin+k-1,
554  $ j ) )
555  l( ibegin+k-1 ) = dble( z( ibegin+k-1,
556  $ j+1 ) )
557  45 CONTINUE
558  d( iend ) = dble( z( iend, j ) )
559  sigma = dble( z( iend, j+1 ) )
560 
561 * Set the corresponding entries in Z to zero
562  CALL zlaset( 'Full', in, 2, czero, czero,
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  $ work( indin1 ), work( indin2 ),
726  $ work( indwrk ), iinfo )
727 * In the complex case, DLARRF cannot write
728 * the new RRR directly into Z and needs an intermediate
729 * workspace
730  DO 56 k = 1, in-1
731  z( ibegin+k-1, newftt ) =
732  $ dcmplx( work( indin1+k-1 ), zero )
733  z( ibegin+k-1, newftt+1 ) =
734  $ dcmplx( work( indin2+k-1 ), zero )
735  56 CONTINUE
736  z( iend, newftt ) =
737  $ dcmplx( work( indin1+in-1 ), zero )
738  IF( iinfo.EQ.0 ) THEN
739 * a new RRR for the cluster was found by DLARRF
740 * update shift and store it
741  ssigma = sigma + tau
742  z( iend, newftt+1 ) = dcmplx( ssigma, zero )
743 * WORK() are the midpoints and WERR() the semi-width
744 * Note that the entries in W are unchanged.
745  DO 116 k = newfst, newlst
746  fudge =
747  $ three*eps*abs(work(wbegin+k-1))
748  work( wbegin + k - 1 ) =
749  $ work( wbegin + k - 1) - tau
750  fudge = fudge +
751  $ four*eps*abs(work(wbegin+k-1))
752 * Fudge errors
753  werr( wbegin + k - 1 ) =
754  $ werr( wbegin + k - 1 ) + fudge
755 * Gaps are not fudged. Provided that WERR is small
756 * when eigenvalues are close, a zero gap indicates
757 * that a new representation is needed for resolving
758 * the cluster. A fudge could lead to a wrong decision
759 * of judging eigenvalues 'separated' which in
760 * reality are not. This could have a negative impact
761 * on the orthogonality of the computed eigenvectors.
762  116 CONTINUE
763 
764  nclus = nclus + 1
765  k = newcls + 2*nclus
766  iwork( k-1 ) = newfst
767  iwork( k ) = newlst
768  ELSE
769  info = -2
770  RETURN
771  ENDIF
772  ELSE
773 *
774 * Compute eigenvector of singleton
775 *
776  iter = 0
777 *
778  tol = four * log(dble(in)) * eps
779 *
780  k = newfst
781  windex = wbegin + k - 1
782  windmn = max(windex - 1,1)
783  windpl = min(windex + 1,m)
784  lambda = work( windex )
785  done = done + 1
786 * Check if eigenvector computation is to be skipped
787  IF((windex.LT.dol).OR.
788  $ (windex.GT.dou)) THEN
789  eskip = .true.
790  GOTO 125
791  ELSE
792  eskip = .false.
793  ENDIF
794  left = work( windex ) - werr( windex )
795  right = work( windex ) + werr( windex )
796  indeig = indexw( windex )
797 * Note that since we compute the eigenpairs for a child,
798 * all eigenvalue approximations are w.r.t the same shift.
799 * In this case, the entries in WORK should be used for
800 * computing the gaps since they exhibit even very small
801 * differences in the eigenvalues, as opposed to the
802 * entries in W which might "look" the same.
803 
804  IF( k .EQ. 1) THEN
805 * In the case RANGE='I' and with not much initial
806 * accuracy in LAMBDA and VL, the formula
807 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
808 * can lead to an overestimation of the left gap and
809 * thus to inadequately early RQI 'convergence'.
810 * Prevent this by forcing a small left gap.
811  lgap = eps*max(abs(left),abs(right))
812  ELSE
813  lgap = wgap(windmn)
814  ENDIF
815  IF( k .EQ. im) THEN
816 * In the case RANGE='I' and with not much initial
817 * accuracy in LAMBDA and VU, the formula
818 * can lead to an overestimation of the right gap and
819 * thus to inadequately early RQI 'convergence'.
820 * Prevent this by forcing a small right gap.
821  rgap = eps*max(abs(left),abs(right))
822  ELSE
823  rgap = wgap(windex)
824  ENDIF
825  gap = min( lgap, rgap )
826  IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
827 * The eigenvector support can become wrong
828 * because significant entries could be cut off due to a
829 * large GAPTOL parameter in LAR1V. Prevent this.
830  gaptol = zero
831  ELSE
832  gaptol = gap * eps
833  ENDIF
834  isupmn = in
835  isupmx = 1
836 * Update WGAP so that it holds the minimum gap
837 * to the left or the right. This is crucial in the
838 * case where bisection is used to ensure that the
839 * eigenvalue is refined up to the required precision.
840 * The correct value is restored afterwards.
841  savgap = wgap(windex)
842  wgap(windex) = gap
843 * We want to use the Rayleigh Quotient Correction
844 * as often as possible since it converges quadratically
845 * when we are close enough to the desired eigenvalue.
846 * However, the Rayleigh Quotient can have the wrong sign
847 * and lead us away from the desired eigenvalue. In this
848 * case, the best we can do is to use bisection.
849  usedbs = .false.
850  usedrq = .false.
851 * Bisection is initially turned off unless it is forced
852  needbs = .NOT.tryrqc
853  120 CONTINUE
854 * Check if bisection should be used to refine eigenvalue
855  IF(needbs) THEN
856 * Take the bisection as new iterate
857  usedbs = .true.
858  itmp1 = iwork( iindr+windex )
859  offset = indexw( wbegin ) - 1
860  CALL dlarrb( in, d(ibegin),
861  $ work(indlld+ibegin-1),indeig,indeig,
862  $ zero, two*eps, offset,
863  $ work(wbegin),wgap(wbegin),
864  $ werr(wbegin),work( indwrk ),
865  $ iwork( iindwk ), pivmin, spdiam,
866  $ itmp1, iinfo )
867  IF( iinfo.NE.0 ) THEN
868  info = -3
869  RETURN
870  ENDIF
871  lambda = work( windex )
872 * Reset twist index from inaccurate LAMBDA to
873 * force computation of true MINGMA
874  iwork( iindr+windex ) = 0
875  ENDIF
876 * Given LAMBDA, compute the eigenvector.
877  CALL zlar1v( in, 1, in, lambda, d( ibegin ),
878  $ l( ibegin ), work(indld+ibegin-1),
879  $ work(indlld+ibegin-1),
880  $ pivmin, gaptol, z( ibegin, windex ),
881  $ .NOT.usedbs, negcnt, ztz, mingma,
882  $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
883  $ nrminv, resid, rqcorr, work( indwrk ) )
884  IF(iter .EQ. 0) THEN
885  bstres = resid
886  bstw = lambda
887  ELSEIF(resid.LT.bstres) THEN
888  bstres = resid
889  bstw = lambda
890  ENDIF
891  isupmn = min(isupmn,isuppz( 2*windex-1 ))
892  isupmx = max(isupmx,isuppz( 2*windex ))
893  iter = iter + 1
894 
895 * sin alpha <= |resid|/gap
896 * Note that both the residual and the gap are
897 * proportional to the matrix, so ||T|| doesn't play
898 * a role in the quotient
899 
900 *
901 * Convergence test for Rayleigh-Quotient iteration
902 * (omitted when Bisection has been used)
903 *
904  IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
905  $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
906  $ THEN
907 * We need to check that the RQCORR update doesn't
908 * move the eigenvalue away from the desired one and
909 * towards a neighbor. -> protection with bisection
910  IF(indeig.LE.negcnt) THEN
911 * The wanted eigenvalue lies to the left
912  sgndef = -one
913  ELSE
914 * The wanted eigenvalue lies to the right
915  sgndef = one
916  ENDIF
917 * We only use the RQCORR if it improves the
918 * the iterate reasonably.
919  IF( ( rqcorr*sgndef.GE.zero )
920  $ .AND.( lambda + rqcorr.LE. right)
921  $ .AND.( lambda + rqcorr.GE. left)
922  $ ) THEN
923  usedrq = .true.
924 * Store new midpoint of bisection interval in WORK
925  IF(sgndef.EQ.one) THEN
926 * The current LAMBDA is on the left of the true
927 * eigenvalue
928  left = lambda
929 * We prefer to assume that the error estimate
930 * is correct. We could make the interval not
931 * as a bracket but to be modified if the RQCORR
932 * chooses to. In this case, the RIGHT side should
933 * be modified as follows:
934 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
935  ELSE
936 * The current LAMBDA is on the right of the true
937 * eigenvalue
938  right = lambda
939 * See comment about assuming the error estimate is
940 * correct above.
941 * LEFT = MIN(LEFT, LAMBDA + RQCORR)
942  ENDIF
943  work( windex ) =
944  $ half * (right + left)
945 * Take RQCORR since it has the correct sign and
946 * improves the iterate reasonably
947  lambda = lambda + rqcorr
948 * Update width of error interval
949  werr( windex ) =
950  $ half * (right-left)
951  ELSE
952  needbs = .true.
953  ENDIF
954  IF(right-left.LT.rqtol*abs(lambda)) THEN
955 * The eigenvalue is computed to bisection accuracy
956 * compute eigenvector and stop
957  usedbs = .true.
958  GOTO 120
959  ELSEIF( iter.LT.maxitr ) THEN
960  GOTO 120
961  ELSEIF( iter.EQ.maxitr ) THEN
962  needbs = .true.
963  GOTO 120
964  ELSE
965  info = 5
966  RETURN
967  END IF
968  ELSE
969  stp2ii = .false.
970  IF(usedrq .AND. usedbs .AND.
971  $ bstres.LE.resid) THEN
972  lambda = bstw
973  stp2ii = .true.
974  ENDIF
975  IF (stp2ii) THEN
976 * improve error angle by second step
977  CALL zlar1v( in, 1, in, lambda,
978  $ d( ibegin ), l( ibegin ),
979  $ work(indld+ibegin-1),
980  $ work(indlld+ibegin-1),
981  $ pivmin, gaptol, z( ibegin, windex ),
982  $ .NOT.usedbs, negcnt, ztz, mingma,
983  $ iwork( iindr+windex ),
984  $ isuppz( 2*windex-1 ),
985  $ nrminv, resid, rqcorr, work( indwrk ) )
986  ENDIF
987  work( windex ) = lambda
988  END IF
989 *
990 * Compute FP-vector support w.r.t. whole matrix
991 *
992  isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
993  isuppz( 2*windex ) = isuppz( 2*windex )+oldien
994  zfrom = isuppz( 2*windex-1 )
995  zto = isuppz( 2*windex )
996  isupmn = isupmn + oldien
997  isupmx = isupmx + oldien
998 * Ensure vector is ok if support in the RQI has changed
999  IF(isupmn.LT.zfrom) THEN
1000  DO 122 ii = isupmn,zfrom-1
1001  z( ii, windex ) = zero
1002  122 CONTINUE
1003  ENDIF
1004  IF(isupmx.GT.zto) THEN
1005  DO 123 ii = zto+1,isupmx
1006  z( ii, windex ) = zero
1007  123 CONTINUE
1008  ENDIF
1009  CALL zdscal( zto-zfrom+1, nrminv,
1010  $ z( zfrom, windex ), 1 )
1011  125 CONTINUE
1012 * Update W
1013  w( windex ) = lambda+sigma
1014 * Recompute the gaps on the left and right
1015 * But only allow them to become larger and not
1016 * smaller (which can only happen through "bad"
1017 * cancellation and doesn't reflect the theory
1018 * where the initial gaps are underestimated due
1019 * to WERR being too crude.)
1020  IF(.NOT.eskip) THEN
1021  IF( k.GT.1) THEN
1022  wgap( windmn ) = max( wgap(windmn),
1023  $ w(windex)-werr(windex)
1024  $ - w(windmn)-werr(windmn) )
1025  ENDIF
1026  IF( windex.LT.wend ) THEN
1027  wgap( windex ) = max( savgap,
1028  $ w( windpl )-werr( windpl )
1029  $ - w( windex )-werr( windex) )
1030  ENDIF
1031  ENDIF
1032  idone = idone + 1
1033  ENDIF
1034 * here ends the code for the current child
1035 *
1036  139 CONTINUE
1037 * Proceed to any remaining child nodes
1038  newfst = j + 1
1039  140 CONTINUE
1040  150 CONTINUE
1041  ndepth = ndepth + 1
1042  GO TO 40
1043  END IF
1044  ibegin = iend + 1
1045  wbegin = wend + 1
1046  170 CONTINUE
1047 *
1048 
1049  RETURN
1050 *
1051 * End of ZLARRV
1052 *
1053  END
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
zlar1v
subroutine zlar1v(N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)
ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the...
Definition: zlar1v.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
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
zlarrv
subroutine zlarrv(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)
ZLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition: zlarrv.f:288
zdscal
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54