LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slarrv()

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
!>          Lower bound 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]VU
!>          VU is REAL
!>          Upper bound of the interval that contains the desired
!>          eigenvalues. VL < VU. 
!>          Note: VU is currently not used by this implementation of SLARRV, VU is
!>          passed to SLARRV because it could be used compute gaps on the right end
!>          of the extremal eigenvalues. However, with not much initial accuracy in
!>          LAMBDA and VU, the formula can lead to an overestimation of the right gap
!>          and thus to inadequately early RQI 'convergence'. This is currently
!>          prevented this by forcing a small right gap. And so it turns out that VU
!>          is currently not used by this implementation of SLARRV.
!>
[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 split.) 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 < 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 occurred 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.
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.

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