ctrsna.f.html

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

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
      COMPLEX            DUMMY( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.179"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            ICAMAX
      REAL               SCNRM2, <a name="SLAMCH.181"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      COMPLEX            CDOTC
      EXTERNAL           <a name="LSAME.183"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, ICAMAX, SCNRM2, <a name="SLAMCH.183"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, CDOTC
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="CLACN2.186"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>, <a name="CLACPY.186"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLATRS.186"></a><a href="clatrs.f.html#CLATRS.1">CLATRS</a>, <a name="CSRSCL.186"></a><a href="csrscl.f.html#CSRSCL.1">CSRSCL</a>, <a name="CTREXC.186"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a>, <a name="SLABAD.186"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>,
     $                   <a name="XERBLA.187"></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>      WANTBH = <a name="LSAME.202"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.203"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTSP = <a name="LSAME.204"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'V'</span> ) .OR. WANTBH
<span class="comment">*</span><span class="comment">
</span>      SOMCON = <a name="LSAME.206"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'S'</span> )
<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 are
</span><span class="comment">*</span><span class="comment">     to be computed.
</span><span class="comment">*</span><span class="comment">
</span>      IF( SOMCON ) THEN
         M = 0
         DO 10 J = 1, N
            IF( SELECT( J ) )
     $         M = M + 1
   10    CONTINUE
      ELSE
         M = N
      END IF
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.224"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'A'</span> ) .AND. .NOT.SOMCON ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
         INFO = -10
      ELSE IF( MM.LT.M ) THEN
         INFO = -13
      ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
         INFO = -16
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.240"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CTRSNA.240"></a><a href="ctrsna.f.html#CTRSNA.1">CTRSNA</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>      IF( N.EQ.1 ) THEN
         IF( SOMCON ) THEN
            IF( .NOT.SELECT( 1 ) )
     $         RETURN
         END IF
         IF( WANTS )
     $      S( 1 ) = ONE
         IF( WANTSP )
     $      SEP( 1 ) = ABS( T( 1, 1 ) )
         RETURN
      END IF
<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.263"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="SLAMCH.264"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> ) / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="SLABAD.266"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
<span class="comment">*</span><span class="comment">
</span>      KS = 1
      DO 50 K = 1, N
<span class="comment">*</span><span class="comment">
</span>         IF( SOMCON ) THEN
            IF( .NOT.SELECT( K ) )
     $         GO TO 50
         END IF
<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>            PROD = CDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
            RNRM = SCNRM2( N, VR( 1, KS ), 1 )
            LNRM = SCNRM2( N, VL( 1, KS ), 1 )
            S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( WANTSP ) THEN
<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">           eigenvector.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy the matrix T to the array WORK and swap the k-th
</span><span class="comment">*</span><span class="comment">           diagonal element to the (1,1) position.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLACPY.296"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'Full'</span>, N, N, T, LDT, WORK, LDWORK )
            CALL <a name="CTREXC.297"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a>( <span class="string">'No Q'</span>, N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Form  C = T22 - lambda*I in WORK(2:N,2:N).
</span><span class="comment">*</span><span class="comment">
</span>            DO 20 I = 2, N
               WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
   20       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Estimate a lower bound for the 1-norm of inv(C'). The 1st
</span><span class="comment">*</span><span class="comment">           and (N+1)th columns of WORK are used to store work vectors.
</span><span class="comment">*</span><span class="comment">
</span>            SEP( KS ) = ZERO
            EST = ZERO
            KASE = 0
            NORMIN = <span class="string">'N'</span>
   30       CONTINUE
            CALL <a name="CLACN2.313"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
<span class="comment">*</span><span class="comment">
</span>            IF( KASE.NE.0 ) THEN
               IF( KASE.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve C'*x = scale*b
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLATRS.320"></a><a href="clatrs.f.html#CLATRS.1">CLATRS</a>( <span class="string">'Upper'</span>, <span class="string">'Conjugate transpose'</span>,
     $                         <span class="string">'Nonunit'</span>, NORMIN, N-1, WORK( 2, 2 ),
     $                         LDWORK, WORK, SCALE, RWORK, IERR )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve C*x = scale*b
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLATRS.327"></a><a href="clatrs.f.html#CLATRS.1">CLATRS</a>( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Nonunit'</span>,
     $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
     $                         SCALE, RWORK, IERR )
               END IF
               NORMIN = <span class="string">'Y'</span>
               IF( SCALE.NE.ONE ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Multiply by 1/SCALE if doing so will not cause
</span><span class="comment">*</span><span class="comment">                 overflow.
</span><span class="comment">*</span><span class="comment">
</span>                  IX = ICAMAX( N-1, WORK, 1 )
                  XNORM = CABS1( WORK( IX, 1 ) )
                  IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
     $               GO TO 40
                  CALL <a name="CSRSCL.341"></a><a href="csrscl.f.html#CSRSCL.1">CSRSCL</a>( N, SCALE, WORK, 1 )
               END IF
               GO TO 30
            END IF
<span class="comment">*</span><span class="comment">
</span>            SEP( KS ) = ONE / MAX( EST, SMLNUM )
         END IF
<span class="comment">*</span><span class="comment">
</span>   40    CONTINUE
         KS = KS + 1
   50 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CTRSNA.354"></a><a href="ctrsna.f.html#CTRSNA.1">CTRSNA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>