dtgsna.f.html

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

</span><span class="comment">*</span><span class="comment">              Complex eigenvalue pair.
</span><span class="comment">*</span><span class="comment">
</span>               RNRM = <a name="DLAPY2.440"></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.442"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( DNRM2( N, VL( 1, KS ), 1 ),
     $                DNRM2( N, VL( 1, KS+1 ), 1 ) )
               CALL DGEMV( <span class="string">'N'</span>, N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
     $                     WORK, 1 )
               TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
               TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
               CALL DGEMV( <span class="string">'N'</span>, N, N, ONE, A, LDA, VR( 1, KS+1 ), 1,
     $                     ZERO, WORK, 1 )
               TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
               TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
               UHAV = TMPRR + TMPII
               UHAVI = TMPIR - TMPRI
               CALL DGEMV( <span class="string">'N'</span>, N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
     $                     WORK, 1 )
               TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
               TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
               CALL DGEMV( <span class="string">'N'</span>, N, N, ONE, B, LDB, VR( 1, KS+1 ), 1,
     $                     ZERO, WORK, 1 )
               TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
               TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
               UHBV = TMPRR + TMPII
               UHBVI = TMPIR - TMPRI
               UHAV = <a name="DLAPY2.464"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( UHAV, UHAVI )
               UHBV = <a name="DLAPY2.465"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( UHBV, UHBVI )
               COND = <a name="DLAPY2.466"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( UHAV, UHBV )
               S( KS ) = COND / ( RNRM*LNRM )
               S( KS+1 ) = S( KS )
<span class="comment">*</span><span class="comment">
</span>            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Real eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>               RNRM = DNRM2( N, VR( 1, KS ), 1 )
               LNRM = DNRM2( N, VL( 1, KS ), 1 )
               CALL DGEMV( <span class="string">'N'</span>, N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
     $                     WORK, 1 )
               UHAV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
               CALL DGEMV( <span class="string">'N'</span>, N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
     $                     WORK, 1 )
               UHBV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
               COND = <a name="DLAPY2.482"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( UHAV, UHBV )
               IF( COND.EQ.ZERO ) THEN
                  S( KS ) = -ONE
               ELSE
                  S( KS ) = COND / ( RNRM*LNRM )
               END IF
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( WANTDF ) THEN
            IF( N.EQ.1 ) THEN
               DIF( KS ) = <a name="DLAPY2.493"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( A( 1, 1 ), B( 1, 1 ) )
               GO TO 20
            END IF
<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>            IF( PAIR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy the  2-by 2 pencil beginning at (A(k,k), B(k, k)).
</span><span class="comment">*</span><span class="comment">              Compute the eigenvalue(s) at position K.
</span><span class="comment">*</span><span class="comment">
</span>               WORK( 1 ) = A( K, K )
               WORK( 2 ) = A( K+1, K )
               WORK( 3 ) = A( K, K+1 )
               WORK( 4 ) = A( K+1, K+1 )
               WORK( 5 ) = B( K, K )
               WORK( 6 ) = B( K+1, K )
               WORK( 7 ) = B( K, K+1 )
               WORK( 8 ) = B( K+1, K+1 )
               CALL <a name="DLAG2.512"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA,
     $                     DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI )
               ALPRQT = ONE
               C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA )
               C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI
               ROOT1 = C1 + SQRT( C1*C1-4.0D0*C2 )
               ROOT2 = C2 / ROOT1
               ROOT1 = ROOT1 / TWO
               COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) )
            END IF
<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 swap the
</span><span class="comment">*</span><span class="comment">           diagonal block beginning at A(k,k) to the (1,1) position.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLACPY.526"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, N, N, A, LDA, WORK, N )
            CALL <a name="DLACPY.527"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</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="DTGEXC.531"></a><a href="dtgexc.f.html#DTGEXC.1">DTGEXC</a>( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N,
     $                   DUMMY, 1, DUMMY1, 1, IFST, ILST,
     $                   WORK( N*N*2+1 ), LWORK-2*N*N, 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
               IF( WORK( 2 ).NE.ZERO )
     $            N1 = 2
               N2 = N - N1
               IF( N2.EQ.0 ) THEN
                  DIF( KS ) = COND
               ELSE
                  I = N*N + 1
                  IZ = 2*N*N + 1
                  CALL <a name="DTGSYL.557"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</a>( <span class="string">'N'</span>, DIFDRI, 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 ),
     $                         WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR )
<span class="comment">*</span><span class="comment">
</span>                  IF( PAIR )
     $               DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ),
     $                           COND )
               END IF
            END IF
            IF( PAIR )
     $         DIF( KS+1 ) = DIF( KS )
         END IF
         IF( PAIR )
     $      KS = KS + 1
<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="DTGSNA.578"></a><a href="dtgsna.f.html#DTGSNA.1">DTGSNA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>