dstebz.f.html

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

<span class="comment">*</span><span class="comment">
</span>         IF( IN.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Special Case -- IN=1
</span><span class="comment">*</span><span class="comment">
</span>            IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
     $         NWL = NWL + 1
            IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
     $         NWU = NWU + 1
            IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
     $          D( IBEGIN )-PIVMIN ) ) THEN
               M = M + 1
               W( M ) = D( IBEGIN )
               IBLOCK( M ) = JB
            END IF
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           General Case -- IN &gt; 1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute Gershgorin Interval
</span><span class="comment">*</span><span class="comment">           and use it as the initial interval
</span><span class="comment">*</span><span class="comment">
</span>            GU = D( IBEGIN )
            GL = D( IBEGIN )
            TMP1 = ZERO
<span class="comment">*</span><span class="comment">
</span>            DO 40 J = IBEGIN, IEND - 1
               TMP2 = ABS( E( J ) )
               GU = MAX( GU, D( J )+TMP1+TMP2 )
               GL = MIN( GL, D( J )-TMP1-TMP2 )
               TMP1 = TMP2
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            GU = MAX( GU, D( IEND )+TMP1 )
            GL = MIN( GL, D( IEND )-TMP1 )
            BNORM = MAX( ABS( GL ), ABS( GU ) )
            GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
            GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute ATOLI for the current submatrix
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABSTOL.LE.ZERO ) THEN
               ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
            ELSE
               ATOLI = ABSTOL
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( IRANGE.GT.1 ) THEN
               IF( GU.LT.WL ) THEN
                  NWL = NWL + IN
                  NWU = NWU + IN
                  GO TO 70
               END IF
               GL = MAX( GL, WL )
               GU = MIN( GU, WU )
               IF( GL.GE.GU )
     $            GO TO 70
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Set Up Initial Interval
</span><span class="comment">*</span><span class="comment">
</span>            WORK( N+1 ) = GL
            WORK( N+IN+1 ) = GU
            CALL <a name="DLAEBZ.497"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a>( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
     $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
     $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
     $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
<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">
</span><span class="comment">*</span><span class="comment">           Compute Eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>            ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
     $              LOG( TWO ) ) + 2
            CALL <a name="DLAEBZ.510"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a>( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
     $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
     $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
     $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy Eigenvalues Into W and IBLOCK
</span><span class="comment">*</span><span class="comment">           Use -JB for block number for unconverged eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>            DO 60 J = 1, IOUT
               TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Flag non-convergence.
</span><span class="comment">*</span><span class="comment">
</span>               IF( J.GT.IOUT-IINFO ) THEN
                  NCNVRG = .TRUE.
                  IB = -JB
               ELSE
                  IB = JB
               END IF
               DO 50 JE = IWORK( J ) + 1 + IWOFF,
     $                 IWORK( J+IN ) + IWOFF
                  W( JE ) = TMP1
                  IBLOCK( JE ) = IB
   50          CONTINUE
   60       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            M = M + IM
         END IF
   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
</span><span class="comment">*</span><span class="comment">     If NWL+1 &lt; IL or NWU &gt; IU, discard extra eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.EQ.3 ) THEN
         IM = 0
         IDISCL = IL - 1 - NWL
         IDISCU = NWU - IU
<span class="comment">*</span><span class="comment">
</span>         IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
            DO 80 JE = 1, M
               IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
                  IDISCL = IDISCL - 1
               ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
                  IDISCU = IDISCU - 1
               ELSE
                  IM = IM + 1
                  W( IM ) = W( JE )
                  IBLOCK( IM ) = IBLOCK( JE )
               END IF
   80       CONTINUE
            M = IM
         END IF
         IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Code to deal with effects of bad arithmetic:
</span><span class="comment">*</span><span class="comment">           Some low eigenvalues to be discarded are not in (WL,WLU],
</span><span class="comment">*</span><span class="comment">           or high eigenvalues to be discarded are not in (WUL,WU]
</span><span class="comment">*</span><span class="comment">           so just kill off the smallest IDISCL/largest IDISCU
</span><span class="comment">*</span><span class="comment">           eigenvalues, by simply finding the smallest/largest
</span><span class="comment">*</span><span class="comment">           eigenvalue(s).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           (If N(w) is monotone non-decreasing, this should never
</span><span class="comment">*</span><span class="comment">               happen.)
</span><span class="comment">*</span><span class="comment">
</span>            IF( IDISCL.GT.0 ) THEN
               WKILL = WU
               DO 100 JDISC = 1, IDISCL
                  IW = 0
                  DO 90 JE = 1, M
                     IF( IBLOCK( JE ).NE.0 .AND.
     $                   ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
                        IW = JE
                        WKILL = W( JE )
                     END IF
   90             CONTINUE
                  IBLOCK( IW ) = 0
  100          CONTINUE
            END IF
            IF( IDISCU.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span>               WKILL = WL
               DO 120 JDISC = 1, IDISCU
                  IW = 0
                  DO 110 JE = 1, M
                     IF( IBLOCK( JE ).NE.0 .AND.
     $                   ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
                        IW = JE
                        WKILL = W( JE )
                     END IF
  110             CONTINUE
                  IBLOCK( IW ) = 0
  120          CONTINUE
            END IF
            IM = 0
            DO 130 JE = 1, M
               IF( IBLOCK( JE ).NE.0 ) THEN
                  IM = IM + 1
                  W( IM ) = W( JE )
                  IBLOCK( IM ) = IBLOCK( JE )
               END IF
  130       CONTINUE
            M = IM
         END IF
         IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
            TOOFEW = .TRUE.
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If ORDER='B', do nothing -- the eigenvalues are already sorted
</span><span class="comment">*</span><span class="comment">        by block.
</span><span class="comment">*</span><span class="comment">     If ORDER='E', sort the eigenvalues from smallest to largest
</span><span class="comment">*</span><span class="comment">
</span>      IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
         DO 150 JE = 1, M - 1
            IE = 0
            TMP1 = W( JE )
            DO 140 J = JE + 1, M
               IF( W( J ).LT.TMP1 ) THEN
                  IE = J
                  TMP1 = W( J )
               END IF
  140       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( IE.NE.0 ) THEN
               ITMP1 = IBLOCK( IE )
               W( IE ) = W( JE )
               IBLOCK( IE ) = IBLOCK( JE )
               W( JE ) = TMP1
               IBLOCK( JE ) = ITMP1
            END IF
  150    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( NCNVRG )
     $   INFO = INFO + 1
      IF( TOOFEW )
     $   INFO = INFO + 2
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DSTEBZ.650"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>