LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dlarrv()

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

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

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

Purpose:
 DLARRV 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 DLARRE.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in]VL
          VL is DOUBLE PRECISION
          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 DOUBLE PRECISION
          Upper bound of the interval that contains the desired
          eigenvalues. VL < VU. 
          Note: VU is currently not used by this implementation of DLARRV, VU is
          passed to DLARRV 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 DLARRV.
[in,out]D
          D is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLARRE.
          On exit, L is overwritten.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          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 DOUBLE PRECISION
[in]RTOL1
          RTOL1 is DOUBLE PRECISION
[in]RTOL2
          RTOL2 is DOUBLE PRECISION
           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 DOUBLE PRECISION 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 DLARRE 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLARRV.
          < 0:  One of the called subroutines signaled an internal problem.
                Needs inspection of the corresponding parameter IINFO
                for further information.

          =-1:  Problem in DLARRB when refining a child's eigenvalues.
          =-2:  Problem in DLARRF 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 DLARRB 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
June 2016
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 294 of file dlarrv.f.

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