ctgsna.f.html

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

      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set M to the number of eigenpairs for which condition numbers
</span><span class="comment">*</span><span class="comment">        are required, and test MM.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SOMCON ) THEN
            M = 0
            DO 10 K = 1, N
               IF( SELECT( K ) )
     $            M = M + 1
   10       CONTINUE
         ELSE
            M = N
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( N.EQ.0 ) THEN
            LWMIN = 1
         ELSE IF( <a name="LSAME.281"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'V'</span> ) .OR. <a name="LSAME.281"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> ) ) THEN
            LWMIN = 2*N*N
         ELSE
            LWMIN = N
         END IF
         WORK( 1 ) = LWMIN
<span class="comment">*</span><span class="comment">
</span>         IF( MM.LT.M ) THEN
            INFO = -15
         ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -18
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.296"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CTGSNA.296"></a><a href="ctgsna.f.html#CTGSNA.1">CTGSNA</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         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">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.309"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="SLAMCH.310"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> ) / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="SLABAD.312"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
      KS = 0
      DO 20 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine whether condition numbers are required for the k-th
</span><span class="comment">*</span><span class="comment">        eigenpair.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SOMCON ) THEN
            IF( .NOT.SELECT( K ) )
     $         GO TO 20
         END IF
<span class="comment">*</span><span class="comment">
</span>         KS = KS + 1
<span class="comment">*</span><span class="comment">
</span>         IF( WANTS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute the reciprocal condition number of the k-th
</span><span class="comment">*</span><span class="comment">           eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>            RNRM = SCNRM2( N, VR( 1, KS ), 1 )
            LNRM = SCNRM2( N, VL( 1, KS ), 1 )
            CALL CGEMV( <span class="string">'N'</span>, N, N, CMPLX( ONE, ZERO ), A, LDA,
     $                  VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 )
            YHAX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 )
            CALL CGEMV( <span class="string">'N'</span>, N, N, CMPLX( ONE, ZERO ), B, LDB,
     $                  VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 )
            YHBX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 )
            COND = <a name="SLAPY2.339"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( ABS( YHAX ), ABS( YHBX ) )
            IF( COND.EQ.ZERO ) THEN
               S( KS ) = -ONE
            ELSE
               S( KS ) = COND / ( RNRM*LNRM )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( WANTDF ) THEN
            IF( N.EQ.1 ) THEN
               DIF( KS ) = <a name="SLAPY2.349"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Estimate the reciprocal condition number of the k-th
</span><span class="comment">*</span><span class="comment">              eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy the matrix (A, B) to the array WORK and move the
</span><span class="comment">*</span><span class="comment">              (k,k)th pair to the (1,1) position.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLACPY.358"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'Full'</span>, N, N, A, LDA, WORK, N )
               CALL <a name="CLACPY.359"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'Full'</span>, N, N, B, LDB, WORK( N*N+1 ), N )
               IFST = K
               ILST = 1
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="CTGEXC.363"></a><a href="ctgexc.f.html#CTGEXC.1">CTGEXC</a>( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
     $                      N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
<span class="comment">*</span><span class="comment">
</span>               IF( IERR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Ill-conditioned problem - swap rejected.
</span><span class="comment">*</span><span class="comment">
</span>                  DIF( KS ) = ZERO
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Reordering successful, solve generalized Sylvester
</span><span class="comment">*</span><span class="comment">                 equation for R and L,
</span><span class="comment">*</span><span class="comment">                            A22 * R - L * A11 = A12
</span><span class="comment">*</span><span class="comment">                            B22 * R - L * B11 = B12,
</span><span class="comment">*</span><span class="comment">                 and compute estimate of Difl[(A11,B11), (A22, B22)].
</span><span class="comment">*</span><span class="comment">
</span>                  N1 = 1
                  N2 = N - N1
                  I = N*N + 1
                  CALL <a name="CTGSYL.382"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'N'</span>, IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
     $                         N, WORK, N, WORK( N1+1 ), N,
     $                         WORK( N*N1+N1+I ), N, WORK( I ), N,
     $                         WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
     $                         1, IWORK, IERR )
               END IF
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>   20 CONTINUE
      WORK( 1 ) = LWMIN
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CTGSNA.395"></a><a href="ctgsna.f.html#CTGSNA.1">CTGSNA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>