LAPACK 3.3.1
Linear Algebra PACKage

dlarrb.f

Go to the documentation of this file.
00001       SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
00002      $                   RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
00003      $                   PIVMIN, SPDIAM, TWIST, INFO )
00004 *
00005 *  -- LAPACK auxiliary routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       INTEGER            IFIRST, ILAST, INFO, N, OFFSET, TWIST
00012       DOUBLE PRECISION   PIVMIN, RTOL1, RTOL2, SPDIAM
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            IWORK( * )
00016       DOUBLE PRECISION   D( * ), LLD( * ), W( * ),
00017      $                   WERR( * ), WGAP( * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  Given the relatively robust representation(RRR) L D L^T, DLARRB
00024 *  does "limited" bisection to refine the eigenvalues of L D L^T,
00025 *  W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
00026 *  guesses for these eigenvalues are input in W, the corresponding estimate
00027 *  of the error in these guesses and their gaps are input in WERR
00028 *  and WGAP, respectively. During bisection, intervals
00029 *  [left, right] are maintained by storing their mid-points and
00030 *  semi-widths in the arrays W and WERR respectively.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  N       (input) INTEGER
00036 *          The order of the matrix.
00037 *
00038 *  D       (input) DOUBLE PRECISION array, dimension (N)
00039 *          The N diagonal elements of the diagonal matrix D.
00040 *
00041 *  LLD     (input) DOUBLE PRECISION array, dimension (N-1)
00042 *          The (N-1) elements L(i)*L(i)*D(i).
00043 *
00044 *  IFIRST  (input) INTEGER
00045 *          The index of the first eigenvalue to be computed.
00046 *
00047 *  ILAST   (input) INTEGER
00048 *          The index of the last eigenvalue to be computed.
00049 *
00050 *  RTOL1   (input) DOUBLE PRECISION
00051 *  RTOL2   (input) DOUBLE PRECISION
00052 *          Tolerance for the convergence of the bisection intervals.
00053 *          An interval [LEFT,RIGHT] has converged if
00054 *          RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
00055 *          where GAP is the (estimated) distance to the nearest
00056 *          eigenvalue.
00057 *
00058 *  OFFSET  (input) INTEGER
00059 *          Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
00060 *          through ILAST-OFFSET elements of these arrays are to be used.
00061 *
00062 *  W       (input/output) DOUBLE PRECISION array, dimension (N)
00063 *          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
00064 *          estimates of the eigenvalues of L D L^T indexed IFIRST throug
00065 *          ILAST.
00066 *          On output, these estimates are refined.
00067 *
00068 *  WGAP    (input/output) DOUBLE PRECISION array, dimension (N-1)
00069 *          On input, the (estimated) gaps between consecutive
00070 *          eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
00071 *          eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
00072 *          then WGAP(IFIRST-OFFSET) must be set to ZERO.
00073 *          On output, these gaps are refined.
00074 *
00075 *  WERR    (input/output) DOUBLE PRECISION array, dimension (N)
00076 *          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
00077 *          the errors in the estimates of the corresponding elements in W.
00078 *          On output, these errors are refined.
00079 *
00080 *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
00081 *          Workspace.
00082 *
00083 *  IWORK   (workspace) INTEGER array, dimension (2*N)
00084 *          Workspace.
00085 *
00086 *  PIVMIN  (input) DOUBLE PRECISION
00087 *          The minimum pivot in the Sturm sequence.
00088 *
00089 *  SPDIAM  (input) DOUBLE PRECISION
00090 *          The spectral diameter of the matrix.
00091 *
00092 *  TWIST   (input) INTEGER
00093 *          The twist index for the twisted factorization that is used
00094 *          for the negcount.
00095 *          TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
00096 *          TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
00097 *          TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
00098 *
00099 *  INFO    (output) INTEGER
00100 *          Error flag.
00101 *
00102 *  Further Details
00103 *  ===============
00104 *
00105 *  Based on contributions by
00106 *     Beresford Parlett, University of California, Berkeley, USA
00107 *     Jim Demmel, University of California, Berkeley, USA
00108 *     Inderjit Dhillon, University of Texas, Austin, USA
00109 *     Osni Marques, LBNL/NERSC, USA
00110 *     Christof Voemel, University of California, Berkeley, USA
00111 *
00112 *  =====================================================================
00113 *
00114 *     .. Parameters ..
00115       DOUBLE PRECISION   ZERO, TWO, HALF
00116       PARAMETER        ( ZERO = 0.0D0, TWO = 2.0D0,
00117      $                   HALF = 0.5D0 )
00118       INTEGER   MAXITR
00119 *     ..
00120 *     .. Local Scalars ..
00121       INTEGER            I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
00122      $                   OLNINT, PREV, R
00123       DOUBLE PRECISION   BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
00124      $                   RGAP, RIGHT, TMP, WIDTH
00125 *     ..
00126 *     .. External Functions ..
00127       INTEGER            DLANEG
00128       EXTERNAL           DLANEG
00129 *
00130 *     ..
00131 *     .. Intrinsic Functions ..
00132       INTRINSIC          ABS, MAX, MIN
00133 *     ..
00134 *     .. Executable Statements ..
00135 *
00136       INFO = 0
00137 *
00138       MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
00139      $           LOG( TWO ) ) + 2
00140       MNWDTH = TWO * PIVMIN
00141 *
00142       R = TWIST
00143       IF((R.LT.1).OR.(R.GT.N)) R = N
00144 *
00145 *     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
00146 *     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
00147 *     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
00148 *     for an unconverged interval is set to the index of the next unconverged
00149 *     interval, and is -1 or 0 for a converged interval. Thus a linked
00150 *     list of unconverged intervals is set up.
00151 *
00152       I1 = IFIRST
00153 *     The number of unconverged intervals
00154       NINT = 0
00155 *     The last unconverged interval found
00156       PREV = 0
00157 
00158       RGAP = WGAP( I1-OFFSET )
00159       DO 75 I = I1, ILAST
00160          K = 2*I
00161          II = I - OFFSET
00162          LEFT = W( II ) - WERR( II )
00163          RIGHT = W( II ) + WERR( II )
00164          LGAP = RGAP
00165          RGAP = WGAP( II )
00166          GAP = MIN( LGAP, RGAP )
00167 
00168 *        Make sure that [LEFT,RIGHT] contains the desired eigenvalue
00169 *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
00170 *
00171 *        Do while( NEGCNT(LEFT).GT.I-1 )
00172 *
00173          BACK = WERR( II )
00174  20      CONTINUE
00175          NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
00176          IF( NEGCNT.GT.I-1 ) THEN
00177             LEFT = LEFT - BACK
00178             BACK = TWO*BACK
00179             GO TO 20
00180          END IF
00181 *
00182 *        Do while( NEGCNT(RIGHT).LT.I )
00183 *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
00184 *
00185          BACK = WERR( II )
00186  50      CONTINUE
00187 
00188          NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
00189           IF( NEGCNT.LT.I ) THEN
00190              RIGHT = RIGHT + BACK
00191              BACK = TWO*BACK
00192              GO TO 50
00193           END IF
00194          WIDTH = HALF*ABS( LEFT - RIGHT )
00195          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
00196          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
00197          IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
00198 *           This interval has already converged and does not need refinement.
00199 *           (Note that the gaps might change through refining the
00200 *            eigenvalues, however, they can only get bigger.)
00201 *           Remove it from the list.
00202             IWORK( K-1 ) = -1
00203 *           Make sure that I1 always points to the first unconverged interval
00204             IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
00205             IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
00206          ELSE
00207 *           unconverged interval found
00208             PREV = I
00209             NINT = NINT + 1
00210             IWORK( K-1 ) = I + 1
00211             IWORK( K ) = NEGCNT
00212          END IF
00213          WORK( K-1 ) = LEFT
00214          WORK( K ) = RIGHT
00215  75   CONTINUE
00216 
00217 *
00218 *     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
00219 *     and while (ITER.LT.MAXITR)
00220 *
00221       ITER = 0
00222  80   CONTINUE
00223       PREV = I1 - 1
00224       I = I1
00225       OLNINT = NINT
00226 
00227       DO 100 IP = 1, OLNINT
00228          K = 2*I
00229          II = I - OFFSET
00230          RGAP = WGAP( II )
00231          LGAP = RGAP
00232          IF(II.GT.1) LGAP = WGAP( II-1 )
00233          GAP = MIN( LGAP, RGAP )
00234          NEXT = IWORK( K-1 )
00235          LEFT = WORK( K-1 )
00236          RIGHT = WORK( K )
00237          MID = HALF*( LEFT + RIGHT )
00238 
00239 *        semiwidth of interval
00240          WIDTH = RIGHT - MID
00241          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
00242          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
00243          IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
00244      $       ( ITER.EQ.MAXITR ) )THEN
00245 *           reduce number of unconverged intervals
00246             NINT = NINT - 1
00247 *           Mark interval as converged.
00248             IWORK( K-1 ) = 0
00249             IF( I1.EQ.I ) THEN
00250                I1 = NEXT
00251             ELSE
00252 *              Prev holds the last unconverged interval previously examined
00253                IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
00254             END IF
00255             I = NEXT
00256             GO TO 100
00257          END IF
00258          PREV = I
00259 *
00260 *        Perform one bisection step
00261 *
00262          NEGCNT = DLANEG( N, D, LLD, MID, PIVMIN, R )
00263          IF( NEGCNT.LE.I-1 ) THEN
00264             WORK( K-1 ) = MID
00265          ELSE
00266             WORK( K ) = MID
00267          END IF
00268          I = NEXT
00269  100  CONTINUE
00270       ITER = ITER + 1
00271 *     do another loop if there are still unconverged intervals
00272 *     However, in the last iteration, all intervals are accepted
00273 *     since this is the best we can do.
00274       IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
00275 *
00276 *
00277 *     At this point, all the intervals have converged
00278       DO 110 I = IFIRST, ILAST
00279          K = 2*I
00280          II = I - OFFSET
00281 *        All intervals marked by '0' have been refined.
00282          IF( IWORK( K-1 ).EQ.0 ) THEN
00283             W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
00284             WERR( II ) = WORK( K ) - W( II )
00285          END IF
00286  110  CONTINUE
00287 *
00288       DO 111 I = IFIRST+1, ILAST
00289          K = 2*I
00290          II = I - OFFSET
00291          WGAP( II-1 ) = MAX( ZERO,
00292      $                     W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))
00293  111  CONTINUE
00294 
00295       RETURN
00296 *
00297 *     End of DLARRB
00298 *
00299       END
 All Files Functions