slarre.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 777 行 · 第 1/4 页

</span><span class="comment">*</span><span class="comment">     Christof Voemel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
     $                   MAXGROWTH, ONE, PERT, TWO, ZERO
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0,
     $                     TWO = 2.0E0, FOUR=4.0E0,
     $                     HNDRD = 100.0E0,
     $                     PERT = 4.0E0,
     $                     HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
     $                     MAXGROWTH = 64.0E0, FUDGE = 2.0E0 )
      INTEGER            MAXTRY, ALLRNG, INDRNG, VALRNG
      PARAMETER          ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
     $                     VALRNG = 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            FORCEB, NOREP, USEDQD
      INTEGER            CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
     $                   IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
     $                   WBEGIN, WEND
      REAL               AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
     $                   EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
     $                   RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
     $                   TAU, TMP, TMP1


<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISEED( 4 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.210"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL                        <a name="SLAMCH.211"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="SLAMCH.212"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="LSAME.212"></a><a href="lsame.f.html#LSAME.1">LSAME</a>

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, <a name="SLARNV.216"></a><a href="slarnv.f.html#SLARNV.1">SLARNV</a>, <a name="SLARRA.216"></a><a href="slarra.f.html#SLARRA.1">SLARRA</a>, <a name="SLARRB.216"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a>, <a name="SLARRC.216"></a><a href="slarrc.f.html#SLARRC.1">SLARRC</a>, <a name="SLARRD.216"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a>,
     $                   <a name="SLASQ2.217"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>
      INFO = 0

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode RANGE
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.231"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'A'</span> ) ) THEN
         IRANGE = ALLRNG
      ELSE IF( <a name="LSAME.233"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'V'</span> ) ) THEN
         IRANGE = VALRNG
      ELSE IF( <a name="LSAME.235"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'I'</span> ) ) THEN
         IRANGE = INDRNG
      END IF

      M = 0

<span class="comment">*</span><span class="comment">     Get machine constants
</span>      SAFMIN = <a name="SLAMCH.242"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
      EPS = <a name="SLAMCH.243"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )

<span class="comment">*</span><span class="comment">     Set parameters
</span>      RTL = HNDRD*EPS
<span class="comment">*</span><span class="comment">     If one were ever to ask for less initial precision in BSRTOL,
</span><span class="comment">*</span><span class="comment">     one should keep in mind that for the subset case, the extremal
</span><span class="comment">*</span><span class="comment">     eigenvalues must be at least as accurate as the current setting
</span><span class="comment">*</span><span class="comment">     (eigenvalues in the middle need not as much accuracy)
</span>      BSRTOL = SQRT(EPS)*(0.5E-3)

<span class="comment">*</span><span class="comment">     Treat case of 1x1 matrix for quick return
</span>      IF( N.EQ.1 ) THEN
         IF( (IRANGE.EQ.ALLRNG).OR.
     $       ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
     $       ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
            M = 1
            W(1) = D(1)
<span class="comment">*</span><span class="comment">           The computation error of the eigenvalue is zero
</span>            WERR(1) = ZERO
            WGAP(1) = ZERO
            IBLOCK( 1 ) = 1
            INDEXW( 1 ) = 1
            GERS(1) = D( 1 )
            GERS(2) = D( 1 )
         ENDIF
<span class="comment">*</span><span class="comment">        store the shift for the initial RRR, which is zero in this case
</span>         E(1) = ZERO
         RETURN
      END IF

<span class="comment">*</span><span class="comment">     General case: tridiagonal matrix of order &gt; 1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
</span><span class="comment">*</span><span class="comment">     Compute maximum off-diagonal entry and pivmin.
</span>      GL = D(1)
      GU = D(1)
      EOLD = ZERO
      EMAX = ZERO
      E(N) = ZERO
      DO 5 I = 1,N
         WERR(I) = ZERO
         WGAP(I) = ZERO
         EABS = ABS( E(I) )
         IF( EABS .GE. EMAX ) THEN
            EMAX = EABS
         END IF
         TMP1 = EABS + EOLD
         GERS( 2*I-1) = D(I) - TMP1
         GL =  MIN( GL, GERS( 2*I - 1))
         GERS( 2*I ) = D(I) + TMP1
         GU = MAX( GU, GERS(2*I) )
         EOLD  = EABS
 5    CONTINUE
<span class="comment">*</span><span class="comment">     The minimum pivot allowed in the Sturm sequence for T
</span>      PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
<span class="comment">*</span><span class="comment">     Compute spectral diameter. The Gerschgorin bounds give an
</span><span class="comment">*</span><span class="comment">     estimate that is wrong by at most a factor of SQRT(2)
</span>      SPDIAM = GU - GL

<span class="comment">*</span><span class="comment">     Compute splitting points
</span>      CALL <a name="SLARRA.303"></a><a href="slarra.f.html#SLARRA.1">SLARRA</a>( N, D, E, E2, SPLTOL, SPDIAM,
     $                    NSPLIT, ISPLIT, IINFO )

<span class="comment">*</span><span class="comment">     Can force use of bisection instead of faster DQDS.
</span><span class="comment">*</span><span class="comment">     Option left in the code for future multisection work.
</span>      FORCEB = .FALSE.

      IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
<span class="comment">*</span><span class="comment">        Set interval [VL,VU] that contains all eigenvalues
</span>         VL = GL
         VU = GU
      ELSE
<span class="comment">*</span><span class="comment">        We call <a name="SLARRD.315"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a> to find crude approximations to the eigenvalues
</span><span class="comment">*</span><span class="comment">        in the desired range. In case IRANGE = INDRNG, we also obtain the
</span><span class="comment">*</span><span class="comment">        interval (VL,VU] that contains all the wanted eigenvalues.
</span><span class="comment">*</span><span class="comment">        An interval [LEFT,RIGHT] has converged if
</span><span class="comment">*</span><span class="comment">        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
</span><span class="comment">*</span><span class="comment">        <a name="SLARRD.320"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a> needs a WORK of size 4*N, IWORK of size 3*N
</span>         CALL <a name="SLARRD.321"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a>( RANGE, <span class="string">'B'</span>, N, VL, VU, IL, IU, GERS,
     $                    BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
     $                    MM, W, WERR, VL, VU, IBLOCK, INDEXW,
     $                    WORK, IWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            INFO = -1
            RETURN
         ENDIF
<span class="comment">*</span><span class="comment">        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
</span>         DO 14 I = MM+1,N
            W( I ) = ZERO
            WERR( I ) = ZERO
            IBLOCK( I ) = 0
            INDEXW( I ) = 0
 14      CONTINUE
      END IF


<span class="comment">*</span><span class="comment">**
</span><span class="comment">*</span><span class="comment">     Loop over unreduced blocks
</span>      IBEGIN = 1
      WBEGIN = 1
      DO 170 JBLK = 1, NSPLIT
         IEND = ISPLIT( JBLK )
         IN = IEND - IBEGIN + 1

<span class="comment">*</span><span class="comment">        1 X 1 block
</span>         IF( IN.EQ.1 ) THEN
            IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
     $         ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
     $        .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
     $        ) THEN
               M = M + 1
               W( M ) = D( IBEGIN )
               WERR(M) = ZERO
<span class="comment">*</span><span class="comment">              The gap for a single block doesn't matter for the later
</span><span class="comment">*</span><span class="comment">              algorithm and is assigned an arbitrary large value
</span>               WGAP(M) = ZERO
               IBLOCK( M ) = JBLK
               INDEXW( M ) = 1
               WBEGIN = WBEGIN + 1
            ENDIF
<span class="comment">*</span><span class="comment">           E( IEND ) holds the shift for the initial RRR
</span>            E( IEND ) = ZERO
            IBEGIN = IEND + 1
            GO TO 170
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Blocks of size larger than 1x1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        E( IEND ) will hold the shift for the initial RRR, for now set it =0