LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
zlarrv.f
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 < 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 *> \ingroup larrv
270 *
271 *> \par Contributors:
272 * ==================
273 *>
274 *> Beresford Parlett, University of California, Berkeley, USA \n
275 *> Jim Demmel, University of California, Berkeley, USA \n
276 *> Inderjit Dhillon, University of Texas, Austin, USA \n
277 *> Osni Marques, LBNL/NERSC, USA \n
278 *> Christof Voemel, University of California, Berkeley, USA
279 *
280 * =====================================================================
281  SUBROUTINE zlarrv( N, VL, VU, D, L, PIVMIN,
282  $ ISPLIT, M, DOL, DOU, MINRGP,
283  $ RTOL1, RTOL2, W, WERR, WGAP,
284  $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
285  $ WORK, IWORK, INFO )
286 *
287 * -- LAPACK auxiliary routine --
288 * -- LAPACK is a software package provided by Univ. of Tennessee, --
289 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290 *
291 * .. Scalar Arguments ..
292  INTEGER DOL, DOU, INFO, LDZ, M, N
293  DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
294 * ..
295 * .. Array Arguments ..
296  INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
297  $ isuppz( * ), iwork( * )
298  DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
299  $ wgap( * ), work( * )
300  COMPLEX*16 Z( ldz, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  INTEGER MAXITR
307  parameter( maxitr = 10 )
308  COMPLEX*16 CZERO
309  parameter( czero = ( 0.0d0, 0.0d0 ) )
310  DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, HALF
311  parameter( zero = 0.0d0, one = 1.0d0,
312  $ two = 2.0d0, three = 3.0d0,
313  $ four = 4.0d0, half = 0.5d0)
314 * ..
315 * .. Local Scalars ..
316  LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
317  INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
318  $ iindc2, iindr, iindwk, iinfo, im, in, indeig,
319  $ indld, indlld, indwrk, isupmn, isupmx, iter,
320  $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
321  $ ndepth, negcnt, newcls, newfst, newftt, newlst,
322  $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
323  $ oldncl, p, parity, q, wbegin, wend, windex,
324  $ windmn, windpl, zfrom, zto, zusedl, zusedu,
325  $ zusedw
326  INTEGER INDIN1, INDIN2
327  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, dlarrb, dlarrf, zdscal, zlar1v,
338  $ zlaset
339 * ..
340 * .. Intrinsic Functions ..
341  INTRINSIC abs, dble, max, min
342  INTRINSIC dcmplx
343 * ..
344 * .. Executable Statements ..
345 * ..
346 
347  info = 0
348 *
349 * Quick return if possible
350 *
351  IF( (n.LE.0).OR.(m.LE.0) ) THEN
352  RETURN
353  END IF
354 *
355 * The first N entries of WORK are reserved for the eigenvalues
356  indld = n+1
357  indlld= 2*n+1
358  indin1 = 3*n + 1
359  indin2 = 4*n + 1
360  indwrk = 5*n + 1
361  minwsize = 12 * n
362 
363  DO 5 i= 1,minwsize
364  work( i ) = zero
365  5 CONTINUE
366 
367 * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
368 * factorization used to compute the FP vector
369  iindr = 0
370 * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
371 * layer and the one above.
372  iindc1 = n
373  iindc2 = 2*n
374  iindwk = 3*n + 1
375 
376  miniwsize = 7 * n
377  DO 10 i= 1,miniwsize
378  iwork( i ) = 0
379  10 CONTINUE
380 
381  zusedl = 1
382  IF(dol.GT.1) THEN
383 * Set lower bound for use of Z
384  zusedl = dol-1
385  ENDIF
386  zusedu = m
387  IF(dou.LT.m) THEN
388 * Set lower bound for use of Z
389  zusedu = dou+1
390  ENDIF
391 * The width of the part of Z that is used
392  zusedw = zusedu - zusedl + 1
393 
394 
395  CALL zlaset( 'Full', n, zusedw, czero, czero,
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 ) = dcmplx( one, zero )
465  isuppz( 2*wbegin-1 ) = ibegin
466  isuppz( 2*wbegin ) = ibegin
467  w( wbegin ) = w( wbegin ) + sigma
468  work( wbegin ) = w( wbegin )
469  ibegin = iend + 1
470  wbegin = wbegin + 1
471  GO TO 170
472  END IF
473 
474 * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
475 * Note that these can be approximations, in this case, the corresp.
476 * entries of WERR give the size of the uncertainty interval.
477 * The eigenvalue approximations will be refined when necessary as
478 * high relative accuracy is required for the computation of the
479 * corresponding eigenvectors.
480  CALL 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  DO 45 k = 1, in - 1
557  d( ibegin+k-1 ) = dble( z( ibegin+k-1,
558  $ j ) )
559  l( ibegin+k-1 ) = dble( z( ibegin+k-1,
560  $ j+1 ) )
561  45 CONTINUE
562  d( iend ) = dble( z( iend, j ) )
563  sigma = dble( z( iend, j+1 ) )
564 
565 * Set the corresponding entries in Z to zero
566  CALL zlaset( 'Full', in, 2, czero, czero,
567  $ z( ibegin, j), ldz )
568  END IF
569 
570 * Compute DL and DLL of current RRR
571  DO 50 j = ibegin, iend-1
572  tmp = d( j )*l( j )
573  work( indld-1+j ) = tmp
574  work( indlld-1+j ) = tmp*l( j )
575  50 CONTINUE
576 
577  IF( ndepth.GT.0 ) THEN
578 * P and Q are index of the first and last eigenvalue to compute
579 * within the current block
580  p = indexw( wbegin-1+oldfst )
581  q = indexw( wbegin-1+oldlst )
582 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
583 * through the Q-OFFSET elements of these arrays are to be used.
584 * OFFSET = P-OLDFST
585  offset = indexw( wbegin ) - 1
586 * perform limited bisection (if necessary) to get approximate
587 * eigenvalues to the precision needed.
588  CALL dlarrb( in, d( ibegin ),
589  $ work(indlld+ibegin-1),
590  $ p, q, rtol1, rtol2, offset,
591  $ work(wbegin),wgap(wbegin),werr(wbegin),
592  $ work( indwrk ), iwork( iindwk ),
593  $ pivmin, spdiam, in, iinfo )
594  IF( iinfo.NE.0 ) THEN
595  info = -1
596  RETURN
597  ENDIF
598 * We also recompute the extremal gaps. W holds all eigenvalues
599 * of the unshifted matrix and must be used for computation
600 * of WGAP, the entries of WORK might stem from RRRs with
601 * different shifts. The gaps from WBEGIN-1+OLDFST to
602 * WBEGIN-1+OLDLST are correctly computed in DLARRB.
603 * However, we only allow the gaps to become greater since
604 * this is what should happen when we decrease WERR
605  IF( oldfst.GT.1) THEN
606  wgap( wbegin+oldfst-2 ) =
607  $ max(wgap(wbegin+oldfst-2),
608  $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
609  $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
610  ENDIF
611  IF( wbegin + oldlst -1 .LT. wend ) THEN
612  wgap( wbegin+oldlst-1 ) =
613  $ max(wgap(wbegin+oldlst-1),
614  $ w(wbegin+oldlst)-werr(wbegin+oldlst)
615  $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
616  ENDIF
617 * Each time the eigenvalues in WORK get refined, we store
618 * the newly found approximation with all shifts applied in W
619  DO 53 j=oldfst,oldlst
620  w(wbegin+j-1) = work(wbegin+j-1)+sigma
621  53 CONTINUE
622  END IF
623 
624 * Process the current node.
625  newfst = oldfst
626  DO 140 j = oldfst, oldlst
627  IF( j.EQ.oldlst ) THEN
628 * we are at the right end of the cluster, this is also the
629 * boundary of the child cluster
630  newlst = j
631  ELSE IF ( wgap( wbegin + j -1).GE.
632  $ minrgp* abs( work(wbegin + j -1) ) ) THEN
633 * the right relative gap is big enough, the child cluster
634 * (NEWFST,..,NEWLST) is well separated from the following
635  newlst = j
636  ELSE
637 * inside a child cluster, the relative gap is not
638 * big enough.
639  GOTO 140
640  END IF
641 
642 * Compute size of child cluster found
643  newsiz = newlst - newfst + 1
644 
645 * NEWFTT is the place in Z where the new RRR or the computed
646 * eigenvector is to be stored
647  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
648 * Store representation at location of the leftmost evalue
649 * of the cluster
650  newftt = wbegin + newfst - 1
651  ELSE
652  IF(wbegin+newfst-1.LT.dol) THEN
653 * Store representation at the left end of Z array
654  newftt = dol - 1
655  ELSEIF(wbegin+newfst-1.GT.dou) THEN
656 * Store representation at the right end of Z array
657  newftt = dou
658  ELSE
659  newftt = wbegin + newfst - 1
660  ENDIF
661  ENDIF
662 
663  IF( newsiz.GT.1) THEN
664 *
665 * Current child is not a singleton but a cluster.
666 * Compute and store new representation of child.
667 *
668 *
669 * Compute left and right cluster gap.
670 *
671 * LGAP and RGAP are not computed from WORK because
672 * the eigenvalue approximations may stem from RRRs
673 * different shifts. However, W hold all eigenvalues
674 * of the unshifted matrix. Still, the entries in WGAP
675 * have to be computed from WORK since the entries
676 * in W might be of the same order so that gaps are not
677 * exhibited correctly for very close eigenvalues.
678  IF( newfst.EQ.1 ) THEN
679  lgap = max( zero,
680  $ w(wbegin)-werr(wbegin) - vl )
681  ELSE
682  lgap = wgap( wbegin+newfst-2 )
683  ENDIF
684  rgap = wgap( wbegin+newlst-1 )
685 *
686 * Compute left- and rightmost eigenvalue of child
687 * to high precision in order to shift as close
688 * as possible and obtain as large relative gaps
689 * as possible
690 *
691  DO 55 k =1,2
692  IF(k.EQ.1) THEN
693  p = indexw( wbegin-1+newfst )
694  ELSE
695  p = indexw( wbegin-1+newlst )
696  ENDIF
697  offset = indexw( wbegin ) - 1
698  CALL dlarrb( in, d(ibegin),
699  $ work( indlld+ibegin-1 ),p,p,
700  $ rqtol, rqtol, offset,
701  $ work(wbegin),wgap(wbegin),
702  $ werr(wbegin),work( indwrk ),
703  $ iwork( iindwk ), pivmin, spdiam,
704  $ in, iinfo )
705  55 CONTINUE
706 *
707  IF((wbegin+newlst-1.LT.dol).OR.
708  $ (wbegin+newfst-1.GT.dou)) THEN
709 * if the cluster contains no desired eigenvalues
710 * skip the computation of that branch of the rep. tree
711 *
712 * We could skip before the refinement of the extremal
713 * eigenvalues of the child, but then the representation
714 * tree could be different from the one when nothing is
715 * skipped. For this reason we skip at this place.
716  idone = idone + newlst - newfst + 1
717  GOTO 139
718  ENDIF
719 *
720 * Compute RRR of child cluster.
721 * Note that the new RRR is stored in Z
722 *
723 * DLARRF needs LWORK = 2*N
724  CALL dlarrf( in, d( ibegin ), l( ibegin ),
725  $ work(indld+ibegin-1),
726  $ newfst, newlst, work(wbegin),
727  $ wgap(wbegin), werr(wbegin),
728  $ spdiam, lgap, rgap, pivmin, tau,
729  $ work( indin1 ), work( indin2 ),
730  $ work( indwrk ), iinfo )
731 * In the complex case, DLARRF cannot write
732 * the new RRR directly into Z and needs an intermediate
733 * workspace
734  DO 56 k = 1, in-1
735  z( ibegin+k-1, newftt ) =
736  $ dcmplx( work( indin1+k-1 ), zero )
737  z( ibegin+k-1, newftt+1 ) =
738  $ dcmplx( work( indin2+k-1 ), zero )
739  56 CONTINUE
740  z( iend, newftt ) =
741  $ dcmplx( work( indin1+in-1 ), zero )
742  IF( iinfo.EQ.0 ) THEN
743 * a new RRR for the cluster was found by DLARRF
744 * update shift and store it
745  ssigma = sigma + tau
746  z( iend, newftt+1 ) = dcmplx( ssigma, zero )
747 * WORK() are the midpoints and WERR() the semi-width
748 * Note that the entries in W are unchanged.
749  DO 116 k = newfst, newlst
750  fudge =
751  $ three*eps*abs(work(wbegin+k-1))
752  work( wbegin + k - 1 ) =
753  $ work( wbegin + k - 1) - tau
754  fudge = fudge +
755  $ four*eps*abs(work(wbegin+k-1))
756 * Fudge errors
757  werr( wbegin + k - 1 ) =
758  $ werr( wbegin + k - 1 ) + fudge
759 * Gaps are not fudged. Provided that WERR is small
760 * when eigenvalues are close, a zero gap indicates
761 * that a new representation is needed for resolving
762 * the cluster. A fudge could lead to a wrong decision
763 * of judging eigenvalues 'separated' which in
764 * reality are not. This could have a negative impact
765 * on the orthogonality of the computed eigenvectors.
766  116 CONTINUE
767 
768  nclus = nclus + 1
769  k = newcls + 2*nclus
770  iwork( k-1 ) = newfst
771  iwork( k ) = newlst
772  ELSE
773  info = -2
774  RETURN
775  ENDIF
776  ELSE
777 *
778 * Compute eigenvector of singleton
779 *
780  iter = 0
781 *
782  tol = four * log(dble(in)) * eps
783 *
784  k = newfst
785  windex = wbegin + k - 1
786  windmn = max(windex - 1,1)
787  windpl = min(windex + 1,m)
788  lambda = work( windex )
789  done = done + 1
790 * Check if eigenvector computation is to be skipped
791  IF((windex.LT.dol).OR.
792  $ (windex.GT.dou)) THEN
793  eskip = .true.
794  GOTO 125
795  ELSE
796  eskip = .false.
797  ENDIF
798  left = work( windex ) - werr( windex )
799  right = work( windex ) + werr( windex )
800  indeig = indexw( windex )
801 * Note that since we compute the eigenpairs for a child,
802 * all eigenvalue approximations are w.r.t the same shift.
803 * In this case, the entries in WORK should be used for
804 * computing the gaps since they exhibit even very small
805 * differences in the eigenvalues, as opposed to the
806 * entries in W which might "look" the same.
807 
808  IF( k .EQ. 1) THEN
809 * In the case RANGE='I' and with not much initial
810 * accuracy in LAMBDA and VL, the formula
811 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
812 * can lead to an overestimation of the left gap and
813 * thus to inadequately early RQI 'convergence'.
814 * Prevent this by forcing a small left gap.
815  lgap = eps*max(abs(left),abs(right))
816  ELSE
817  lgap = wgap(windmn)
818  ENDIF
819  IF( k .EQ. im) THEN
820 * In the case RANGE='I' and with not much initial
821 * accuracy in LAMBDA and VU, the formula
822 * can lead to an overestimation of the right gap and
823 * thus to inadequately early RQI 'convergence'.
824 * Prevent this by forcing a small right gap.
825  rgap = eps*max(abs(left),abs(right))
826  ELSE
827  rgap = wgap(windex)
828  ENDIF
829  gap = min( lgap, rgap )
830  IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
831 * The eigenvector support can become wrong
832 * because significant entries could be cut off due to a
833 * large GAPTOL parameter in LAR1V. Prevent this.
834  gaptol = zero
835  ELSE
836  gaptol = gap * eps
837  ENDIF
838  isupmn = in
839  isupmx = 1
840 * Update WGAP so that it holds the minimum gap
841 * to the left or the right. This is crucial in the
842 * case where bisection is used to ensure that the
843 * eigenvalue is refined up to the required precision.
844 * The correct value is restored afterwards.
845  savgap = wgap(windex)
846  wgap(windex) = gap
847 * We want to use the Rayleigh Quotient Correction
848 * as often as possible since it converges quadratically
849 * when we are close enough to the desired eigenvalue.
850 * However, the Rayleigh Quotient can have the wrong sign
851 * and lead us away from the desired eigenvalue. In this
852 * case, the best we can do is to use bisection.
853  usedbs = .false.
854  usedrq = .false.
855 * Bisection is initially turned off unless it is forced
856  needbs = .NOT.tryrqc
857  120 CONTINUE
858 * Check if bisection should be used to refine eigenvalue
859  IF(needbs) THEN
860 * Take the bisection as new iterate
861  usedbs = .true.
862  itmp1 = iwork( iindr+windex )
863  offset = indexw( wbegin ) - 1
864  CALL dlarrb( in, d(ibegin),
865  $ work(indlld+ibegin-1),indeig,indeig,
866  $ zero, two*eps, offset,
867  $ work(wbegin),wgap(wbegin),
868  $ werr(wbegin),work( indwrk ),
869  $ iwork( iindwk ), pivmin, spdiam,
870  $ itmp1, iinfo )
871  IF( iinfo.NE.0 ) THEN
872  info = -3
873  RETURN
874  ENDIF
875  lambda = work( windex )
876 * Reset twist index from inaccurate LAMBDA to
877 * force computation of true MINGMA
878  iwork( iindr+windex ) = 0
879  ENDIF
880 * Given LAMBDA, compute the eigenvector.
881  CALL zlar1v( in, 1, in, lambda, d( ibegin ),
882  $ l( ibegin ), work(indld+ibegin-1),
883  $ work(indlld+ibegin-1),
884  $ pivmin, gaptol, z( ibegin, windex ),
885  $ .NOT.usedbs, negcnt, ztz, mingma,
886  $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
887  $ nrminv, resid, rqcorr, work( indwrk ) )
888  IF(iter .EQ. 0) THEN
889  bstres = resid
890  bstw = lambda
891  ELSEIF(resid.LT.bstres) THEN
892  bstres = resid
893  bstw = lambda
894  ENDIF
895  isupmn = min(isupmn,isuppz( 2*windex-1 ))
896  isupmx = max(isupmx,isuppz( 2*windex ))
897  iter = iter + 1
898 
899 * sin alpha <= |resid|/gap
900 * Note that both the residual and the gap are
901 * proportional to the matrix, so ||T|| doesn't play
902 * a role in the quotient
903 
904 *
905 * Convergence test for Rayleigh-Quotient iteration
906 * (omitted when Bisection has been used)
907 *
908  IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
909  $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
910  $ THEN
911 * We need to check that the RQCORR update doesn't
912 * move the eigenvalue away from the desired one and
913 * towards a neighbor. -> protection with bisection
914  IF(indeig.LE.negcnt) THEN
915 * The wanted eigenvalue lies to the left
916  sgndef = -one
917  ELSE
918 * The wanted eigenvalue lies to the right
919  sgndef = one
920  ENDIF
921 * We only use the RQCORR if it improves the
922 * the iterate reasonably.
923  IF( ( rqcorr*sgndef.GE.zero )
924  $ .AND.( lambda + rqcorr.LE. right)
925  $ .AND.( lambda + rqcorr.GE. left)
926  $ ) THEN
927  usedrq = .true.
928 * Store new midpoint of bisection interval in WORK
929  IF(sgndef.EQ.one) THEN
930 * The current LAMBDA is on the left of the true
931 * eigenvalue
932  left = lambda
933 * We prefer to assume that the error estimate
934 * is correct. We could make the interval not
935 * as a bracket but to be modified if the RQCORR
936 * chooses to. In this case, the RIGHT side should
937 * be modified as follows:
938 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
939  ELSE
940 * The current LAMBDA is on the right of the true
941 * eigenvalue
942  right = lambda
943 * See comment about assuming the error estimate is
944 * correct above.
945 * LEFT = MIN(LEFT, LAMBDA + RQCORR)
946  ENDIF
947  work( windex ) =
948  $ half * (right + left)
949 * Take RQCORR since it has the correct sign and
950 * improves the iterate reasonably
951  lambda = lambda + rqcorr
952 * Update width of error interval
953  werr( windex ) =
954  $ half * (right-left)
955  ELSE
956  needbs = .true.
957  ENDIF
958  IF(right-left.LT.rqtol*abs(lambda)) THEN
959 * The eigenvalue is computed to bisection accuracy
960 * compute eigenvector and stop
961  usedbs = .true.
962  GOTO 120
963  ELSEIF( iter.LT.maxitr ) THEN
964  GOTO 120
965  ELSEIF( iter.EQ.maxitr ) THEN
966  needbs = .true.
967  GOTO 120
968  ELSE
969  info = 5
970  RETURN
971  END IF
972  ELSE
973  stp2ii = .false.
974  IF(usedrq .AND. usedbs .AND.
975  $ bstres.LE.resid) THEN
976  lambda = bstw
977  stp2ii = .true.
978  ENDIF
979  IF (stp2ii) THEN
980 * improve error angle by second step
981  CALL zlar1v( in, 1, in, lambda,
982  $ d( ibegin ), l( ibegin ),
983  $ work(indld+ibegin-1),
984  $ work(indlld+ibegin-1),
985  $ pivmin, gaptol, z( ibegin, windex ),
986  $ .NOT.usedbs, negcnt, ztz, mingma,
987  $ iwork( iindr+windex ),
988  $ isuppz( 2*windex-1 ),
989  $ nrminv, resid, rqcorr, work( indwrk ) )
990  ENDIF
991  work( windex ) = lambda
992  END IF
993 *
994 * Compute FP-vector support w.r.t. whole matrix
995 *
996  isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
997  isuppz( 2*windex ) = isuppz( 2*windex )+oldien
998  zfrom = isuppz( 2*windex-1 )
999  zto = isuppz( 2*windex )
1000  isupmn = isupmn + oldien
1001  isupmx = isupmx + oldien
1002 * Ensure vector is ok if support in the RQI has changed
1003  IF(isupmn.LT.zfrom) THEN
1004  DO 122 ii = isupmn,zfrom-1
1005  z( ii, windex ) = zero
1006  122 CONTINUE
1007  ENDIF
1008  IF(isupmx.GT.zto) THEN
1009  DO 123 ii = zto+1,isupmx
1010  z( ii, windex ) = zero
1011  123 CONTINUE
1012  ENDIF
1013  CALL zdscal( zto-zfrom+1, nrminv,
1014  $ z( zfrom, windex ), 1 )
1015  125 CONTINUE
1016 * Update W
1017  w( windex ) = lambda+sigma
1018 * Recompute the gaps on the left and right
1019 * But only allow them to become larger and not
1020 * smaller (which can only happen through "bad"
1021 * cancellation and doesn't reflect the theory
1022 * where the initial gaps are underestimated due
1023 * to WERR being too crude.)
1024  IF(.NOT.eskip) THEN
1025  IF( k.GT.1) THEN
1026  wgap( windmn ) = max( wgap(windmn),
1027  $ w(windex)-werr(windex)
1028  $ - w(windmn)-werr(windmn) )
1029  ENDIF
1030  IF( windex.LT.wend ) THEN
1031  wgap( windex ) = max( savgap,
1032  $ w( windpl )-werr( windpl )
1033  $ - w( windex )-werr( windex) )
1034  ENDIF
1035  ENDIF
1036  idone = idone + 1
1037  ENDIF
1038 * here ends the code for the current child
1039 *
1040  139 CONTINUE
1041 * Proceed to any remaining child nodes
1042  newfst = j + 1
1043  140 CONTINUE
1044  150 CONTINUE
1045  ndepth = ndepth + 1
1046  GO TO 40
1047  END IF
1048  ibegin = iend + 1
1049  wbegin = wend + 1
1050  170 CONTINUE
1051 *
1052 
1053  RETURN
1054 *
1055 * End of ZLARRV
1056 *
1057  END
zdscal
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
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
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:230
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:286
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:106
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
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82