dstemr.f.html

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

      RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = LWMIN
         IWORK( 1 ) = LIWMIN
<span class="comment">*</span><span class="comment">
</span>         IF( WANTZ .AND. ALLEIG ) THEN
            NZCMIN = N
         ELSE IF( WANTZ .AND. VALEIG ) THEN
            CALL <a name="DLARRC.327"></a><a href="dlarrc.f.html#DLARRC.1">DLARRC</a>( <span class="string">'T'</span>, N, VL, VU, D, E, SAFMIN,
     $                            NZCMIN, ITMP, ITMP2, INFO )
         ELSE IF( WANTZ .AND. INDEIG ) THEN
            NZCMIN = IIU-IIL+1
         ELSE
<span class="comment">*</span><span class="comment">           WANTZ .EQ. FALSE.
</span>            NZCMIN = 0
         ENDIF
         IF( ZQUERY .AND. INFO.EQ.0 ) THEN
            Z( 1,1 ) = NZCMIN
         ELSE IF( NZC.LT.NZCMIN .AND. .NOT.ZQUERY ) THEN
            INFO = -14
         END IF
      END IF

      IF( INFO.NE.0 ) THEN
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="XERBLA.344"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DSTEMR.344"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>'</span>, -INFO )
<span class="comment">*</span><span class="comment">
</span>         RETURN
      ELSE IF( LQUERY .OR. ZQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Handle N = 0, 1, and 2 cases immediately
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         IF( ALLEIG .OR. INDEIG ) THEN
            M = 1
            W( 1 ) = D( 1 )
         ELSE
            IF( WL.LT.D( 1 ) .AND. WU.GE.D( 1 ) ) THEN
               M = 1
               W( 1 ) = D( 1 )
            END IF
         END IF
         IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
            Z( 1, 1 ) = ONE
            ISUPPZ(1) = 1
            ISUPPZ(2) = 1
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.2 ) THEN
         IF( .NOT.WANTZ ) THEN
            CALL <a name="DLAE2.377"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>( D(1), E(1), D(2), R1, R2 )
         ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
            CALL <a name="DLAEV2.379"></a><a href="dlaev2.f.html#DLAEV2.1">DLAEV2</a>( D(1), E(1), D(2), R1, R2, CS, SN )
         END IF
         IF( ALLEIG.OR.
     $      (VALEIG.AND.(R2.GT.WL).AND.
     $                  (R2.LE.WU)).OR.
     $      (INDEIG.AND.(IIL.EQ.1)) ) THEN
            M = M+1
            W( M ) = R2
            IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
               Z( 1, M ) = -SN
               Z( 2, M ) = CS
<span class="comment">*</span><span class="comment">              Note: At most one of SN and CS can be zero.
</span>               IF (SN.NE.ZERO) THEN
                  IF (CS.NE.ZERO) THEN
                     ISUPPZ(2*M-1) = 1
                     ISUPPZ(2*M-1) = 2
                  ELSE
                     ISUPPZ(2*M-1) = 1
                     ISUPPZ(2*M-1) = 1
                  END IF
               ELSE
                  ISUPPZ(2*M-1) = 2
                  ISUPPZ(2*M) = 2
               END IF
            ENDIF
         ENDIF
         IF( ALLEIG.OR.
     $      (VALEIG.AND.(R1.GT.WL).AND.
     $                  (R1.LE.WU)).OR.
     $      (INDEIG.AND.(IIU.EQ.2)) ) THEN
            M = M+1
            W( M ) = R1
            IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
               Z( 1, M ) = CS
               Z( 2, M ) = SN
<span class="comment">*</span><span class="comment">              Note: At most one of SN and CS can be zero.
</span>               IF (SN.NE.ZERO) THEN
                  IF (CS.NE.ZERO) THEN
                     ISUPPZ(2*M-1) = 1
                     ISUPPZ(2*M-1) = 2
                  ELSE
                     ISUPPZ(2*M-1) = 1
                     ISUPPZ(2*M-1) = 1
                  END IF
               ELSE
                  ISUPPZ(2*M-1) = 2
                  ISUPPZ(2*M) = 2
               END IF
            ENDIF
         ENDIF
         RETURN
      END IF

<span class="comment">*</span><span class="comment">     Continue with general N
</span>
      INDGRS = 1
      INDERR = 2*N + 1
      INDGP = 3*N + 1
      INDD = 4*N + 1
      INDE2 = 5*N + 1
      INDWRK = 6*N + 1
<span class="comment">*</span><span class="comment">
</span>      IINSPL = 1
      IINDBL = N + 1
      IINDW = 2*N + 1
      IINDWK = 3*N + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale matrix to allowable range, if necessary.
</span><span class="comment">*</span><span class="comment">     The allowable range is related to the PIVMIN parameter; see the
</span><span class="comment">*</span><span class="comment">     comments in <a name="DLARRD.448"></a><a href="dlarrd.f.html#DLARRD.1">DLARRD</a>.  The preference for scaling small values
</span><span class="comment">*</span><span class="comment">     up is heuristic; we expect users' matrices not to be close to the
</span><span class="comment">*</span><span class="comment">     RMAX threshold.
</span><span class="comment">*</span><span class="comment">
</span>      SCALE = ONE
      TNRM = <a name="DLANST.453"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'M'</span>, N, D, E )
      IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
         SCALE = RMIN / TNRM
      ELSE IF( TNRM.GT.RMAX ) THEN
         SCALE = RMAX / TNRM
      END IF
      IF( SCALE.NE.ONE ) THEN
         CALL DSCAL( N, SCALE, D, 1 )
         CALL DSCAL( N-1, SCALE, E, 1 )
         TNRM = TNRM*SCALE
         IF( VALEIG ) THEN
<span class="comment">*</span><span class="comment">           If eigenvalues in interval have to be found,
</span><span class="comment">*</span><span class="comment">           scale (WL, WU] accordingly
</span>            WL = WL*SCALE
            WU = WU*SCALE
         ENDIF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the desired eigenvalues of the tridiagonal after splitting
</span><span class="comment">*</span><span class="comment">     into smaller subblocks if the corresponding off-diagonal elements
</span><span class="comment">*</span><span class="comment">     are small
</span><span class="comment">*</span><span class="comment">     THRESH is the splitting parameter for <a name="DLARRE.474"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>
</span><span class="comment">*</span><span class="comment">     A negative THRESH forces the old splitting criterion based on the
</span><span class="comment">*</span><span class="comment">     size of the off-diagonal. A positive THRESH switches to splitting
</span><span class="comment">*</span><span class="comment">     which preserves relative accuracy.
</span><span class="comment">*</span><span class="comment">
</span>      IF( TRYRAC ) THEN
<span class="comment">*</span><span class="comment">        Test whether the matrix warrants the more expensive relative approach.
</span>         CALL <a name="DLARRR.481"></a><a href="dlarrr.f.html#DLARRR.1">DLARRR</a>( N, D, E, IINFO )
      ELSE
<span class="comment">*</span><span class="comment">        The user does not care about relative accurately eigenvalues
</span>         IINFO = -1
      ENDIF