LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlarrj.f
Go to the documentation of this file.
00001 *> \brief \b DLARRJ
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLARRJ + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrj.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrj.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrj.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
00022 *                          RTOL, OFFSET, W, WERR, WORK, IWORK,
00023 *                          PIVMIN, SPDIAM, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       INTEGER            IFIRST, ILAST, INFO, N, OFFSET
00027 *       DOUBLE PRECISION   PIVMIN, RTOL, SPDIAM
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IWORK( * )
00031 *       DOUBLE PRECISION   D( * ), E2( * ), W( * ),
00032 *      $                   WERR( * ), WORK( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> Given the initial eigenvalue approximations of T, DLARRJ
00042 *> does  bisection to refine the eigenvalues of T,
00043 *> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
00044 *> guesses for these eigenvalues are input in W, the corresponding estimate
00045 *> of the error in these guesses in WERR. During bisection, intervals
00046 *> [left, right] are maintained by storing their mid-points and
00047 *> semi-widths in the arrays W and WERR respectively.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] N
00054 *> \verbatim
00055 *>          N is INTEGER
00056 *>          The order of the matrix.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] D
00060 *> \verbatim
00061 *>          D is DOUBLE PRECISION array, dimension (N)
00062 *>          The N diagonal elements of T.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] E2
00066 *> \verbatim
00067 *>          E2 is DOUBLE PRECISION array, dimension (N-1)
00068 *>          The Squares of the (N-1) subdiagonal elements of T.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] IFIRST
00072 *> \verbatim
00073 *>          IFIRST is INTEGER
00074 *>          The index of the first eigenvalue to be computed.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] ILAST
00078 *> \verbatim
00079 *>          ILAST is INTEGER
00080 *>          The index of the last eigenvalue to be computed.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] RTOL
00084 *> \verbatim
00085 *>          RTOL is DOUBLE PRECISION
00086 *>          Tolerance for the convergence of the bisection intervals.
00087 *>          An interval [LEFT,RIGHT] has converged if
00088 *>          RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|).
00089 *> \endverbatim
00090 *>
00091 *> \param[in] OFFSET
00092 *> \verbatim
00093 *>          OFFSET is INTEGER
00094 *>          Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
00095 *>          through ILAST-OFFSET elements of these arrays are to be used.
00096 *> \endverbatim
00097 *>
00098 *> \param[in,out] W
00099 *> \verbatim
00100 *>          W is DOUBLE PRECISION array, dimension (N)
00101 *>          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
00102 *>          estimates of the eigenvalues of L D L^T indexed IFIRST through
00103 *>          ILAST.
00104 *>          On output, these estimates are refined.
00105 *> \endverbatim
00106 *>
00107 *> \param[in,out] WERR
00108 *> \verbatim
00109 *>          WERR is DOUBLE PRECISION array, dimension (N)
00110 *>          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
00111 *>          the errors in the estimates of the corresponding elements in W.
00112 *>          On output, these errors are refined.
00113 *> \endverbatim
00114 *>
00115 *> \param[out] WORK
00116 *> \verbatim
00117 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
00118 *>          Workspace.
00119 *> \endverbatim
00120 *>
00121 *> \param[out] IWORK
00122 *> \verbatim
00123 *>          IWORK is INTEGER array, dimension (2*N)
00124 *>          Workspace.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] PIVMIN
00128 *> \verbatim
00129 *>          PIVMIN is DOUBLE PRECISION
00130 *>          The minimum pivot in the Sturm sequence for T.
00131 *> \endverbatim
00132 *>
00133 *> \param[in] SPDIAM
00134 *> \verbatim
00135 *>          SPDIAM is DOUBLE PRECISION
00136 *>          The spectral diameter of T.
00137 *> \endverbatim
00138 *>
00139 *> \param[out] INFO
00140 *> \verbatim
00141 *>          INFO is INTEGER
00142 *>          Error flag.
00143 *> \endverbatim
00144 *
00145 *  Authors:
00146 *  ========
00147 *
00148 *> \author Univ. of Tennessee 
00149 *> \author Univ. of California Berkeley 
00150 *> \author Univ. of Colorado Denver 
00151 *> \author NAG Ltd. 
00152 *
00153 *> \date November 2011
00154 *
00155 *> \ingroup auxOTHERauxiliary
00156 *
00157 *> \par Contributors:
00158 *  ==================
00159 *>
00160 *> Beresford Parlett, University of California, Berkeley, USA \n
00161 *> Jim Demmel, University of California, Berkeley, USA \n
00162 *> Inderjit Dhillon, University of Texas, Austin, USA \n
00163 *> Osni Marques, LBNL/NERSC, USA \n
00164 *> Christof Voemel, University of California, Berkeley, USA
00165 *
00166 *  =====================================================================
00167       SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
00168      $                   RTOL, OFFSET, W, WERR, WORK, IWORK,
00169      $                   PIVMIN, SPDIAM, INFO )
00170 *
00171 *  -- LAPACK auxiliary routine (version 3.4.0) --
00172 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00173 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00174 *     November 2011
00175 *
00176 *     .. Scalar Arguments ..
00177       INTEGER            IFIRST, ILAST, INFO, N, OFFSET
00178       DOUBLE PRECISION   PIVMIN, RTOL, SPDIAM
00179 *     ..
00180 *     .. Array Arguments ..
00181       INTEGER            IWORK( * )
00182       DOUBLE PRECISION   D( * ), E2( * ), W( * ),
00183      $                   WERR( * ), WORK( * )
00184 *     ..
00185 *
00186 *  =====================================================================
00187 *
00188 *     .. Parameters ..
00189       DOUBLE PRECISION   ZERO, ONE, TWO, HALF
00190       PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00191      $                   HALF = 0.5D0 )
00192       INTEGER   MAXITR
00193 *     ..
00194 *     .. Local Scalars ..
00195       INTEGER            CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
00196      $                   OLNINT, P, PREV, SAVI1
00197       DOUBLE PRECISION   DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
00198 *
00199 *     ..
00200 *     .. Intrinsic Functions ..
00201       INTRINSIC          ABS, MAX
00202 *     ..
00203 *     .. Executable Statements ..
00204 *
00205       INFO = 0
00206 *
00207       MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
00208      $           LOG( TWO ) ) + 2
00209 *
00210 *     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
00211 *     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
00212 *     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
00213 *     for an unconverged interval is set to the index of the next unconverged
00214 *     interval, and is -1 or 0 for a converged interval. Thus a linked
00215 *     list of unconverged intervals is set up.
00216 *
00217 
00218       I1 = IFIRST
00219       I2 = ILAST
00220 *     The number of unconverged intervals
00221       NINT = 0
00222 *     The last unconverged interval found
00223       PREV = 0
00224       DO 75 I = I1, I2
00225          K = 2*I
00226          II = I - OFFSET
00227          LEFT = W( II ) - WERR( II )
00228          MID = W(II)
00229          RIGHT = W( II ) + WERR( II )
00230          WIDTH = RIGHT - MID
00231          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
00232 
00233 *        The following test prevents the test of converged intervals
00234          IF( WIDTH.LT.RTOL*TMP ) THEN
00235 *           This interval has already converged and does not need refinement.
00236 *           (Note that the gaps might change through refining the
00237 *            eigenvalues, however, they can only get bigger.)
00238 *           Remove it from the list.
00239             IWORK( K-1 ) = -1
00240 *           Make sure that I1 always points to the first unconverged interval
00241             IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1
00242             IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1
00243          ELSE
00244 *           unconverged interval found
00245             PREV = I
00246 *           Make sure that [LEFT,RIGHT] contains the desired eigenvalue
00247 *
00248 *           Do while( CNT(LEFT).GT.I-1 )
00249 *
00250             FAC = ONE
00251  20         CONTINUE
00252             CNT = 0
00253             S = LEFT
00254             DPLUS = D( 1 ) - S
00255             IF( DPLUS.LT.ZERO ) CNT = CNT + 1
00256             DO 30 J = 2, N
00257                DPLUS = D( J ) - S - E2( J-1 )/DPLUS
00258                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
00259  30         CONTINUE
00260             IF( CNT.GT.I-1 ) THEN
00261                LEFT = LEFT - WERR( II )*FAC
00262                FAC = TWO*FAC
00263                GO TO 20
00264             END IF
00265 *
00266 *           Do while( CNT(RIGHT).LT.I )
00267 *
00268             FAC = ONE
00269  50         CONTINUE
00270             CNT = 0
00271             S = RIGHT
00272             DPLUS = D( 1 ) - S
00273             IF( DPLUS.LT.ZERO ) CNT = CNT + 1
00274             DO 60 J = 2, N
00275                DPLUS = D( J ) - S - E2( J-1 )/DPLUS
00276                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
00277  60         CONTINUE
00278             IF( CNT.LT.I ) THEN
00279                RIGHT = RIGHT + WERR( II )*FAC
00280                FAC = TWO*FAC
00281                GO TO 50
00282             END IF
00283             NINT = NINT + 1
00284             IWORK( K-1 ) = I + 1
00285             IWORK( K ) = CNT
00286          END IF
00287          WORK( K-1 ) = LEFT
00288          WORK( K ) = RIGHT
00289  75   CONTINUE
00290 
00291 
00292       SAVI1 = I1
00293 *
00294 *     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
00295 *     and while (ITER.LT.MAXITR)
00296 *
00297       ITER = 0
00298  80   CONTINUE
00299       PREV = I1 - 1
00300       I = I1
00301       OLNINT = NINT
00302 
00303       DO 100 P = 1, OLNINT
00304          K = 2*I
00305          II = I - OFFSET
00306          NEXT = IWORK( K-1 )
00307          LEFT = WORK( K-1 )
00308          RIGHT = WORK( K )
00309          MID = HALF*( LEFT + RIGHT )
00310 
00311 *        semiwidth of interval
00312          WIDTH = RIGHT - MID
00313          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
00314 
00315          IF( ( WIDTH.LT.RTOL*TMP ) .OR.
00316      $      (ITER.EQ.MAXITR) )THEN
00317 *           reduce number of unconverged intervals
00318             NINT = NINT - 1
00319 *           Mark interval as converged.
00320             IWORK( K-1 ) = 0
00321             IF( I1.EQ.I ) THEN
00322                I1 = NEXT
00323             ELSE
00324 *              Prev holds the last unconverged interval previously examined
00325                IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
00326             END IF
00327             I = NEXT
00328             GO TO 100
00329          END IF
00330          PREV = I
00331 *
00332 *        Perform one bisection step
00333 *
00334          CNT = 0
00335          S = MID
00336          DPLUS = D( 1 ) - S
00337          IF( DPLUS.LT.ZERO ) CNT = CNT + 1
00338          DO 90 J = 2, N
00339             DPLUS = D( J ) - S - E2( J-1 )/DPLUS
00340             IF( DPLUS.LT.ZERO ) CNT = CNT + 1
00341  90      CONTINUE
00342          IF( CNT.LE.I-1 ) THEN
00343             WORK( K-1 ) = MID
00344          ELSE
00345             WORK( K ) = MID
00346          END IF
00347          I = NEXT
00348 
00349  100  CONTINUE
00350       ITER = ITER + 1
00351 *     do another loop if there are still unconverged intervals
00352 *     However, in the last iteration, all intervals are accepted
00353 *     since this is the best we can do.
00354       IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
00355 *
00356 *
00357 *     At this point, all the intervals have converged
00358       DO 110 I = SAVI1, ILAST
00359          K = 2*I
00360          II = I - OFFSET
00361 *        All intervals marked by '0' have been refined.
00362          IF( IWORK( K-1 ).EQ.0 ) THEN
00363             W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
00364             WERR( II ) = WORK( K ) - W( II )
00365          END IF
00366  110  CONTINUE
00367 *
00368 
00369       RETURN
00370 *
00371 *     End of DLARRJ
00372 *
00373       END
 All Files Functions