dstedc.f.html

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

</span>         EPS = <a name="DLAMCH.270"></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>   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="DLANST.307"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'M'</span>, M, D( START ), E( START ) )
               CALL <a name="DLASCL.308"></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.310"></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>               IF( ICOMPZ.EQ.1 ) THEN
                  STRTRW = 1
               ELSE
                  STRTRW = START
               END IF
               CALL <a name="DLAED0.318"></a><a href="dlaed0.f.html#DLAED0.1">DLAED0</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="DLASCL.329"></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
               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="DSTEQR.339"></a><a href="dsteqr.f.html#DSTEQR.1">DSTEQR</a>( <span class="string">'I'</span>, M, D( START ), E( START ), WORK, M,
     $                         WORK( M*M+1 ), INFO )
                  CALL <a name="DLACPY.341"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'A'</span>, N, M, Z( 1, START ), LDZ,
     $                         WORK( STOREZ ), N )
                  CALL DGEMM( <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="DSTEQR.347"></a><a href="dsteqr.f.html#DSTEQR.1">DSTEQR</a>( <span class="string">'I'</span>, M, D( START ), E( START ),
     $                         Z( START, START ), LDZ, WORK, INFO )
               ELSE
                  CALL <a name="DSTERF.350"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</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="DLASRT.373"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</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 DSWAP( 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="DSTEDC.405"></a><a href="dstedc.f.html#DSTEDC.1">DSTEDC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>