sgees.f.html

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

</span>      SDIM = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Perform QR iteration, accumulating Schur vectors in VS if desired
</span><span class="comment">*</span><span class="comment">     (Workspace: need N+1, prefer N+HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment">
</span>      IWRK = ITAU
      CALL <a name="SHSEQR.294"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>( <span class="string">'S'</span>, JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
     $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
      IF( IEVAL.GT.0 )
     $   INFO = IEVAL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort eigenvalues if desired
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTST .AND. INFO.EQ.0 ) THEN
         IF( SCALEA ) THEN
            CALL <a name="SLASCL.303"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
            CALL <a name="SLASCL.304"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
         END IF
         DO 10 I = 1, N
            BWORK( I ) = SELECT( WR( I ), WI( I ) )
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reorder eigenvalues and transform Schur vectors
</span><span class="comment">*</span><span class="comment">        (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="STRSEN.313"></a><a href="strsen.f.html#STRSEN.1">STRSEN</a>( <span class="string">'N'</span>, JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
     $                SDIM, S, SEP, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
     $                ICOND )
         IF( ICOND.GT.0 )
     $      INFO = N + ICOND
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( WANTVS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Undo balancing
</span><span class="comment">*</span><span class="comment">        (Workspace: need N)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEBAK.325"></a><a href="sgebak.f.html#SGEBAK.1">SGEBAK</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
     $                IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( SCALEA ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Undo scaling for the Schur form of A
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASCL.333"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'H'</span>, 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
         CALL SCOPY( N, A, LDA+1, WR, 1 )
         IF( CSCALE.EQ.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If scaling back towards underflow, adjust WI if an
</span><span class="comment">*</span><span class="comment">           offdiagonal element of a 2-by-2 block in the Schur form
</span><span class="comment">*</span><span class="comment">           underflows.
</span><span class="comment">*</span><span class="comment">
</span>            IF( IEVAL.GT.0 ) THEN
               I1 = IEVAL + 1
               I2 = IHI - 1
               CALL <a name="SLASCL.344"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, ILO-1, 1, WI,
     $                      MAX( ILO-1, 1 ), IERR )
            ELSE IF( WANTST ) THEN
               I1 = 1
               I2 = N - 1
            ELSE
               I1 = ILO
               I2 = IHI - 1
            END IF
            INXT = I1 - 1
            DO 20 I = I1, I2
               IF( I.LT.INXT )
     $            GO TO 20
               IF( WI( I ).EQ.ZERO ) THEN
                  INXT = I + 1
               ELSE
                  IF( A( I+1, I ).EQ.ZERO ) THEN
                     WI( I ) = ZERO
                     WI( I+1 ) = ZERO
                  ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
     $                     ZERO ) THEN
                     WI( I ) = ZERO
                     WI( I+1 ) = ZERO
                     IF( I.GT.1 )
     $                  CALL SSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
                     IF( N.GT.I+1 )
     $                  CALL SSWAP( N-I-1, A( I, I+2 ), LDA,
     $                              A( I+1, I+2 ), LDA )
                     IF( WANTVS ) THEN
                        CALL SSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
                     END IF
                     A( I, I+1 ) = A( I+1, I )
                     A( I+1, I ) = ZERO
                  END IF
                  INXT = I + 2
               END IF
   20       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Undo scaling for the imaginary part of the eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASCL.385"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N-IEVAL, 1,
     $                WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( WANTST .AND. INFO.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check if reordering successful
</span><span class="comment">*</span><span class="comment">
</span>         LASTSL = .TRUE.
         LST2SL = .TRUE.
         SDIM = 0
         IP = 0
         DO 30 I = 1, N
            CURSL = SELECT( WR( I ), WI( I ) )
            IF( WI( I ).EQ.ZERO ) THEN
               IF( CURSL )
     $            SDIM = SDIM + 1
               IP = 0
               IF( CURSL .AND. .NOT.LASTSL )
     $            INFO = N + 2
            ELSE
               IF( IP.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Last eigenvalue of conjugate pair
</span><span class="comment">*</span><span class="comment">
</span>                  CURSL = CURSL .OR. LASTSL
                  LASTSL = CURSL
                  IF( CURSL )
     $               SDIM = SDIM + 2
                  IP = -1
                  IF( CURSL .AND. .NOT.LST2SL )
     $               INFO = N + 2
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 First eigenvalue of conjugate pair
</span><span class="comment">*</span><span class="comment">
</span>                  IP = 1
               END IF
            END IF
            LST2SL = LASTSL
            LASTSL = CURSL
   30    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = MAXWRK
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SGEES.432"></a><a href="sgees.f.html#SGEES.1">SGEES</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>