shsein.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>            IF( FROMQR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              If affiliation of eigenvalues is known, check whether
</span><span class="comment">*</span><span class="comment">              the matrix splits.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Determine KL and KR such that 1 &lt;= KL &lt;= K &lt;= KR &lt;= N
</span><span class="comment">*</span><span class="comment">              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
</span><span class="comment">*</span><span class="comment">              KR = N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Then inverse iteration can be performed with the
</span><span class="comment">*</span><span class="comment">              submatrix H(KL:N,KL:N) for a left eigenvector, and with
</span><span class="comment">*</span><span class="comment">              the submatrix H(1:KR,1:KR) for a right eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>               DO 20 I = K, KL + 1, -1
                  IF( H( I, I-1 ).EQ.ZERO )
     $               GO TO 30
   20          CONTINUE
   30          CONTINUE
               KL = I
               IF( K.GT.KR ) THEN
                  DO 40 I = K, N - 1
                     IF( H( I+1, I ).EQ.ZERO )
     $                  GO TO 50
   40             CONTINUE
   50             CONTINUE
                  KR = I
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( KL.NE.KLN ) THEN
               KLN = KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
</span><span class="comment">*</span><span class="comment">              has not ben computed before.
</span><span class="comment">*</span><span class="comment">
</span>               HNORM = <a name="SLANHS.310"></a><a href="slanhs.f.html#SLANHS.1">SLANHS</a>( <span class="string">'I'</span>, KR-KL+1, H( KL, KL ), LDH, WORK )
               IF( HNORM.GT.ZERO ) THEN
                  EPS3 = HNORM*ULP
               ELSE
                  EPS3 = SMLNUM
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Perturb eigenvalue if it is close to any previous
</span><span class="comment">*</span><span class="comment">           selected eigenvalues affiliated to the submatrix
</span><span class="comment">*</span><span class="comment">           H(KL:KR,KL:KR). Close roots are modified by EPS3.
</span><span class="comment">*</span><span class="comment">
</span>            WKR = WR( K )
            WKI = WI( K )
   60       CONTINUE
            DO 70 I = K - 1, KL, -1
               IF( SELECT( I ) .AND. ABS( WR( I )-WKR )+
     $             ABS( WI( I )-WKI ).LT.EPS3 ) THEN
                  WKR = WKR + EPS3
                  GO TO 60
               END IF
   70       CONTINUE
            WR( K ) = WKR
<span class="comment">*</span><span class="comment">
</span>            PAIR = WKI.NE.ZERO
            IF( PAIR ) THEN
               KSI = KSR + 1
            ELSE
               KSI = KSR
            END IF
            IF( LEFTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute left eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLAEIN.344"></a><a href="slaein.f.html#SLAEIN.1">SLAEIN</a>( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
     $                      WKR, WKI, VL( KL, KSR ), VL( KL, KSI ),
     $                      WORK, LDWORK, WORK( N*N+N+1 ), EPS3, SMLNUM,
     $                      BIGNUM, IINFO )
               IF( IINFO.GT.0 ) THEN
                  IF( PAIR ) THEN
                     INFO = INFO + 2
                  ELSE
                     INFO = INFO + 1
                  END IF
                  IFAILL( KSR ) = K
                  IFAILL( KSI ) = K
               ELSE
                  IFAILL( KSR ) = 0
                  IFAILL( KSI ) = 0
               END IF
               DO 80 I = 1, KL - 1
                  VL( I, KSR ) = ZERO
   80          CONTINUE
               IF( PAIR ) THEN
                  DO 90 I = 1, KL - 1
                     VL( I, KSI ) = ZERO
   90             CONTINUE
               END IF
            END IF
            IF( RIGHTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute right eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLAEIN.373"></a><a href="slaein.f.html#SLAEIN.1">SLAEIN</a>( .TRUE., NOINIT, KR, H, LDH, WKR, WKI,
     $                      VR( 1, KSR ), VR( 1, KSI ), WORK, LDWORK,
     $                      WORK( N*N+N+1 ), EPS3, SMLNUM, BIGNUM,
     $                      IINFO )
               IF( IINFO.GT.0 ) THEN
                  IF( PAIR ) THEN
                     INFO = INFO + 2
                  ELSE
                     INFO = INFO + 1
                  END IF
                  IFAILR( KSR ) = K
                  IFAILR( KSI ) = K
               ELSE
                  IFAILR( KSR ) = 0
                  IFAILR( KSI ) = 0
               END IF
               DO 100 I = KR + 1, N
                  VR( I, KSR ) = ZERO
  100          CONTINUE
               IF( PAIR ) THEN
                  DO 110 I = KR + 1, N
                     VR( I, KSI ) = ZERO
  110             CONTINUE
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( PAIR ) THEN
               KSR = KSR + 2
            ELSE
               KSR = KSR + 1
            END IF
         END IF
  120 CONTINUE
<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="SHSEIN.409"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>