Annotation of rpl/lapack/lapack/dlarrj.f, revision 1.18

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

CVSweb interface <joel.bertrand@systella.fr>