LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
slarrv.f File Reference

Go to the source code of this file.

Functions/Subroutines

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 of L D LT. More...
 

Function/Subroutine Documentation

subroutine slarrv ( integer  N,
real  VL,
real  VU,
real, dimension( * )  D,
real, dimension( * )  L,
real  PIVMIN,
integer, dimension( * )  ISPLIT,
integer  M,
integer  DOL,
integer  DOU,
real  MINRGP,
real  RTOL1,
real  RTOL2,
real, dimension( * )  W,
real, dimension( * )  WERR,
real, dimension( * )  WGAP,
integer, dimension( * )  IBLOCK,
integer, dimension( * )  INDEXW,
real, dimension( * )  GERS,
real, dimension( ldz, * )  Z,
integer  LDZ,
integer, dimension( * )  ISUPPZ,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.

Download SLARRV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SLARRV computes the eigenvectors of the tridiagonal matrix
 T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
 The input eigenvalues should have been computed by SLARRE.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in]VL
          VL is REAL
[in]VU
          VU is REAL
          Lower and upper bounds of the interval that contains the desired
          eigenvalues. VL < VU. Needed to compute gaps on the left or right
          end of the extremal eigenvalues in the desired RANGE.
[in,out]D
          D is REAL array, dimension (N)
          On entry, the N diagonal elements of the diagonal matrix D.
          On exit, D may be overwritten.
[in,out]L
          L is REAL array, dimension (N)
          On entry, the (N-1) subdiagonal elements of the unit
          bidiagonal matrix L are in elements 1 to N-1 of L
          (if the matrix is not splitted.) At the end of each block
          is stored the corresponding shift as given by SLARRE.
          On exit, L is overwritten.
[in]PIVMIN
          PIVMIN is REAL
          The minimum pivot allowed in the Sturm sequence.
[in]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into blocks.
          The first block consists of rows/columns 1 to
          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
          through ISPLIT( 2 ), etc.
[in]M
          M is INTEGER
          The total number of input eigenvalues.  0 <= M <= N.
[in]DOL
          DOL is INTEGER
[in]DOU
          DOU is INTEGER
          If the user wants to compute only selected eigenvectors from all
          the eigenvalues supplied, he can specify an index range DOL:DOU.
          Or else the setting DOL=1, DOU=M should be applied.
          Note that DOL and DOU refer to the order in which the eigenvalues
          are stored in W.
          If the user wants to compute only selected eigenpairs, then
          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
          computed eigenvectors. All other columns of Z are set to zero.
[in]MINRGP
          MINRGP is REAL
[in]RTOL1
          RTOL1 is REAL
[in]RTOL2
          RTOL2 is REAL
           Parameters for bisection.
           An interval [LEFT,RIGHT] has converged if
           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
[in,out]W
          W is REAL array, dimension (N)
          The first M elements of W contain the APPROXIMATE eigenvalues for
          which eigenvectors are to be computed.  The eigenvalues
          should be grouped by split-off block and ordered from
          smallest to largest within the block ( The output array
          W from SLARRE is expected here ). Furthermore, they are with
          respect to the shift of the corresponding root representation
          for their block. On exit, W holds the eigenvalues of the
          UNshifted matrix.
[in,out]WERR
          WERR is REAL array, dimension (N)
          The first M elements contain the semiwidth of the uncertainty
          interval of the corresponding eigenvalue in W
[in,out]WGAP
          WGAP is REAL array, dimension (N)
          The separation from the right neighbor eigenvalue in W.
[in]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          The indices of the blocks (submatrices) associated with the
          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
          W(i) belongs to the first block from the top, =2 if W(i)
          belongs to the second block, etc.
[in]INDEXW
          INDEXW is INTEGER array, dimension (N)
          The indices of the eigenvalues within each block (submatrix);
          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
[in]GERS
          GERS is REAL array, dimension (2*N)
          The N Gerschgorin intervals (the i-th Gerschgorin interval
          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
          be computed from the original UNshifted matrix.
[out]Z
          Z is REAL array, dimension (LDZ, max(1,M) )
          If INFO = 0, the first M columns of Z contain the
          orthonormal eigenvectors of the matrix T
          corresponding to the input eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).
[out]ISUPPZ
          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
          The support of the eigenvectors in Z, i.e., the indices
          indicating the nonzero elements in Z. The I-th eigenvector
          is nonzero only in elements ISUPPZ( 2*I-1 ) through
          ISUPPZ( 2*I ).
[out]WORK
          WORK is REAL array, dimension (12*N)
[out]IWORK
          IWORK is INTEGER array, dimension (7*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit

          > 0:  A problem occured in SLARRV.
          < 0:  One of the called subroutines signaled an internal problem.
                Needs inspection of the corresponding parameter IINFO
                for further information.

          =-1:  Problem in SLARRB when refining a child's eigenvalues.
          =-2:  Problem in SLARRF when computing the RRR of a child.
                When a child is inside a tight cluster, it can be difficult
                to find an RRR. A partial remedy from the user's point of
                view is to make the parameter MINRGP smaller and recompile.
                However, as the orthogonality of the computed vectors is
                proportional to 1/MINRGP, the user should be aware that
                he might be trading in precision when he decreases MINRGP.
          =-3:  Problem in SLARRB when refining a single eigenvalue
                after the Rayleigh correction was rejected.
          = 5:  The Rayleigh Quotient Iteration failed to converge to
                full accuracy in MAXITR steps.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 285 of file slarrv.f.

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 *
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:112
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:195
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:198
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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:232
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: