slarrd.f.html

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

         WORK( N+5 ) = GL
         WORK( N+6 ) = GU
         IWORK( 1 ) = -1
         IWORK( 2 ) = -1
         IWORK( 3 ) = N + 1
         IWORK( 4 ) = N + 1
         IWORK( 5 ) = IL - 1
         IWORK( 6 ) = IU
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLAEBZ.361"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN,
     $         D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
     $                IWORK, W, IBLOCK, IINFO )
         IF( IINFO .NE. 0 ) THEN
            INFO = IINFO
            RETURN
         END IF
<span class="comment">*</span><span class="comment">        On exit, output intervals may not be ordered by ascending negcount
</span>         IF( IWORK( 6 ).EQ.IU ) THEN
            WL = WORK( N+1 )
            WLU = WORK( N+3 )
            NWL = IWORK( 1 )
            WU = WORK( N+4 )
            WUL = WORK( N+2 )
            NWU = IWORK( 4 )
         ELSE
            WL = WORK( N+2 )
            WLU = WORK( N+4 )
            NWL = IWORK( 2 )
            WU = WORK( N+3 )
            WUL = WORK( N+1 )
            NWU = IWORK( 3 )
         END IF
<span class="comment">*</span><span class="comment">        On exit, the interval [WL, WLU] contains a value with negcount NWL,
</span><span class="comment">*</span><span class="comment">        and [WUL, WU] contains a value with negcount NWU.
</span>         IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
            INFO = 4
            RETURN
         END IF

      ELSEIF( IRANGE.EQ.VALRNG ) THEN
         WL = VL
         WU = VU

      ELSEIF( IRANGE.EQ.ALLRNG ) THEN
         WL = GL
         WU = GU
      ENDIF



<span class="comment">*</span><span class="comment">     Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
</span><span class="comment">*</span><span class="comment">     NWL accumulates the number of eigenvalues .le. WL,
</span><span class="comment">*</span><span class="comment">     NWU accumulates the number of eigenvalues .le. WU
</span>      M = 0
      IEND = 0
      INFO = 0
      NWL = 0
      NWU = 0
<span class="comment">*</span><span class="comment">
</span>      DO 70 JBLK = 1, NSPLIT
         IOFF = IEND
         IBEGIN = IOFF + 1
         IEND = ISPLIT( JBLK )
         IN = IEND - IOFF
<span class="comment">*</span><span class="comment">
</span>         IF( IN.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">           1x1 block
</span>            IF( WL.GE.D( IBEGIN )-PIVMIN )
     $         NWL = NWL + 1
            IF( WU.GE.D( IBEGIN )-PIVMIN )
     $         NWU = NWU + 1
            IF( IRANGE.EQ.ALLRNG .OR.
     $           ( WL.LT.D( IBEGIN )-PIVMIN
     $             .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) 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>               IBLOCK( M ) = JBLK
               INDEXW( M ) = 1
            END IF

<span class="comment">*</span><span class="comment">        Disabled 2x2 case because of a failure on the following matrix
</span><span class="comment">*</span><span class="comment">        RANGE = 'I', IL = IU = 4
</span><span class="comment">*</span><span class="comment">          Original Tridiagonal, d = [
</span><span class="comment">*</span><span class="comment">           -0.150102010615740E+00
</span><span class="comment">*</span><span class="comment">           -0.849897989384260E+00
</span><span class="comment">*</span><span class="comment">           -0.128208148052635E-15
</span><span class="comment">*</span><span class="comment">            0.128257718286320E-15
</span><span class="comment">*</span><span class="comment">          ];
</span><span class="comment">*</span><span class="comment">          e = [
</span><span class="comment">*</span><span class="comment">           -0.357171383266986E+00
</span><span class="comment">*</span><span class="comment">           -0.180411241501588E-15
</span><span class="comment">*</span><span class="comment">           -0.175152352710251E-15
</span><span class="comment">*</span><span class="comment">          ];
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         ELSE IF( IN.EQ.2 ) THEN
</span><span class="comment">*</span><span class="comment">*           2x2 block
</span><span class="comment">*</span><span class="comment">            DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
</span><span class="comment">*</span><span class="comment">            TMP1 = HALF*(D(IBEGIN)+D(IEND))
</span><span class="comment">*</span><span class="comment">            L1 = TMP1 - DISC
</span><span class="comment">*</span><span class="comment">            IF( WL.GE. L1-PIVMIN )
</span><span class="comment">*</span><span class="comment">     $         NWL = NWL + 1
</span><span class="comment">*</span><span class="comment">            IF( WU.GE. L1-PIVMIN )
</span><span class="comment">*</span><span class="comment">     $         NWU = NWU + 1
</span><span class="comment">*</span><span class="comment">            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
</span><span class="comment">*</span><span class="comment">     $          L1-PIVMIN ) ) THEN
</span><span class="comment">*</span><span class="comment">               M = M + 1
</span><span class="comment">*</span><span class="comment">               W( M ) = L1
</span><span class="comment">*</span><span class="comment">*              The uncertainty of eigenvalues of a 2x2 matrix is very small
</span><span class="comment">*</span><span class="comment">               WERR( M ) = EPS * ABS( W( M ) ) * TWO
</span><span class="comment">*</span><span class="comment">               IBLOCK( M ) = JBLK
</span><span class="comment">*</span><span class="comment">               INDEXW( M ) = 1
</span><span class="comment">*</span><span class="comment">            ENDIF
</span><span class="comment">*</span><span class="comment">            L2 = TMP1 + DISC
</span><span class="comment">*</span><span class="comment">            IF( WL.GE. L2-PIVMIN )
</span><span class="comment">*</span><span class="comment">     $         NWL = NWL + 1
</span><span class="comment">*</span><span class="comment">            IF( WU.GE. L2-PIVMIN )
</span><span class="comment">*</span><span class="comment">     $         NWU = NWU + 1
</span><span class="comment">*</span><span class="comment">            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
</span><span class="comment">*</span><span class="comment">     $          L2-PIVMIN ) ) THEN
</span><span class="comment">*</span><span class="comment">               M = M + 1
</span><span class="comment">*</span><span class="comment">               W( M ) = L2
</span><span class="comment">*</span><span class="comment">*              The uncertainty of eigenvalues of a 2x2 matrix is very small
</span><span class="comment">*</span><span class="comment">               WERR( M ) = EPS * ABS( W( M ) ) * TWO
</span><span class="comment">*</span><span class="comment">               IBLOCK( M ) = JBLK
</span><span class="comment">*</span><span class="comment">               INDEXW( M ) = 2
</span><span class="comment">*</span><span class="comment">            ENDIF
</span>         ELSE
<span class="comment">*</span><span class="comment">           General Case - block of size IN &gt;= 2
</span><span class="comment">*</span><span class="comment">           Compute local Gerschgorin interval and use it as the initial
</span><span class="comment">*</span><span class="comment">           interval for <a name="SLAEBZ.484"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>
</span>            GU = D( IBEGIN )
            GL = D( IBEGIN )
            TMP1 = ZERO

            DO 40 J = IBEGIN, IEND
               GL =  MIN( GL, GERS( 2*J - 1))
               GU = MAX( GU, GERS(2*J) )
   40       CONTINUE
            SPDIAM = GU - GL
            GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
            GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
<span class="comment">*</span><span class="comment">
</span>            IF( IRANGE.GT.1 ) THEN
               IF( GU.LT.WL ) THEN
<span class="comment">*</span><span class="comment">                 the local block contains none of the wanted eigenvalues
</span>                  NWL = NWL + IN
                  NWU = NWU + IN
                  GO TO 70
               END IF
<span class="comment">*</span><span class="comment">              refine search interval if possible, only range (WL,WU] matters
</span>               GL = MAX( GL, WL )
               GU = MIN( GU, WU )
               IF( GL.GE.GU )
     $            GO TO 70
            END IF

<span class="comment">*</span><span class="comment">           Find negcount of initial interval boundaries GL and GU
</span>            WORK( N+1 ) = GL
            WORK( N+IN+1 ) = GU
            CALL <a name="SLAEBZ.514"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
     $                   D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
     $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
     $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
            IF( IINFO .NE. 0 ) THEN
               INFO = IINFO
               RETURN
            END IF
<span class="comment">*</span><span class="comment">
</span>            NWL = NWL + IWORK( 1 )
            NWU = NWU + IWORK( IN+1 )
            IWOFF = M - IWORK( 1 )

<span class="comment">*</span><span class="comment">           Compute Eigenvalues
</span>            ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
     $              LOG( TWO ) ) + 2
            CALL <a name="SLAEBZ.530"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
     $                   D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
     $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
     $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
            IF( IINFO .NE. 0 ) THEN
               INFO = IINFO
               RETURN