ztrsna.f.html

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

      COMPLEX*16         CDUM, PROD
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
      COMPLEX*16         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            IZAMAX
      DOUBLE PRECISION   <a name="DLAMCH.181"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, DZNRM2
      COMPLEX*16         ZDOTC
      EXTERNAL           <a name="LSAME.183"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, IZAMAX, <a name="DLAMCH.183"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, DZNRM2, ZDOTC
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="XERBLA.186"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, <a name="ZDRSCL.186"></a><a href="zdrscl.f.html#ZDRSCL.1">ZDRSCL</a>, <a name="ZLACN2.186"></a><a href="zlacn2.f.html#ZLACN2.1">ZLACN2</a>, <a name="ZLACPY.186"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLATRS.186"></a><a href="zlatrs.f.html#ZLATRS.1">ZLATRS</a>, <a name="ZTREXC.186"></a><a href="ztrexc.f.html#ZTREXC.1">ZTREXC</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, DIMAG, MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      DOUBLE PRECISION   CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( 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.201"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.202"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTSP = <a name="LSAME.203"></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.205"></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.223"></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.239"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZTRSNA.239"></a><a href="ztrsna.f.html#ZTRSNA.1">ZTRSNA</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="DLAMCH.262"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="DLAMCH.263"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="DLABAD.265"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</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 = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
            RNRM = DZNRM2( N, VR( 1, KS ), 1 )
            LNRM = DZNRM2( 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="ZLACPY.295"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, N, N, T, LDT, WORK, LDWORK )
            CALL <a name="ZTREXC.296"></a><a href="ztrexc.f.html#ZTREXC.1">ZTREXC</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="ZLACN2.312"></a><a href="zlacn2.f.html#ZLACN2.1">ZLACN2</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="ZLATRS.319"></a><a href="zlatrs.f.html#ZLATRS.1">ZLATRS</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="ZLATRS.326"></a><a href="zlatrs.f.html#ZLATRS.1">ZLATRS</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 = IZAMAX( N-1, WORK, 1 )
                  XNORM = CABS1( WORK( IX, 1 ) )
                  IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
     $               GO TO 40
                  CALL <a name="ZDRSCL.340"></a><a href="zdrscl.f.html#ZDRSCL.1">ZDRSCL</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="ZTRSNA.353"></a><a href="ztrsna.f.html#ZTRSNA.1">ZTRSNA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>