zstedc.f.html

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

</span>      IF( N.LE.SMLSIZ ) THEN
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZSTEQR.270"></a><a href="zsteqr.f.html#ZSTEQR.1">ZSTEQR</a>( COMPZ, N, D, E, Z, LDZ, RWORK, INFO )
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If COMPZ = 'I', we simply call <a name="DSTEDC.274"></a><a href="dstedc.f.html#DSTEDC.1">DSTEDC</a> instead.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ICOMPZ.EQ.2 ) THEN
            CALL <a name="DLASET.277"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, RWORK, N )
            LL = N*N + 1
            CALL <a name="DSTEDC.279"></a><a href="dstedc.f.html#DSTEDC.1">DSTEDC</a>( <span class="string">'I'</span>, N, D, E, RWORK, N,
     $                   RWORK( LL ), LRWORK-LL+1, IWORK, LIWORK, INFO )
            DO 20 J = 1, N
               DO 10 I = 1, N
                  Z( I, J ) = RWORK( ( J-1 )*N+I )
   10          CONTINUE
   20       CONTINUE
            GO TO 70
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        From now on, only option left to be handled is COMPZ = 'V',
</span><span class="comment">*</span><span class="comment">        i.e. ICOMPZ = 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale.
</span><span class="comment">*</span><span class="comment">
</span>         ORGNRM = <a name="DLANST.294"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'M'</span>, N, D, E )
         IF( ORGNRM.EQ.ZERO )
     $      GO TO 70
<span class="comment">*</span><span class="comment">
</span>         EPS = <a name="DLAMCH.298"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
<span class="comment">*</span><span class="comment">
</span>         START = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        while ( START &lt;= N )
</span><span class="comment">*</span><span class="comment">
</span>   30    CONTINUE
         IF( START.LE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Let FINISH be the position of the next subdiagonal entry
</span><span class="comment">*</span><span class="comment">           such that E( FINISH ) &lt;= TINY or FINISH = N if no such
</span><span class="comment">*</span><span class="comment">           subdiagonal exists.  The matrix identified by the elements
</span><span class="comment">*</span><span class="comment">           between START and FINISH constitutes an independent
</span><span class="comment">*</span><span class="comment">           sub-problem.
</span><span class="comment">*</span><span class="comment">
</span>            FINISH = START
   40       CONTINUE
            IF( FINISH.LT.N ) THEN
               TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
     $                    SQRT( ABS( D( FINISH+1 ) ) )
               IF( ABS( E( FINISH ) ).GT.TINY ) THEN
                  FINISH = FINISH + 1
                  GO TO 40
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           (Sub) Problem determined.  Compute its size and solve it.
</span><span class="comment">*</span><span class="comment">
</span>            M = FINISH - START + 1
            IF( M.GT.SMLSIZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Scale.
</span><span class="comment">*</span><span class="comment">
</span>               ORGNRM = <a name="DLANST.331"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'M'</span>, M, D( START ), E( START ) )
               CALL <a name="DLASCL.332"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
     $                      INFO )
               CALL <a name="DLASCL.334"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
     $                      M-1, INFO )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLAED0.337"></a><a href="zlaed0.f.html#ZLAED0.1">ZLAED0</a>( N, M, D( START ), E( START ), Z( 1, START ),
     $                      LDZ, WORK, N, RWORK, IWORK, INFO )
               IF( INFO.GT.0 ) THEN
                  INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
     $                   MOD( INFO, ( M+1 ) ) + START - 1
                  GO TO 70
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Scale back.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLASCL.347"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
     $                      INFO )
<span class="comment">*</span><span class="comment">
</span>            ELSE
               CALL <a name="DSTEQR.351"></a><a href="dsteqr.f.html#DSTEQR.1">DSTEQR</a>( <span class="string">'I'</span>, M, D( START ), E( START ), RWORK, M,
     $                      RWORK( M*M+1 ), INFO )
               CALL <a name="ZLACRM.353"></a><a href="zlacrm.f.html#ZLACRM.1">ZLACRM</a>( N, M, Z( 1, START ), LDZ, RWORK, M, WORK, N,
     $                      RWORK( M*M+1 ) )
               CALL <a name="ZLACPY.355"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'A'</span>, N, M, WORK, N, Z( 1, START ), LDZ )
               IF( INFO.GT.0 ) THEN
                  INFO = START*( N+1 ) + FINISH
                  GO TO 70
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            START = FINISH + 1
            GO TO 30
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        endwhile
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If the problem split any number of times, then the eigenvalues
</span><span class="comment">*</span><span class="comment">        will not be properly ordered.  Here we permute the eigenvalues
</span><span class="comment">*</span><span class="comment">        (and the associated eigenvectors) into ascending order.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.NE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Use Selection Sort to minimize swaps of eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>            DO 60 II = 2, N
               I = II - 1
               K = I
               P = D( I )
               DO 50 J = II, N
                  IF( D( J ).LT.P ) THEN
                     K = J
                     P = D( J )
                  END IF
   50          CONTINUE
               IF( K.NE.I ) THEN
                  D( K ) = D( I )
                  D( I ) = P
                  CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
               END IF
   60       CONTINUE
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>   70 CONTINUE
      WORK( 1 ) = LWMIN
      RWORK( 1 ) = LRWMIN
      IWORK( 1 ) = LIWMIN
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZSTEDC.402"></a><a href="zstedc.f.html#ZSTEDC.1">ZSTEDC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>