sgeesx.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORGHR.351"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
     $                LWORK-IWRK+1, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</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">     (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment">
</span>      IWRK = ITAU
      CALL <a name="SHSEQR.361"></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.370"></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.371"></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, transform Schur vectors, and compute
</span><span class="comment">*</span><span class="comment">        reciprocal condition numbers
</span><span class="comment">*</span><span class="comment">        (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
</span><span class="comment">*</span><span class="comment">                     otherwise, need N )
</span><span class="comment">*</span><span class="comment">        (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
</span><span class="comment">*</span><span class="comment">                     otherwise, need 0 )
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="STRSEN.384"></a><a href="strsen.f.html#STRSEN.1">STRSEN</a>( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
     $                SDIM, RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
     $                IWORK, LIWORK, ICOND )
         IF( .NOT.WANTSN )
     $      MAXWRK = MAX( MAXWRK, N+2*SDIM*( N-SDIM ) )
         IF( ICOND.EQ.-15 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Not enough real workspace
</span><span class="comment">*</span><span class="comment">
</span>            INFO = -16
         ELSE IF( ICOND.EQ.-17 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Not enough integer workspace
</span><span class="comment">*</span><span class="comment">
</span>            INFO = -18
         ELSE IF( ICOND.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           <a name="STRSEN.401"></a><a href="strsen.f.html#STRSEN.1">STRSEN</a> failed to reorder or to restore standard Schur form
</span><span class="comment">*</span><span class="comment">
</span>            INFO = ICOND + N
         END IF
      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">        (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEBAK.412"></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.420"></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( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
            DUM( 1 ) = RCONDV
            CALL <a name="SLASCL.424"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
            RCONDV = DUM( 1 )
         END IF
         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.436"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
     $                      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 )
                     CALL SSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
                     A( I, I+1 ) = A( I+1, I )
                     A( I+1, I ) = ZERO
                  END IF
                  INXT = I + 2
               END IF
   20       CONTINUE
         END IF
         CALL <a name="SLASCL.472"></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
      IF( WANTSV .OR. WANTSB ) THEN
         IWORK( 1 ) = SDIM*(N-SDIM)
      ELSE
         IWORK( 1 ) = 1
      END IF
<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="SGEESX.525"></a><a href="sgeesx.f.html#SGEESX.1">SGEESX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>