chsein.f.html

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

</span>      EXTERNAL           <a name="CLAEIN.170"></a><a href="claein.f.html#CLAEIN.1">CLAEIN</a>, <a name="XERBLA.170"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, MAX, REAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      REAL               CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode and test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      BOTHV = <a name="LSAME.185"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'B'</span> )
      RIGHTV = <a name="LSAME.186"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'R'</span> ) .OR. BOTHV
      LEFTV = <a name="LSAME.187"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'L'</span> ) .OR. BOTHV
<span class="comment">*</span><span class="comment">
</span>      FROMQR = <a name="LSAME.189"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EIGSRC, <span class="string">'Q'</span> )
<span class="comment">*</span><span class="comment">
</span>      NOINIT = <a name="LSAME.191"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( INITV, <span class="string">'N'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set M to the number of columns required to store the selected
</span><span class="comment">*</span><span class="comment">     eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      DO 10 K = 1, N
         IF( SELECT( K ) )
     $      M = M + 1
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
         INFO = -1
      ELSE IF( .NOT.FROMQR .AND. .NOT.<a name="LSAME.205"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EIGSRC, <span class="string">'N'</span> ) ) THEN
         INFO = -2
      ELSE IF( .NOT.NOINIT .AND. .NOT.<a name="LSAME.207"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( INITV, <span class="string">'U'</span> ) ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
         INFO = -10
      ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
         INFO = -12
      ELSE IF( MM.LT.M ) THEN
         INFO = -13
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.221"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CHSEIN.221"></a><a href="chsein.f.html#CHSEIN.1">CHSEIN</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set machine-dependent constants.
</span><span class="comment">*</span><span class="comment">
</span>      UNFL = <a name="SLAMCH.232"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
      ULP = <a name="SLAMCH.233"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = UNFL*( N / ULP )
<span class="comment">*</span><span class="comment">
</span>      LDWORK = N
<span class="comment">*</span><span class="comment">
</span>      KL = 1
      KLN = 0
      IF( FROMQR ) THEN
         KR = 0
      ELSE
         KR = N
      END IF
      KS = 1
<span class="comment">*</span><span class="comment">
</span>      DO 100 K = 1, N
         IF( SELECT( K ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute eigenvector(s) corresponding to W(K).
</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="CLANHS.287"></a><a href="clanhs.f.html#CLANHS.1">CLANHS</a>( <span class="string">'I'</span>, KR-KL+1, H( KL, KL ), LDH, RWORK )
               IF( HNORM.GT.RZERO ) 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>            WK = W( K )
   60       CONTINUE
            DO 70 I = K - 1, KL, -1
               IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
                  WK = WK + EPS3
                  GO TO 60
               END IF
   70       CONTINUE
            W( K ) = WK
<span class="comment">*</span><span class="comment">
</span>            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="CLAEIN.313"></a><a href="claein.f.html#CLAEIN.1">CLAEIN</a>( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
     $                      WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
     $                      SMLNUM, IINFO )
               IF( IINFO.GT.0 ) THEN
                  INFO = INFO + 1
                  IFAILL( KS ) = K
               ELSE
                  IFAILL( KS ) = 0
               END IF
               DO 80 I = 1, KL - 1
                  VL( I, KS ) = ZERO
   80          CONTINUE
            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="CLAEIN.330"></a><a href="claein.f.html#CLAEIN.1">CLAEIN</a>( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
     $                      WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
               IF( IINFO.GT.0 ) THEN
                  INFO = INFO + 1
                  IFAILR( KS ) = K
               ELSE
                  IFAILR( KS ) = 0
               END IF
               DO 90 I = KR + 1, N
                  VR( I, KS ) = ZERO
   90          CONTINUE
            END IF
            KS = KS + 1
         END IF
  100 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="CHSEIN.348"></a><a href="chsein.f.html#CHSEIN.1">CHSEIN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>