dtrsna.f.html

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

</span>               PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
               RNRM = DNRM2( N, VR( 1, KS ), 1 )
               LNRM = DNRM2( N, VL( 1, KS ), 1 )
               S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Complex eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>               PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
               PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
     $                 1 )
               PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
               PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
     $                 1 )
               RNRM = <a name="DLAPY2.344"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( DNRM2( N, VR( 1, KS ), 1 ),
     $                DNRM2( N, VR( 1, KS+1 ), 1 ) )
               LNRM = <a name="DLAPY2.346"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( DNRM2( N, VL( 1, KS ), 1 ),
     $                DNRM2( N, VL( 1, KS+1 ), 1 ) )
               COND = <a name="DLAPY2.348"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( PROD1, PROD2 ) / ( RNRM*LNRM )
               S( KS ) = COND
               S( KS+1 ) = COND
            END IF
         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 diagonal
</span><span class="comment">*</span><span class="comment">           block beginning at T(k,k) to the (1,1) position.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLACPY.362"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, N, N, T, LDT, WORK, LDWORK )
            IFST = K
            ILST = 1
            CALL <a name="DTREXC.365"></a><a href="dtrexc.f.html#DTREXC.1">DTREXC</a>( <span class="string">'No Q'</span>, N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
     $                   WORK( 1, N+1 ), IERR )
<span class="comment">*</span><span class="comment">
</span>            IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Could not swap because blocks not well separated
</span><span class="comment">*</span><span class="comment">
</span>               SCALE = ONE
               EST = BIGNUM
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Reordering successful
</span><span class="comment">*</span><span class="comment">
</span>               IF( WORK( 2, 1 ).EQ.ZERO ) THEN
<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
                  N2 = 1
                  NN = N - 1
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Triangularize the 2 by 2 block by unitary
</span><span class="comment">*</span><span class="comment">                 transformation U = [  cs   i*ss ]
</span><span class="comment">*</span><span class="comment">                                    [ i*ss   cs  ].
</span><span class="comment">*</span><span class="comment">                 such that the (1,1) position of WORK is complex
</span><span class="comment">*</span><span class="comment">                 eigenvalue lambda with positive imaginary part. (2,2)
</span><span class="comment">*</span><span class="comment">                 position of WORK is the complex eigenvalue lambda
</span><span class="comment">*</span><span class="comment">                 with negative imaginary  part.
</span><span class="comment">*</span><span class="comment">
</span>                  MU = SQRT( ABS( WORK( 1, 2 ) ) )*
     $                 SQRT( ABS( WORK( 2, 1 ) ) )
                  DELTA = <a name="DLAPY2.399"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( MU, WORK( 2, 1 ) )
                  CS = MU / DELTA
                  SN = -WORK( 2, 1 ) / DELTA
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Form
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
</span><span class="comment">*</span><span class="comment">                                        [   mu                     ]
</span><span class="comment">*</span><span class="comment">                                        [         ..               ]
</span><span class="comment">*</span><span class="comment">                                        [             ..           ]
</span><span class="comment">*</span><span class="comment">                                        [                  mu      ]
</span><span class="comment">*</span><span class="comment">                 where C' is conjugate transpose of complex matrix C,
</span><span class="comment">*</span><span class="comment">                 and RWORK is stored starting in the N+1-st column of
</span><span class="comment">*</span><span class="comment">                 WORK.
</span><span class="comment">*</span><span class="comment">
</span>                  DO 30 J = 3, N
                     WORK( 2, J ) = CS*WORK( 2, J )
                     WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
   30             CONTINUE
                  WORK( 2, 2 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>                  WORK( 1, N+1 ) = TWO*MU
                  DO 40 I = 2, N - 1
                     WORK( I, N+1 ) = SN*WORK( 1, I+1 )
   40             CONTINUE
                  N2 = 2
                  NN = 2*( N-1 )
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Estimate norm(inv(C'))
</span><span class="comment">*</span><span class="comment">
</span>               EST = ZERO
               KASE = 0
   50          CONTINUE
               CALL <a name="DLACN2.433"></a><a href="dlacn2.f.html#DLACN2.1">DLACN2</a>( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
     $                      EST, KASE, ISAVE )
               IF( KASE.NE.0 ) THEN
                  IF( KASE.EQ.1 ) THEN
                     IF( N2.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Real eigenvalue: solve C'*x = scale*c.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="DLAQTR.441"></a><a href="dlaqtr.f.html#DLAQTR.1">DLAQTR</a>( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
     $                               LDWORK, DUMMY, DUMM, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
                     ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Complex eigenvalue: solve
</span><span class="comment">*</span><span class="comment">                       C'*(p+iq) = scale*(c+id) in real arithmetic.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="DLAQTR.450"></a><a href="dlaqtr.f.html#DLAQTR.1">DLAQTR</a>( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
     $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
                     END IF
                  ELSE
                     IF( N2.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Real eigenvalue: solve C*x = scale*c.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="DLAQTR.460"></a><a href="dlaqtr.f.html#DLAQTR.1">DLAQTR</a>( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
     $                               LDWORK, DUMMY, DUMM, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
                     ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Complex eigenvalue: solve
</span><span class="comment">*</span><span class="comment">                       C*(p+iq) = scale*(c+id) in real arithmetic.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="DLAQTR.469"></a><a href="dlaqtr.f.html#DLAQTR.1">DLAQTR</a>( .FALSE., .FALSE., N-1,
     $                               WORK( 2, 2 ), LDWORK,
     $                               WORK( 1, N+1 ), MU, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
<span class="comment">*</span><span class="comment">
</span>                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  GO TO 50
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            SEP( KS ) = SCALE / MAX( EST, SMLNUM )
            IF( PAIR )
     $         SEP( KS+1 ) = SEP( KS )
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( PAIR )
     $      KS = KS + 1
<span class="comment">*</span><span class="comment">
</span>   60 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DTRSNA.493"></a><a href="dtrsna.f.html#DTRSNA.1">DTRSNA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>