sstedc.f.html

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

<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>   10    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
   20       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 20
               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.EQ.1 ) THEN
               START = FINISH + 1
               GO TO 10
            END IF
            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.306"></a><a href="slanst.f.html#SLANST.1">SLANST</a>( <span class="string">'M'</span>, M, D( START ), E( START ) )
               CALL <a name="SLASCL.307"></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.309"></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>               IF( ICOMPZ.EQ.1 ) THEN
                  STRTRW = 1
               ELSE
                  STRTRW = START
               END IF
               CALL <a name="SLAED0.317"></a><a href="slaed0.f.html#SLAED0.1">SLAED0</a>( ICOMPZ, N, M, D( START ), E( START ),
     $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
     $                      WORK( STOREZ ), IWORK, INFO )
               IF( INFO.NE.0 ) THEN
                  INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
     $                   MOD( INFO, ( M+1 ) ) + START - 1
                  GO TO 50
               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.328"></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
               IF( ICOMPZ.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Since QR won't update a Z matrix which is larger than
</span><span class="comment">*</span><span class="comment">                 the length of D, we must solve the sub-problem in a
</span><span class="comment">*</span><span class="comment">                 workspace and then multiply back into Z.
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SSTEQR.338"></a><a href="ssteqr.f.html#SSTEQR.1">SSTEQR</a>( <span class="string">'I'</span>, M, D( START ), E( START ), WORK, M,
     $                         WORK( M*M+1 ), INFO )
                  CALL <a name="SLACPY.340"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, N, M, Z( 1, START ), LDZ,
     $                         WORK( STOREZ ), N )
                  CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, N, M, M, ONE,
     $                        WORK( STOREZ ), N, WORK, M, ZERO,
     $                        Z( 1, START ), LDZ )
               ELSE IF( ICOMPZ.EQ.2 ) THEN
                  CALL <a name="SSTEQR.346"></a><a href="ssteqr.f.html#SSTEQR.1">SSTEQR</a>( <span class="string">'I'</span>, M, D( START ), E( START ),
     $                         Z( START, START ), LDZ, WORK, INFO )
               ELSE
                  CALL <a name="SSTERF.349"></a><a href="ssterf.f.html#SSTERF.1">SSTERF</a>( M, D( START ), E( START ), INFO )
               END IF
               IF( INFO.NE.0 ) THEN
                  INFO = START*( N+1 ) + FINISH
                  GO TO 50
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            START = FINISH + 1
            GO TO 10
         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
            IF( ICOMPZ.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Use Quick Sort
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLASRT.372"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'I'</span>, N, D, 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">              Use Selection Sort to minimize swaps of eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>               DO 40 II = 2, N
                  I = II - 1
                  K = I
                  P = D( I )
                  DO 30 J = II, N
                     IF( D( J ).LT.P ) THEN
                        K = J
                        P = D( J )
                     END IF
   30             CONTINUE
                  IF( K.NE.I ) THEN
                     D( K ) = D( I )
                     D( I ) = P
                     CALL SSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
                  END IF
   40          CONTINUE
            END IF
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>   50 CONTINUE
      WORK( 1 ) = LWMIN
      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="SSTEDC.404"></a><a href="sstedc.f.html#SSTEDC.1">SSTEDC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>