Return to dlarrb.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b DLARRB provides limited bisection to locate eigenvalues for more accuracy.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLARRB + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrb.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrb.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrb.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
   22: *                          RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
   23: *                          PIVMIN, SPDIAM, TWIST, INFO )
   24: * 
   25: *       .. Scalar Arguments ..
   26: *       INTEGER            IFIRST, ILAST, INFO, N, OFFSET, TWIST
   27: *       DOUBLE PRECISION   PIVMIN, RTOL1, RTOL2, SPDIAM
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       INTEGER            IWORK( * )
   31: *       DOUBLE PRECISION   D( * ), LLD( * ), W( * ),
   32: *      $                   WERR( * ), WGAP( * ), WORK( * )
   33: *       ..
   34: *  
   35: *
   36: *> \par Purpose:
   37: *  =============
   38: *>
   39: *> \verbatim
   40: *>
   41: *> Given the relatively robust representation(RRR) L D L^T, DLARRB
   42: *> does "limited" bisection to refine the eigenvalues of L D L^T,
   43: *> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
   44: *> guesses for these eigenvalues are input in W, the corresponding estimate
   45: *> of the error in these guesses and their gaps are input in WERR
   46: *> and WGAP, respectively. During bisection, intervals
   47: *> [left, right] are maintained by storing their mid-points and
   48: *> semi-widths in the arrays W and WERR respectively.
   49: *> \endverbatim
   50: *
   51: *  Arguments:
   52: *  ==========
   53: *
   54: *> \param[in] N
   55: *> \verbatim
   56: *>          N is INTEGER
   57: *>          The order of the matrix.
   58: *> \endverbatim
   59: *>
   60: *> \param[in] D
   61: *> \verbatim
   62: *>          D is DOUBLE PRECISION array, dimension (N)
   63: *>          The N diagonal elements of the diagonal matrix D.
   64: *> \endverbatim
   65: *>
   66: *> \param[in] LLD
   67: *> \verbatim
   68: *>          LLD is DOUBLE PRECISION array, dimension (N-1)
   69: *>          The (N-1) elements L(i)*L(i)*D(i).
   70: *> \endverbatim
   71: *>
   72: *> \param[in] IFIRST
   73: *> \verbatim
   74: *>          IFIRST is INTEGER
   75: *>          The index of the first eigenvalue to be computed.
   76: *> \endverbatim
   77: *>
   78: *> \param[in] ILAST
   79: *> \verbatim
   80: *>          ILAST is INTEGER
   81: *>          The index of the last eigenvalue to be computed.
   82: *> \endverbatim
   83: *>
   84: *> \param[in] RTOL1
   85: *> \verbatim
   86: *>          RTOL1 is DOUBLE PRECISION
   87: *> \endverbatim
   88: *>
   89: *> \param[in] RTOL2
   90: *> \verbatim
   91: *>          RTOL2 is DOUBLE PRECISION
   92: *>          Tolerance for the convergence of the bisection intervals.
   93: *>          An interval [LEFT,RIGHT] has converged if
   94: *>          RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
   95: *>          where GAP is the (estimated) distance to the nearest
   96: *>          eigenvalue.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] OFFSET
  100: *> \verbatim
  101: *>          OFFSET is INTEGER
  102: *>          Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
  103: *>          through ILAST-OFFSET elements of these arrays are to be used.
  104: *> \endverbatim
  105: *>
  106: *> \param[in,out] W
  107: *> \verbatim
  108: *>          W is DOUBLE PRECISION array, dimension (N)
  109: *>          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
  110: *>          estimates of the eigenvalues of L D L^T indexed IFIRST throug
  111: *>          ILAST.
  112: *>          On output, these estimates are refined.
  113: *> \endverbatim
  114: *>
  115: *> \param[in,out] WGAP
  116: *> \verbatim
  117: *>          WGAP is DOUBLE PRECISION array, dimension (N-1)
  118: *>          On input, the (estimated) gaps between consecutive
  119: *>          eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
  120: *>          eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
  121: *>          then WGAP(IFIRST-OFFSET) must be set to ZERO.
  122: *>          On output, these gaps are refined.
  123: *> \endverbatim
  124: *>
  125: *> \param[in,out] WERR
  126: *> \verbatim
  127: *>          WERR is DOUBLE PRECISION array, dimension (N)
  128: *>          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
  129: *>          the errors in the estimates of the corresponding elements in W.
  130: *>          On output, these errors are refined.
  131: *> \endverbatim
  132: *>
  133: *> \param[out] WORK
  134: *> \verbatim
  135: *>          WORK is DOUBLE PRECISION array, dimension (2*N)
  136: *>          Workspace.
  137: *> \endverbatim
  138: *>
  139: *> \param[out] IWORK
  140: *> \verbatim
  141: *>          IWORK is INTEGER array, dimension (2*N)
  142: *>          Workspace.
  143: *> \endverbatim
  144: *>
  145: *> \param[in] PIVMIN
  146: *> \verbatim
  147: *>          PIVMIN is DOUBLE PRECISION
  148: *>          The minimum pivot in the Sturm sequence.
  149: *> \endverbatim
  150: *>
  151: *> \param[in] SPDIAM
  152: *> \verbatim
  153: *>          SPDIAM is DOUBLE PRECISION
  154: *>          The spectral diameter of the matrix.
  155: *> \endverbatim
  156: *>
  157: *> \param[in] TWIST
  158: *> \verbatim
  159: *>          TWIST is INTEGER
  160: *>          The twist index for the twisted factorization that is used
  161: *>          for the negcount.
  162: *>          TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
  163: *>          TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
  164: *>          TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
  165: *> \endverbatim
  166: *>
  167: *> \param[out] INFO
  168: *> \verbatim
  169: *>          INFO is INTEGER
  170: *>          Error flag.
  171: *> \endverbatim
  172: *
  173: *  Authors:
  174: *  ========
  175: *
  176: *> \author Univ. of Tennessee 
  177: *> \author Univ. of California Berkeley 
  178: *> \author Univ. of Colorado Denver 
  179: *> \author NAG Ltd. 
  180: *
  181: *> \date September 2012
  182: *
  183: *> \ingroup auxOTHERauxiliary
  184: *
  185: *> \par Contributors:
  186: *  ==================
  187: *>
  188: *> Beresford Parlett, University of California, Berkeley, USA \n
  189: *> Jim Demmel, University of California, Berkeley, USA \n
  190: *> Inderjit Dhillon, University of Texas, Austin, USA \n
  191: *> Osni Marques, LBNL/NERSC, USA \n
  192: *> Christof Voemel, University of California, Berkeley, USA
  193: *
  194: *  =====================================================================
  195:       SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
  196:      $                   RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
  197:      $                   PIVMIN, SPDIAM, TWIST, INFO )
  198: *
  199: *  -- LAPACK auxiliary routine (version 3.4.2) --
  200: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  201: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  202: *     September 2012
  203: *
  204: *     .. Scalar Arguments ..
  205:       INTEGER            IFIRST, ILAST, INFO, N, OFFSET, TWIST
  206:       DOUBLE PRECISION   PIVMIN, RTOL1, RTOL2, SPDIAM
  207: *     ..
  208: *     .. Array Arguments ..
  209:       INTEGER            IWORK( * )
  210:       DOUBLE PRECISION   D( * ), LLD( * ), W( * ),
  211:      $                   WERR( * ), WGAP( * ), WORK( * )
  212: *     ..
  213: *
  214: *  =====================================================================
  215: *
  216: *     .. Parameters ..
  217:       DOUBLE PRECISION   ZERO, TWO, HALF
  218:       PARAMETER        ( ZERO = 0.0D0, TWO = 2.0D0,
  219:      $                   HALF = 0.5D0 )
  220:       INTEGER   MAXITR
  221: *     ..
  222: *     .. Local Scalars ..
  223:       INTEGER            I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
  224:      $                   OLNINT, PREV, R
  225:       DOUBLE PRECISION   BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
  226:      $                   RGAP, RIGHT, TMP, WIDTH
  227: *     ..
  228: *     .. External Functions ..
  229:       INTEGER            DLANEG
  230:       EXTERNAL           DLANEG
  231: *
  232: *     ..
  233: *     .. Intrinsic Functions ..
  234:       INTRINSIC          ABS, MAX, MIN
  235: *     ..
  236: *     .. Executable Statements ..
  237: *
  238:       INFO = 0
  239: *
  240:       MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
  241:      $           LOG( TWO ) ) + 2
  242:       MNWDTH = TWO * PIVMIN
  243: *
  244:       R = TWIST
  245:       IF((R.LT.1).OR.(R.GT.N)) R = N
  246: *
  247: *     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
  248: *     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
  249: *     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
  250: *     for an unconverged interval is set to the index of the next unconverged
  251: *     interval, and is -1 or 0 for a converged interval. Thus a linked
  252: *     list of unconverged intervals is set up.
  253: *
  254:       I1 = IFIRST
  255: *     The number of unconverged intervals
  256:       NINT = 0
  257: *     The last unconverged interval found
  258:       PREV = 0
  259: 
  260:       RGAP = WGAP( I1-OFFSET )
  261:       DO 75 I = I1, ILAST
  262:          K = 2*I
  263:          II = I - OFFSET
  264:          LEFT = W( II ) - WERR( II )
  265:          RIGHT = W( II ) + WERR( II )
  266:          LGAP = RGAP
  267:          RGAP = WGAP( II )
  268:          GAP = MIN( LGAP, RGAP )
  269: 
  270: *        Make sure that [LEFT,RIGHT] contains the desired eigenvalue
  271: *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
  272: *
  273: *        Do while( NEGCNT(LEFT).GT.I-1 )
  274: *
  275:          BACK = WERR( II )
  276:  20      CONTINUE
  277:          NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
  278:          IF( NEGCNT.GT.I-1 ) THEN
  279:             LEFT = LEFT - BACK
  280:             BACK = TWO*BACK
  281:             GO TO 20
  282:          END IF
  283: *
  284: *        Do while( NEGCNT(RIGHT).LT.I )
  285: *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
  286: *
  287:          BACK = WERR( II )
  288:  50      CONTINUE
  289: 
  290:          NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
  291:           IF( NEGCNT.LT.I ) THEN
  292:              RIGHT = RIGHT + BACK
  293:              BACK = TWO*BACK
  294:              GO TO 50
  295:           END IF
  296:          WIDTH = HALF*ABS( LEFT - RIGHT )
  297:          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
  298:          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
  299:          IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
  300: *           This interval has already converged and does not need refinement.
  301: *           (Note that the gaps might change through refining the
  302: *            eigenvalues, however, they can only get bigger.)
  303: *           Remove it from the list.
  304:             IWORK( K-1 ) = -1
  305: *           Make sure that I1 always points to the first unconverged interval
  306:             IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
  307:             IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
  308:          ELSE
  309: *           unconverged interval found
  310:             PREV = I
  311:             NINT = NINT + 1
  312:             IWORK( K-1 ) = I + 1
  313:             IWORK( K ) = NEGCNT
  314:          END IF
  315:          WORK( K-1 ) = LEFT
  316:          WORK( K ) = RIGHT
  317:  75   CONTINUE
  318: 
  319: *
  320: *     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
  321: *     and while (ITER.LT.MAXITR)
  322: *
  323:       ITER = 0
  324:  80   CONTINUE
  325:       PREV = I1 - 1
  326:       I = I1
  327:       OLNINT = NINT
  328: 
  329:       DO 100 IP = 1, OLNINT
  330:          K = 2*I
  331:          II = I - OFFSET
  332:          RGAP = WGAP( II )
  333:          LGAP = RGAP
  334:          IF(II.GT.1) LGAP = WGAP( II-1 )
  335:          GAP = MIN( LGAP, RGAP )
  336:          NEXT = IWORK( K-1 )
  337:          LEFT = WORK( K-1 )
  338:          RIGHT = WORK( K )
  339:          MID = HALF*( LEFT + RIGHT )
  340: 
  341: *        semiwidth of interval
  342:          WIDTH = RIGHT - MID
  343:          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
  344:          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
  345:          IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
  346:      $       ( ITER.EQ.MAXITR ) )THEN
  347: *           reduce number of unconverged intervals
  348:             NINT = NINT - 1
  349: *           Mark interval as converged.
  350:             IWORK( K-1 ) = 0
  351:             IF( I1.EQ.I ) THEN
  352:                I1 = NEXT
  353:             ELSE
  354: *              Prev holds the last unconverged interval previously examined
  355:                IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
  356:             END IF
  357:             I = NEXT
  358:             GO TO 100
  359:          END IF
  360:          PREV = I
  361: *
  362: *        Perform one bisection step
  363: *
  364:          NEGCNT = DLANEG( N, D, LLD, MID, PIVMIN, R )
  365:          IF( NEGCNT.LE.I-1 ) THEN
  366:             WORK( K-1 ) = MID
  367:          ELSE
  368:             WORK( K ) = MID
  369:          END IF
  370:          I = NEXT
  371:  100  CONTINUE
  372:       ITER = ITER + 1
  373: *     do another loop if there are still unconverged intervals
  374: *     However, in the last iteration, all intervals are accepted
  375: *     since this is the best we can do.
  376:       IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
  377: *
  378: *
  379: *     At this point, all the intervals have converged
  380:       DO 110 I = IFIRST, ILAST
  381:          K = 2*I
  382:          II = I - OFFSET
  383: *        All intervals marked by '0' have been refined.
  384:          IF( IWORK( K-1 ).EQ.0 ) THEN
  385:             W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
  386:             WERR( II ) = WORK( K ) - W( II )
  387:          END IF
  388:  110  CONTINUE
  389: *
  390:       DO 111 I = IFIRST+1, ILAST
  391:          K = 2*I
  392:          II = I - OFFSET
  393:          WGAP( II-1 ) = MAX( ZERO,
  394:      $                     W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))
  395:  111  CONTINUE
  396: 
  397:       RETURN
  398: *
  399: *     End of DLARRB
  400: *
  401:       END