dstevr.f.html

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

     $      Z( 1, 1 ) = ONE
         RETURN
      END IF
<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">
</span>      SAFMIN = <a name="DLAMCH.314"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      EPS = <a name="DLAMCH.315"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = SAFMIN / EPS
      BIGNUM = ONE / SMLNUM
      RMIN = SQRT( SMLNUM )
      RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
<span class="comment">*</span><span class="comment">
</span><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">
</span>      ISCALE = 0
      VLL = VL
      VUU = VU
<span class="comment">*</span><span class="comment">
</span>      TNRM = <a name="DLANST.328"></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
         ISCALE = 1
         SIGMA = RMIN / TNRM
      ELSE IF( TNRM.GT.RMAX ) THEN
         ISCALE = 1
         SIGMA = RMAX / TNRM
      END IF
      IF( ISCALE.EQ.1 ) THEN
         CALL DSCAL( N, SIGMA, D, 1 )
         CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
         IF( VALEIG ) THEN
            VLL = VL*SIGMA
            VUU = VU*SIGMA
         END IF
      END IF

<span class="comment">*</span><span class="comment">     Initialize indices into workspaces.  Note: These indices are used only
</span><span class="comment">*</span><span class="comment">     if <a name="DSTERF.346"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> or <a name="DSTEMR.346"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> fail.
</span>
<span class="comment">*</span><span class="comment">     IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in <a name="DSTEBZ.348"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and
</span><span class="comment">*</span><span class="comment">     stores the block indices of each of the M&lt;=N eigenvalues.
</span>      INDIBL = 1
<span class="comment">*</span><span class="comment">     IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in <a name="DSTEBZ.351"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and
</span><span class="comment">*</span><span class="comment">     stores the starting and finishing indices of each block.
</span>      INDISP = INDIBL + N
<span class="comment">*</span><span class="comment">     IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
</span><span class="comment">*</span><span class="comment">     that corresponding to eigenvectors that fail to converge in
</span><span class="comment">*</span><span class="comment">     <a name="DSTEIN.356"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.  This information is discarded; if any fail, the driver
</span><span class="comment">*</span><span class="comment">     returns INFO &gt; 0.
</span>      INDIFL = INDISP + N
<span class="comment">*</span><span class="comment">     INDIWO is the offset of the remaining integer workspace.
</span>      INDIWO = INDISP + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If all eigenvalues are desired, then
</span><span class="comment">*</span><span class="comment">     call <a name="DSTERF.363"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> or <a name="DSTEMR.363"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>.  If this fails for some eigenvalue, then
</span><span class="comment">*</span><span class="comment">     try <a name="DSTEBZ.364"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>      TEST = .FALSE.
      IF( INDEIG ) THEN
         IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
            TEST = .TRUE.
         END IF
      END IF
      IF( ( ALLEIG .OR. TEST ) .AND. IEEEOK.EQ.1 ) THEN
         CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
         IF( .NOT.WANTZ ) THEN
            CALL DCOPY( N, D, 1, W, 1 )
            CALL <a name="DSTERF.377"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>( N, W, WORK, INFO )
         ELSE
            CALL DCOPY( N, D, 1, WORK( N+1 ), 1 )
            IF (ABSTOL .LE. TWO*N*EPS) THEN
               TRYRAC = .TRUE.
            ELSE
               TRYRAC = .FALSE.
            END IF
            CALL <a name="DSTEMR.385"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>( JOBZ, <span class="string">'A'</span>, N, WORK( N+1 ), WORK, VL, VU, IL,
     $                   IU, M, W, Z, LDZ, N, ISUPPZ, TRYRAC,
     $                   WORK( 2*N+1 ), LWORK-2*N, IWORK, LIWORK, INFO )
<span class="comment">*</span><span class="comment">
</span>         END IF
         IF( INFO.EQ.0 ) THEN
            M = N
            GO TO 10
         END IF
         INFO = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Otherwise, call <a name="DSTEBZ.397"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and, if eigenvectors are desired, <a name="DSTEIN.397"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         ORDER = <span class="string">'B'</span>
      ELSE
         ORDER = <span class="string">'E'</span>
      END IF

      CALL <a name="DSTEBZ.405"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
     $             NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK,
     $             IWORK( INDIWO ), INFO )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         CALL <a name="DSTEIN.410"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
     $                Z, LDZ, WORK, IWORK( INDIWO ), IWORK( INDIFL ),
     $                INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If matrix was scaled, then rescale eigenvalues appropriately.
</span><span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
      IF( ISCALE.EQ.1 ) THEN
         IF( INFO.EQ.0 ) THEN
            IMAX = M
         ELSE
            IMAX = INFO - 1
         END IF
         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If eigenvalues are not in order, then sort them, along with
</span><span class="comment">*</span><span class="comment">     eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         DO 30 J = 1, M - 1
            I = 0
            TMP1 = W( J )
            DO 20 JJ = J + 1, M
               IF( W( JJ ).LT.TMP1 ) THEN
                  I = JJ
                  TMP1 = W( JJ )
               END IF
   20       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( I.NE.0 ) THEN
               ITMP1 = IWORK( I )
               W( I ) = W( J )
               IWORK( I ) = IWORK( J )
               W( J ) = TMP1
               IWORK( J ) = ITMP1
               CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
            END IF
   30    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      Causes problems with tests 19 &amp; 20:
</span><span class="comment">*</span><span class="comment">      IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
      IWORK( 1 ) = LIWMIN
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DSTEVR.460"></a><a href="dstevr.f.html#DSTEVR.1">DSTEVR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>