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