cstedc.f.html

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

<span class="comment">*</span><span class="comment">
</span>         CALL <a name="CSTEQR.269"></a><a href="csteqr.f.html#CSTEQR.1">CSTEQR</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="SSTEDC.273"></a><a href="sstedc.f.html#SSTEDC.1">SSTEDC</a> instead.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ICOMPZ.EQ.2 ) THEN
            CALL <a name="SLASET.276"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, RWORK, N )
            LL = N*N + 1
            CALL <a name="SSTEDC.278"></a><a href="sstedc.f.html#SSTEDC.1">SSTEDC</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="SLANST.293"></a><a href="slanst.f.html#SLANST.1">SLANST</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="SLAMCH.297"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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="SLANST.330"></a><a href="slanst.f.html#SLANST.1">SLANST</a>( <span class="string">'M'</span>, M, D( START ), E( START ) )
               CALL <a name="SLASCL.331"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
     $                      INFO )
               CALL <a name="SLASCL.333"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="CLAED0.336"></a><a href="claed0.f.html#CLAED0.1">CLAED0</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="SLASCL.346"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="SSTEQR.350"></a><a href="ssteqr.f.html#SSTEQR.1">SSTEQR</a>( <span class="string">'I'</span>, M, D( START ), E( START ), RWORK, M,
     $                      RWORK( M*M+1 ), INFO )
               CALL <a name="CLACRM.352"></a><a href="clacrm.f.html#CLACRM.1">CLACRM</a>( N, M, Z( 1, START ), LDZ, RWORK, M, WORK, N,
     $                      RWORK( M*M+1 ) )
               CALL <a name="CLACPY.354"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</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 CSWAP( 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="CSTEDC.401"></a><a href="cstedc.f.html#CSTEDC.1">CSTEDC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>