LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 *> \endverbatim
63 *>
64 *> \param[in] VU
65 *> \verbatim
66 *> VU is DOUBLE PRECISION
67 *> Lower and upper bounds of the interval that contains the desired
68 *> eigenvalues. VL < VU. Needed to compute gaps on the left or right
69 *> end of the extremal eigenvalues in the desired RANGE.
70 *> \endverbatim
71 *>
72 *> \param[in,out] D
73 *> \verbatim
74 *> D is DOUBLE PRECISION array, dimension (N)
75 *> On entry, the N diagonal elements of the diagonal matrix D.
76 *> On exit, D may be overwritten.
77 *> \endverbatim
78 *>
79 *> \param[in,out] L
80 *> \verbatim
81 *> L is DOUBLE PRECISION array, dimension (N)
82 *> On entry, the (N-1) subdiagonal elements of the unit
83 *> bidiagonal matrix L are in elements 1 to N-1 of L
84 *> (if the matrix is not splitted.) At the end of each block
85 *> is stored the corresponding shift as given by DLARRE.
86 *> On exit, L is overwritten.
87 *> \endverbatim
88 *>
89 *> \param[in] PIVMIN
90 *> \verbatim
91 *> PIVMIN is DOUBLE PRECISION
92 *> The minimum pivot allowed in the Sturm sequence.
93 *> \endverbatim
94 *>
95 *> \param[in] ISPLIT
96 *> \verbatim
97 *> ISPLIT is INTEGER array, dimension (N)
98 *> The splitting points, at which T breaks up into blocks.
99 *> The first block consists of rows/columns 1 to
100 *> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
101 *> through ISPLIT( 2 ), etc.
102 *> \endverbatim
103 *>
104 *> \param[in] M
105 *> \verbatim
106 *> M is INTEGER
107 *> The total number of input eigenvalues. 0 <= M <= N.
108 *> \endverbatim
109 *>
110 *> \param[in] DOL
111 *> \verbatim
112 *> DOL is INTEGER
113 *> \endverbatim
114 *>
115 *> \param[in] DOU
116 *> \verbatim
117 *> DOU is INTEGER
118 *> If the user wants to compute only selected eigenvectors from all
119 *> the eigenvalues supplied, he can specify an index range DOL:DOU.
120 *> Or else the setting DOL=1, DOU=M should be applied.
121 *> Note that DOL and DOU refer to the order in which the eigenvalues
122 *> are stored in W.
123 *> If the user wants to compute only selected eigenpairs, then
124 *> the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
125 *> computed eigenvectors. All other columns of Z are set to zero.
126 *> \endverbatim
127 *>
128 *> \param[in] MINRGP
129 *> \verbatim
130 *> MINRGP is DOUBLE PRECISION
131 *> \endverbatim
132 *>
133 *> \param[in] RTOL1
134 *> \verbatim
135 *> RTOL1 is DOUBLE PRECISION
136 *> \endverbatim
137 *>
138 *> \param[in] RTOL2
139 *> \verbatim
140 *> RTOL2 is DOUBLE PRECISION
141 *> Parameters for bisection.
142 *> An interval [LEFT,RIGHT] has converged if
143 *> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
144 *> \endverbatim
145 *>
146 *> \param[in,out] W
147 *> \verbatim
148 *> W is DOUBLE PRECISION array, dimension (N)
149 *> The first M elements of W contain the APPROXIMATE eigenvalues for
150 *> which eigenvectors are to be computed. The eigenvalues
151 *> should be grouped by split-off block and ordered from
152 *> smallest to largest within the block ( The output array
153 *> W from DLARRE is expected here ). Furthermore, they are with
154 *> respect to the shift of the corresponding root representation
155 *> for their block. On exit, W holds the eigenvalues of the
156 *> UNshifted matrix.
157 *> \endverbatim
158 *>
159 *> \param[in,out] WERR
160 *> \verbatim
161 *> WERR is DOUBLE PRECISION array, dimension (N)
162 *> The first M elements contain the semiwidth of the uncertainty
163 *> interval of the corresponding eigenvalue in W
164 *> \endverbatim
165 *>
166 *> \param[in,out] WGAP
167 *> \verbatim
168 *> WGAP is DOUBLE PRECISION array, dimension (N)
169 *> The separation from the right neighbor eigenvalue in W.
170 *> \endverbatim
171 *>
172 *> \param[in] IBLOCK
173 *> \verbatim
174 *> IBLOCK is INTEGER array, dimension (N)
175 *> The indices of the blocks (submatrices) associated with the
176 *> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
177 *> W(i) belongs to the first block from the top, =2 if W(i)
178 *> belongs to the second block, etc.
179 *> \endverbatim
180 *>
181 *> \param[in] INDEXW
182 *> \verbatim
183 *> INDEXW is INTEGER array, dimension (N)
184 *> The indices of the eigenvalues within each block (submatrix);
185 *> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
186 *> i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
187 *> \endverbatim
188 *>
189 *> \param[in] GERS
190 *> \verbatim
191 *> GERS is DOUBLE PRECISION array, dimension (2*N)
192 *> The N Gerschgorin intervals (the i-th Gerschgorin interval
193 *> is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
194 *> be computed from the original UNshifted matrix.
195 *> \endverbatim
196 *>
197 *> \param[out] Z
198 *> \verbatim
199 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M) )
200 *> If INFO = 0, the first M columns of Z contain the
201 *> orthonormal eigenvectors of the matrix T
202 *> corresponding to the input eigenvalues, with the i-th
203 *> column of Z holding the eigenvector associated with W(i).
204 *> Note: the user must ensure that at least max(1,M) columns are
205 *> supplied in the array Z.
206 *> \endverbatim
207 *>
208 *> \param[in] LDZ
209 *> \verbatim
210 *> LDZ is INTEGER
211 *> The leading dimension of the array Z. LDZ >= 1, and if
212 *> JOBZ = 'V', LDZ >= max(1,N).
213 *> \endverbatim
214 *>
215 *> \param[out] ISUPPZ
216 *> \verbatim
217 *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
218 *> The support of the eigenvectors in Z, i.e., the indices
219 *> indicating the nonzero elements in Z. The I-th eigenvector
220 *> is nonzero only in elements ISUPPZ( 2*I-1 ) through
221 *> ISUPPZ( 2*I ).
222 *> \endverbatim
223 *>
224 *> \param[out] WORK
225 *> \verbatim
226 *> WORK is DOUBLE PRECISION array, dimension (12*N)
227 *> \endverbatim
228 *>
229 *> \param[out] IWORK
230 *> \verbatim
231 *> IWORK is INTEGER array, dimension (7*N)
232 *> \endverbatim
233 *>
234 *> \param[out] INFO
235 *> \verbatim
236 *> INFO is INTEGER
237 *> = 0: successful exit
238 *>
239 *> > 0: A problem occured in ZLARRV.
240 *> < 0: One of the called subroutines signaled an internal problem.
241 *> Needs inspection of the corresponding parameter IINFO
242 *> for further information.
243 *>
244 *> =-1: Problem in DLARRB when refining a child's eigenvalues.
245 *> =-2: Problem in DLARRF when computing the RRR of a child.
246 *> When a child is inside a tight cluster, it can be difficult
247 *> to find an RRR. A partial remedy from the user's point of
248 *> view is to make the parameter MINRGP smaller and recompile.
249 *> However, as the orthogonality of the computed vectors is
250 *> proportional to 1/MINRGP, the user should be aware that
251 *> he might be trading in precision when he decreases MINRGP.
252 *> =-3: Problem in DLARRB when refining a single eigenvalue
253 *> after the Rayleigh correction was rejected.
254 *> = 5: The Rayleigh Quotient Iteration failed to converge to
255 *> full accuracy in MAXITR steps.
256 *> \endverbatim
257 *
258 * Authors:
259 * ========
260 *
261 *> \author Univ. of Tennessee
262 *> \author Univ. of California Berkeley
263 *> \author Univ. of Colorado Denver
264 *> \author NAG Ltd.
265 *
266 *> \date September 2012
267 *
268 *> \ingroup complex16OTHERauxiliary
269 *
270 *> \par Contributors:
271 * ==================
272 *>
273 *> Beresford Parlett, University of California, Berkeley, USA \n
274 *> Jim Demmel, University of California, Berkeley, USA \n
275 *> Inderjit Dhillon, University of Texas, Austin, USA \n
276 *> Osni Marques, LBNL/NERSC, USA \n
277 *> Christof Voemel, University of California, Berkeley, USA
278 *
279 * =====================================================================
280  SUBROUTINE zlarrv( N, VL, VU, D, L, PIVMIN,
281  $ isplit, m, dol, dou, minrgp,
282  $ rtol1, rtol2, w, werr, wgap,
283  $ iblock, indexw, gers, z, ldz, isuppz,
284  $ work, iwork, info )
285 *
286 * -- LAPACK auxiliary routine (version 3.4.2) --
287 * -- LAPACK is a software package provided by Univ. of Tennessee, --
288 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
289 * September 2012
290 *
291 * .. Scalar Arguments ..
292  INTEGER dol, dou, info, ldz, m, n
293  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 * The first N entries of WORK are reserved for the eigenvalues
348  indld = n+1
349  indlld= 2*n+1
350  indin1 = 3*n + 1
351  indin2 = 4*n + 1
352  indwrk = 5*n + 1
353  minwsize = 12 * n
354 
355  DO 5 i= 1,minwsize
356  work( i ) = zero
357  5 CONTINUE
358 
359 * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
360 * factorization used to compute the FP vector
361  iindr = 0
362 * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
363 * layer and the one above.
364  iindc1 = n
365  iindc2 = 2*n
366  iindwk = 3*n + 1
367 
368  miniwsize = 7 * n
369  DO 10 i= 1,miniwsize
370  iwork( i ) = 0
371  10 CONTINUE
372 
373  zusedl = 1
374  IF(dol.GT.1) THEN
375 * Set lower bound for use of Z
376  zusedl = dol-1
377  ENDIF
378  zusedu = m
379  IF(dou.LT.m) THEN
380 * Set lower bound for use of Z
381  zusedu = dou+1
382  ENDIF
383 * The width of the part of Z that is used
384  zusedw = zusedu - zusedl + 1
385 
386 
387  CALL zlaset( 'Full', n, zusedw, czero, czero,
388  $ z(1,zusedl), ldz )
389 
390  eps = dlamch( 'Precision' )
391  rqtol = two * eps
392 *
393 * Set expert flags for standard code.
394  tryrqc = .true.
395 
396  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
397  ELSE
398 * Only selected eigenpairs are computed. Since the other evalues
399 * are not refined by RQ iteration, bisection has to compute to full
400 * accuracy.
401  rtol1 = four * eps
402  rtol2 = four * eps
403  ENDIF
404 
405 * The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
406 * desired eigenvalues. The support of the nonzero eigenvector
407 * entries is contained in the interval IBEGIN:IEND.
408 * Remark that if k eigenpairs are desired, then the eigenvectors
409 * are stored in k contiguous columns of Z.
410 
411 * DONE is the number of eigenvectors already computed
412  done = 0
413  ibegin = 1
414  wbegin = 1
415  DO 170 jblk = 1, iblock( m )
416  iend = isplit( jblk )
417  sigma = l( iend )
418 * Find the eigenvectors of the submatrix indexed IBEGIN
419 * through IEND.
420  wend = wbegin - 1
421  15 CONTINUE
422  IF( wend.LT.m ) THEN
423  IF( iblock( wend+1 ).EQ.jblk ) THEN
424  wend = wend + 1
425  go to 15
426  END IF
427  END IF
428  IF( wend.LT.wbegin ) THEN
429  ibegin = iend + 1
430  go to 170
431  ELSEIF( (wend.LT.dol).OR.(wbegin.GT.dou) ) THEN
432  ibegin = iend + 1
433  wbegin = wend + 1
434  go to 170
435  END IF
436 
437 * Find local spectral diameter of the block
438  gl = gers( 2*ibegin-1 )
439  gu = gers( 2*ibegin )
440  DO 20 i = ibegin+1 , iend
441  gl = min( gers( 2*i-1 ), gl )
442  gu = max( gers( 2*i ), gu )
443  20 CONTINUE
444  spdiam = gu - gl
445 
446 * OLDIEN is the last index of the previous block
447  oldien = ibegin - 1
448 * Calculate the size of the current block
449  in = iend - ibegin + 1
450 * The number of eigenvalues in the current block
451  im = wend - wbegin + 1
452 
453 * This is for a 1x1 block
454  IF( ibegin.EQ.iend ) THEN
455  done = done+1
456  z( ibegin, wbegin ) = dcmplx( one, zero )
457  isuppz( 2*wbegin-1 ) = ibegin
458  isuppz( 2*wbegin ) = ibegin
459  w( wbegin ) = w( wbegin ) + sigma
460  work( wbegin ) = w( wbegin )
461  ibegin = iend + 1
462  wbegin = wbegin + 1
463  go to 170
464  END IF
465 
466 * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
467 * Note that these can be approximations, in this case, the corresp.
468 * entries of WERR give the size of the uncertainty interval.
469 * The eigenvalue approximations will be refined when necessary as
470 * high relative accuracy is required for the computation of the
471 * corresponding eigenvectors.
472  CALL dcopy( im, w( wbegin ), 1,
473  $ work( wbegin ), 1 )
474 
475 * We store in W the eigenvalue approximations w.r.t. the original
476 * matrix T.
477  DO 30 i=1,im
478  w(wbegin+i-1) = w(wbegin+i-1)+sigma
479  30 CONTINUE
480 
481 
482 * NDEPTH is the current depth of the representation tree
483  ndepth = 0
484 * PARITY is either 1 or 0
485  parity = 1
486 * NCLUS is the number of clusters for the next level of the
487 * representation tree, we start with NCLUS = 1 for the root
488  nclus = 1
489  iwork( iindc1+1 ) = 1
490  iwork( iindc1+2 ) = im
491 
492 * IDONE is the number of eigenvectors already computed in the current
493 * block
494  idone = 0
495 * loop while( IDONE.LT.IM )
496 * generate the representation tree for the current block and
497 * compute the eigenvectors
498  40 CONTINUE
499  IF( idone.LT.im ) THEN
500 * This is a crude protection against infinitely deep trees
501  IF( ndepth.GT.m ) THEN
502  info = -2
503  RETURN
504  ENDIF
505 * breadth first processing of the current level of the representation
506 * tree: OLDNCL = number of clusters on current level
507  oldncl = nclus
508 * reset NCLUS to count the number of child clusters
509  nclus = 0
510 *
511  parity = 1 - parity
512  IF( parity.EQ.0 ) THEN
513  oldcls = iindc1
514  newcls = iindc2
515  ELSE
516  oldcls = iindc2
517  newcls = iindc1
518  END IF
519 * Process the clusters on the current level
520  DO 150 i = 1, oldncl
521  j = oldcls + 2*i
522 * OLDFST, OLDLST = first, last index of current cluster.
523 * cluster indices start with 1 and are relative
524 * to WBEGIN when accessing W, WGAP, WERR, Z
525  oldfst = iwork( j-1 )
526  oldlst = iwork( j )
527  IF( ndepth.GT.0 ) THEN
528 * Retrieve relatively robust representation (RRR) of cluster
529 * that has been computed at the previous level
530 * The RRR is stored in Z and overwritten once the eigenvectors
531 * have been computed or when the cluster is refined
532 
533  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
534 * Get representation from location of the leftmost evalue
535 * of the cluster
536  j = wbegin + oldfst - 1
537  ELSE
538  IF(wbegin+oldfst-1.LT.dol) THEN
539 * Get representation from the left end of Z array
540  j = dol - 1
541  ELSEIF(wbegin+oldfst-1.GT.dou) THEN
542 * Get representation from the right end of Z array
543  j = dou
544  ELSE
545  j = wbegin + oldfst - 1
546  ENDIF
547  ENDIF
548  DO 45 k = 1, in - 1
549  d( ibegin+k-1 ) = dble( z( ibegin+k-1,
550  $ j ) )
551  l( ibegin+k-1 ) = dble( z( ibegin+k-1,
552  $ j+1 ) )
553  45 CONTINUE
554  d( iend ) = dble( z( iend, j ) )
555  sigma = dble( z( iend, j+1 ) )
556 
557 * Set the corresponding entries in Z to zero
558  CALL zlaset( 'Full', in, 2, czero, czero,
559  $ z( ibegin, j), ldz )
560  END IF
561 
562 * Compute DL and DLL of current RRR
563  DO 50 j = ibegin, iend-1
564  tmp = d( j )*l( j )
565  work( indld-1+j ) = tmp
566  work( indlld-1+j ) = tmp*l( j )
567  50 CONTINUE
568 
569  IF( ndepth.GT.0 ) THEN
570 * P and Q are index of the first and last eigenvalue to compute
571 * within the current block
572  p = indexw( wbegin-1+oldfst )
573  q = indexw( wbegin-1+oldlst )
574 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
575 * through the Q-OFFSET elements of these arrays are to be used.
576 * OFFSET = P-OLDFST
577  offset = indexw( wbegin ) - 1
578 * perform limited bisection (if necessary) to get approximate
579 * eigenvalues to the precision needed.
580  CALL dlarrb( in, d( ibegin ),
581  $ work(indlld+ibegin-1),
582  $ p, q, rtol1, rtol2, offset,
583  $ work(wbegin),wgap(wbegin),werr(wbegin),
584  $ work( indwrk ), iwork( iindwk ),
585  $ pivmin, spdiam, in, iinfo )
586  IF( iinfo.NE.0 ) THEN
587  info = -1
588  RETURN
589  ENDIF
590 * We also recompute the extremal gaps. W holds all eigenvalues
591 * of the unshifted matrix and must be used for computation
592 * of WGAP, the entries of WORK might stem from RRRs with
593 * different shifts. The gaps from WBEGIN-1+OLDFST to
594 * WBEGIN-1+OLDLST are correctly computed in DLARRB.
595 * However, we only allow the gaps to become greater since
596 * this is what should happen when we decrease WERR
597  IF( oldfst.GT.1) THEN
598  wgap( wbegin+oldfst-2 ) =
599  $ max(wgap(wbegin+oldfst-2),
600  $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
601  $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
602  ENDIF
603  IF( wbegin + oldlst -1 .LT. wend ) THEN
604  wgap( wbegin+oldlst-1 ) =
605  $ max(wgap(wbegin+oldlst-1),
606  $ w(wbegin+oldlst)-werr(wbegin+oldlst)
607  $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
608  ENDIF
609 * Each time the eigenvalues in WORK get refined, we store
610 * the newly found approximation with all shifts applied in W
611  DO 53 j=oldfst,oldlst
612  w(wbegin+j-1) = work(wbegin+j-1)+sigma
613  53 CONTINUE
614  END IF
615 
616 * Process the current node.
617  newfst = oldfst
618  DO 140 j = oldfst, oldlst
619  IF( j.EQ.oldlst ) THEN
620 * we are at the right end of the cluster, this is also the
621 * boundary of the child cluster
622  newlst = j
623  ELSE IF ( wgap( wbegin + j -1).GE.
624  $ minrgp* abs( work(wbegin + j -1) ) ) THEN
625 * the right relative gap is big enough, the child cluster
626 * (NEWFST,..,NEWLST) is well separated from the following
627  newlst = j
628  ELSE
629 * inside a child cluster, the relative gap is not
630 * big enough.
631  goto 140
632  END IF
633 
634 * Compute size of child cluster found
635  newsiz = newlst - newfst + 1
636 
637 * NEWFTT is the place in Z where the new RRR or the computed
638 * eigenvector is to be stored
639  IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
640 * Store representation at location of the leftmost evalue
641 * of the cluster
642  newftt = wbegin + newfst - 1
643  ELSE
644  IF(wbegin+newfst-1.LT.dol) THEN
645 * Store representation at the left end of Z array
646  newftt = dol - 1
647  ELSEIF(wbegin+newfst-1.GT.dou) THEN
648 * Store representation at the right end of Z array
649  newftt = dou
650  ELSE
651  newftt = wbegin + newfst - 1
652  ENDIF
653  ENDIF
654 
655  IF( newsiz.GT.1) THEN
656 *
657 * Current child is not a singleton but a cluster.
658 * Compute and store new representation of child.
659 *
660 *
661 * Compute left and right cluster gap.
662 *
663 * LGAP and RGAP are not computed from WORK because
664 * the eigenvalue approximations may stem from RRRs
665 * different shifts. However, W hold all eigenvalues
666 * of the unshifted matrix. Still, the entries in WGAP
667 * have to be computed from WORK since the entries
668 * in W might be of the same order so that gaps are not
669 * exhibited correctly for very close eigenvalues.
670  IF( newfst.EQ.1 ) THEN
671  lgap = max( zero,
672  $ w(wbegin)-werr(wbegin) - vl )
673  ELSE
674  lgap = wgap( wbegin+newfst-2 )
675  ENDIF
676  rgap = wgap( wbegin+newlst-1 )
677 *
678 * Compute left- and rightmost eigenvalue of child
679 * to high precision in order to shift as close
680 * as possible and obtain as large relative gaps
681 * as possible
682 *
683  DO 55 k =1,2
684  IF(k.EQ.1) THEN
685  p = indexw( wbegin-1+newfst )
686  ELSE
687  p = indexw( wbegin-1+newlst )
688  ENDIF
689  offset = indexw( wbegin ) - 1
690  CALL dlarrb( in, d(ibegin),
691  $ work( indlld+ibegin-1 ),p,p,
692  $ rqtol, rqtol, offset,
693  $ work(wbegin),wgap(wbegin),
694  $ werr(wbegin),work( indwrk ),
695  $ iwork( iindwk ), pivmin, spdiam,
696  $ in, iinfo )
697  55 CONTINUE
698 *
699  IF((wbegin+newlst-1.LT.dol).OR.
700  $ (wbegin+newfst-1.GT.dou)) THEN
701 * if the cluster contains no desired eigenvalues
702 * skip the computation of that branch of the rep. tree
703 *
704 * We could skip before the refinement of the extremal
705 * eigenvalues of the child, but then the representation
706 * tree could be different from the one when nothing is
707 * skipped. For this reason we skip at this place.
708  idone = idone + newlst - newfst + 1
709  goto 139
710  ENDIF
711 *
712 * Compute RRR of child cluster.
713 * Note that the new RRR is stored in Z
714 *
715 * DLARRF needs LWORK = 2*N
716  CALL dlarrf( in, d( ibegin ), l( ibegin ),
717  $ work(indld+ibegin-1),
718  $ newfst, newlst, work(wbegin),
719  $ wgap(wbegin), werr(wbegin),
720  $ spdiam, lgap, rgap, pivmin, tau,
721  $ work( indin1 ), work( indin2 ),
722  $ work( indwrk ), iinfo )
723 * In the complex case, DLARRF cannot write
724 * the new RRR directly into Z and needs an intermediate
725 * workspace
726  DO 56 k = 1, in-1
727  z( ibegin+k-1, newftt ) =
728  $ dcmplx( work( indin1+k-1 ), zero )
729  z( ibegin+k-1, newftt+1 ) =
730  $ dcmplx( work( indin2+k-1 ), zero )
731  56 CONTINUE
732  z( iend, newftt ) =
733  $ dcmplx( work( indin1+in-1 ), zero )
734  IF( iinfo.EQ.0 ) THEN
735 * a new RRR for the cluster was found by DLARRF
736 * update shift and store it
737  ssigma = sigma + tau
738  z( iend, newftt+1 ) = dcmplx( ssigma, zero )
739 * WORK() are the midpoints and WERR() the semi-width
740 * Note that the entries in W are unchanged.
741  DO 116 k = newfst, newlst
742  fudge =
743  $ three*eps*abs(work(wbegin+k-1))
744  work( wbegin + k - 1 ) =
745  $ work( wbegin + k - 1) - tau
746  fudge = fudge +
747  $ four*eps*abs(work(wbegin+k-1))
748 * Fudge errors
749  werr( wbegin + k - 1 ) =
750  $ werr( wbegin + k - 1 ) + fudge
751 * Gaps are not fudged. Provided that WERR is small
752 * when eigenvalues are close, a zero gap indicates
753 * that a new representation is needed for resolving
754 * the cluster. A fudge could lead to a wrong decision
755 * of judging eigenvalues 'separated' which in
756 * reality are not. This could have a negative impact
757 * on the orthogonality of the computed eigenvectors.
758  116 CONTINUE
759 
760  nclus = nclus + 1
761  k = newcls + 2*nclus
762  iwork( k-1 ) = newfst
763  iwork( k ) = newlst
764  ELSE
765  info = -2
766  RETURN
767  ENDIF
768  ELSE
769 *
770 * Compute eigenvector of singleton
771 *
772  iter = 0
773 *
774  tol = four * log(dble(in)) * eps
775 *
776  k = newfst
777  windex = wbegin + k - 1
778  windmn = max(windex - 1,1)
779  windpl = min(windex + 1,m)
780  lambda = work( windex )
781  done = done + 1
782 * Check if eigenvector computation is to be skipped
783  IF((windex.LT.dol).OR.
784  $ (windex.GT.dou)) THEN
785  eskip = .true.
786  goto 125
787  ELSE
788  eskip = .false.
789  ENDIF
790  left = work( windex ) - werr( windex )
791  right = work( windex ) + werr( windex )
792  indeig = indexw( windex )
793 * Note that since we compute the eigenpairs for a child,
794 * all eigenvalue approximations are w.r.t the same shift.
795 * In this case, the entries in WORK should be used for
796 * computing the gaps since they exhibit even very small
797 * differences in the eigenvalues, as opposed to the
798 * entries in W which might "look" the same.
799 
800  IF( k .EQ. 1) THEN
801 * In the case RANGE='I' and with not much initial
802 * accuracy in LAMBDA and VL, the formula
803 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
804 * can lead to an overestimation of the left gap and
805 * thus to inadequately early RQI 'convergence'.
806 * Prevent this by forcing a small left gap.
807  lgap = eps*max(abs(left),abs(right))
808  ELSE
809  lgap = wgap(windmn)
810  ENDIF
811  IF( k .EQ. im) THEN
812 * In the case RANGE='I' and with not much initial
813 * accuracy in LAMBDA and VU, the formula
814 * can lead to an overestimation of the right gap and
815 * thus to inadequately early RQI 'convergence'.
816 * Prevent this by forcing a small right gap.
817  rgap = eps*max(abs(left),abs(right))
818  ELSE
819  rgap = wgap(windex)
820  ENDIF
821  gap = min( lgap, rgap )
822  IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
823 * The eigenvector support can become wrong
824 * because significant entries could be cut off due to a
825 * large GAPTOL parameter in LAR1V. Prevent this.
826  gaptol = zero
827  ELSE
828  gaptol = gap * eps
829  ENDIF
830  isupmn = in
831  isupmx = 1
832 * Update WGAP so that it holds the minimum gap
833 * to the left or the right. This is crucial in the
834 * case where bisection is used to ensure that the
835 * eigenvalue is refined up to the required precision.
836 * The correct value is restored afterwards.
837  savgap = wgap(windex)
838  wgap(windex) = gap
839 * We want to use the Rayleigh Quotient Correction
840 * as often as possible since it converges quadratically
841 * when we are close enough to the desired eigenvalue.
842 * However, the Rayleigh Quotient can have the wrong sign
843 * and lead us away from the desired eigenvalue. In this
844 * case, the best we can do is to use bisection.
845  usedbs = .false.
846  usedrq = .false.
847 * Bisection is initially turned off unless it is forced
848  needbs = .NOT.tryrqc
849  120 CONTINUE
850 * Check if bisection should be used to refine eigenvalue
851  IF(needbs) THEN
852 * Take the bisection as new iterate
853  usedbs = .true.
854  itmp1 = iwork( iindr+windex )
855  offset = indexw( wbegin ) - 1
856  CALL dlarrb( in, d(ibegin),
857  $ work(indlld+ibegin-1),indeig,indeig,
858  $ zero, two*eps, offset,
859  $ work(wbegin),wgap(wbegin),
860  $ werr(wbegin),work( indwrk ),
861  $ iwork( iindwk ), pivmin, spdiam,
862  $ itmp1, iinfo )
863  IF( iinfo.NE.0 ) THEN
864  info = -3
865  RETURN
866  ENDIF
867  lambda = work( windex )
868 * Reset twist index from inaccurate LAMBDA to
869 * force computation of true MINGMA
870  iwork( iindr+windex ) = 0
871  ENDIF
872 * Given LAMBDA, compute the eigenvector.
873  CALL zlar1v( in, 1, in, lambda, d( ibegin ),
874  $ l( ibegin ), work(indld+ibegin-1),
875  $ work(indlld+ibegin-1),
876  $ pivmin, gaptol, z( ibegin, windex ),
877  $ .NOT.usedbs, negcnt, ztz, mingma,
878  $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
879  $ nrminv, resid, rqcorr, work( indwrk ) )
880  IF(iter .EQ. 0) THEN
881  bstres = resid
882  bstw = lambda
883  ELSEIF(resid.LT.bstres) THEN
884  bstres = resid
885  bstw = lambda
886  ENDIF
887  isupmn = min(isupmn,isuppz( 2*windex-1 ))
888  isupmx = max(isupmx,isuppz( 2*windex ))
889  iter = iter + 1
890 
891 * sin alpha <= |resid|/gap
892 * Note that both the residual and the gap are
893 * proportional to the matrix, so ||T|| doesn't play
894 * a role in the quotient
895 
896 *
897 * Convergence test for Rayleigh-Quotient iteration
898 * (omitted when Bisection has been used)
899 *
900  IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
901  $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
902  $ THEN
903 * We need to check that the RQCORR update doesn't
904 * move the eigenvalue away from the desired one and
905 * towards a neighbor. -> protection with bisection
906  IF(indeig.LE.negcnt) THEN
907 * The wanted eigenvalue lies to the left
908  sgndef = -one
909  ELSE
910 * The wanted eigenvalue lies to the right
911  sgndef = one
912  ENDIF
913 * We only use the RQCORR if it improves the
914 * the iterate reasonably.
915  IF( ( rqcorr*sgndef.GE.zero )
916  $ .AND.( lambda + rqcorr.LE. right)
917  $ .AND.( lambda + rqcorr.GE. left)
918  $ ) THEN
919  usedrq = .true.
920 * Store new midpoint of bisection interval in WORK
921  IF(sgndef.EQ.one) THEN
922 * The current LAMBDA is on the left of the true
923 * eigenvalue
924  left = lambda
925 * We prefer to assume that the error estimate
926 * is correct. We could make the interval not
927 * as a bracket but to be modified if the RQCORR
928 * chooses to. In this case, the RIGHT side should
929 * be modified as follows:
930 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
931  ELSE
932 * The current LAMBDA is on the right of the true
933 * eigenvalue
934  right = lambda
935 * See comment about assuming the error estimate is
936 * correct above.
937 * LEFT = MIN(LEFT, LAMBDA + RQCORR)
938  ENDIF
939  work( windex ) =
940  $ half * (right + left)
941 * Take RQCORR since it has the correct sign and
942 * improves the iterate reasonably
943  lambda = lambda + rqcorr
944 * Update width of error interval
945  werr( windex ) =
946  $ half * (right-left)
947  ELSE
948  needbs = .true.
949  ENDIF
950  IF(right-left.LT.rqtol*abs(lambda)) THEN
951 * The eigenvalue is computed to bisection accuracy
952 * compute eigenvector and stop
953  usedbs = .true.
954  goto 120
955  ELSEIF( iter.LT.maxitr ) THEN
956  goto 120
957  ELSEIF( iter.EQ.maxitr ) THEN
958  needbs = .true.
959  goto 120
960  ELSE
961  info = 5
962  RETURN
963  END IF
964  ELSE
965  stp2ii = .false.
966  IF(usedrq .AND. usedbs .AND.
967  $ bstres.LE.resid) THEN
968  lambda = bstw
969  stp2ii = .true.
970  ENDIF
971  IF (stp2ii) THEN
972 * improve error angle by second step
973  CALL zlar1v( in, 1, in, lambda,
974  $ d( ibegin ), l( ibegin ),
975  $ work(indld+ibegin-1),
976  $ work(indlld+ibegin-1),
977  $ pivmin, gaptol, z( ibegin, windex ),
978  $ .NOT.usedbs, negcnt, ztz, mingma,
979  $ iwork( iindr+windex ),
980  $ isuppz( 2*windex-1 ),
981  $ nrminv, resid, rqcorr, work( indwrk ) )
982  ENDIF
983  work( windex ) = lambda
984  END IF
985 *
986 * Compute FP-vector support w.r.t. whole matrix
987 *
988  isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
989  isuppz( 2*windex ) = isuppz( 2*windex )+oldien
990  zfrom = isuppz( 2*windex-1 )
991  zto = isuppz( 2*windex )
992  isupmn = isupmn + oldien
993  isupmx = isupmx + oldien
994 * Ensure vector is ok if support in the RQI has changed
995  IF(isupmn.LT.zfrom) THEN
996  DO 122 ii = isupmn,zfrom-1
997  z( ii, windex ) = zero
998  122 CONTINUE
999  ENDIF
1000  IF(isupmx.GT.zto) THEN
1001  DO 123 ii = zto+1,isupmx
1002  z( ii, windex ) = zero
1003  123 CONTINUE
1004  ENDIF
1005  CALL zdscal( zto-zfrom+1, nrminv,
1006  $ z( zfrom, windex ), 1 )
1007  125 CONTINUE
1008 * Update W
1009  w( windex ) = lambda+sigma
1010 * Recompute the gaps on the left and right
1011 * But only allow them to become larger and not
1012 * smaller (which can only happen through "bad"
1013 * cancellation and doesn't reflect the theory
1014 * where the initial gaps are underestimated due
1015 * to WERR being too crude.)
1016  IF(.NOT.eskip) THEN
1017  IF( k.GT.1) THEN
1018  wgap( windmn ) = max( wgap(windmn),
1019  $ w(windex)-werr(windex)
1020  $ - w(windmn)-werr(windmn) )
1021  ENDIF
1022  IF( windex.LT.wend ) THEN
1023  wgap( windex ) = max( savgap,
1024  $ w( windpl )-werr( windpl )
1025  $ - w( windex )-werr( windex) )
1026  ENDIF
1027  ENDIF
1028  idone = idone + 1
1029  ENDIF
1030 * here ends the code for the current child
1031 *
1032  139 CONTINUE
1033 * Proceed to any remaining child nodes
1034  newfst = j + 1
1035  140 CONTINUE
1036  150 CONTINUE
1037  ndepth = ndepth + 1
1038  go to 40
1039  END IF
1040  ibegin = iend + 1
1041  wbegin = wend + 1
1042  170 CONTINUE
1043 *
1044 
1045  RETURN
1046 *
1047 * End of ZLARRV
1048 *
1049  END
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:229
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:52
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:191
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:107
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:195
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
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:280
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:64