dstebz.f.html

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

</span><span class="comment">*</span><span class="comment">     Decode RANGE
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.210"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'A'</span> ) ) THEN
         IRANGE = 1
      ELSE IF( <a name="LSAME.212"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'V'</span> ) ) THEN
         IRANGE = 2
      ELSE IF( <a name="LSAME.214"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'I'</span> ) ) THEN
         IRANGE = 3
      ELSE
         IRANGE = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode ORDER
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.222"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( ORDER, <span class="string">'B'</span> ) ) THEN
         IORDER = 2
      ELSE IF( <a name="LSAME.224"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( ORDER, <span class="string">'E'</span> ) ) THEN
         IORDER = 1
      ELSE
         IORDER = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for Errors
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.LE.0 ) THEN
         INFO = -1
      ELSE IF( IORDER.LE.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( IRANGE.EQ.2 ) THEN
         IF( VL.GE.VU )
     $      INFO = -5
      ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
     $          THEN
         INFO = -6
      ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
     $          THEN
         INFO = -7
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.250"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DSTEBZ.250"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize error flags
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      NCNVRG = .FALSE.
      TOOFEW = .FALSE.
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Simplifications:
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
     $   IRANGE = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">     NB is the minimum vector length for vector bisection, or 0
</span><span class="comment">*</span><span class="comment">     if only scalar is to be done.
</span><span class="comment">*</span><span class="comment">
</span>      SAFEMN = <a name="DLAMCH.275"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      ULP = <a name="DLAMCH.276"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      RTOLI = ULP*RELFAC
      NB = <a name="ILAENV.278"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DSTEBZ.278"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>'</span>, <span class="string">' '</span>, N, -1, -1, -1 )
      IF( NB.LE.1 )
     $   NB = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Special Case when N=1
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         NSPLIT = 1
         ISPLIT( 1 ) = 1
         IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
            M = 0
         ELSE
            W( 1 ) = D( 1 )
            IBLOCK( 1 ) = 1
            M = 1
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute Splitting Points
</span><span class="comment">*</span><span class="comment">
</span>      NSPLIT = 1
      WORK( N ) = ZERO
      PIVMIN = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">DIR$ NOVECTOR
</span>      DO 10 J = 2, N
         TMP1 = E( J-1 )**2
         IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
            ISPLIT( NSPLIT ) = J - 1
            NSPLIT = NSPLIT + 1
            WORK( J-1 ) = ZERO
         ELSE
            WORK( J-1 ) = TMP1
            PIVMIN = MAX( PIVMIN, TMP1 )
         END IF
   10 CONTINUE
      ISPLIT( NSPLIT ) = N
      PIVMIN = PIVMIN*SAFEMN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute Interval and ATOLI
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.EQ.3 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        RANGE='I': Compute the interval containing eigenvalues
</span><span class="comment">*</span><span class="comment">                   IL through IU.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute Gershgorin interval for entire (split) matrix
</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( 1 )
         GL = D( 1 )
         TMP1 = ZERO
<span class="comment">*</span><span class="comment">
</span>         DO 20 J = 1, N - 1
            TMP2 = SQRT( WORK( J ) )
            GU = MAX( GU, D( J )+TMP1+TMP2 )
            GL = MIN( GL, D( J )-TMP1-TMP2 )
            TMP1 = TMP2
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         GU = MAX( GU, D( N )+TMP1 )
         GL = MIN( GL, D( N )-TMP1 )
         TNORM = MAX( ABS( GL ), ABS( GU ) )
         GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
         GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute Iteration parameters
</span><span class="comment">*</span><span class="comment">
</span>         ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
     $           LOG( TWO ) ) + 2
         IF( ABSTOL.LE.ZERO ) THEN
            ATOLI = ULP*TNORM
         ELSE
            ATOLI = ABSTOL
         END IF
<span class="comment">*</span><span class="comment">
</span>         WORK( N+1 ) = GL
         WORK( N+2 ) = GL
         WORK( N+3 ) = GU
         WORK( N+4 ) = GU
         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="DLAEBZ.368"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a>( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
     $                WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
     $                IWORK, W, IBLOCK, IINFO )
<span class="comment">*</span><span class="comment">
</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">
</span>         IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
            INFO = 4
            RETURN
         END IF
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        RANGE='A' or 'V' -- Set ATOLI
</span><span class="comment">*</span><span class="comment">
</span>         TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
     $           ABS( D( N ) )+ABS( E( N-1 ) ) )
<span class="comment">*</span><span class="comment">
</span>         DO 30 J = 2, N - 1
            TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
     $              ABS( E( J ) ) )
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( ABSTOL.LE.ZERO ) THEN
            ATOLI = ULP*TNORM
         ELSE
            ATOLI = ABSTOL
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( IRANGE.EQ.2 ) THEN
            WL = VL
            WU = VU
         ELSE
            WL = ZERO
            WU = ZERO
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><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><span class="comment">*</span><span class="comment">
</span>      M = 0
      IEND = 0
      INFO = 0
      NWL = 0
      NWU = 0
<span class="comment">*</span><span class="comment">
</span>      DO 70 JB = 1, NSPLIT
         IOFF = IEND
         IBEGIN = IOFF + 1
         IEND = ISPLIT( JB )
         IN = IEND - IOFF