ztgsna.f.html

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

</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0: Successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0: If INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The reciprocal of the condition number of the i-th generalized
</span><span class="comment">*</span><span class="comment">  eigenvalue w = (a, b) is defined as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v))
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where u and v are the right and left eigenvectors of (A, B)
</span><span class="comment">*</span><span class="comment">  corresponding to w; |z| denotes the absolute value of the complex
</span><span class="comment">*</span><span class="comment">  number, and norm(u) denotes the 2-norm of the vector u. The pair
</span><span class="comment">*</span><span class="comment">  (a, b) corresponds to an eigenvalue w = a/b (= v'Au/v'Bu) of the
</span><span class="comment">*</span><span class="comment">  matrix pair (A, B). If both a and b equal zero, then (A,B) is
</span><span class="comment">*</span><span class="comment">  singular and S(I) = -1 is returned.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound on the chordal distance between the i-th
</span><span class="comment">*</span><span class="comment">  computed generalized eigenvalue w and the corresponding exact
</span><span class="comment">*</span><span class="comment">  eigenvalue lambda is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          chord(w, lambda) &lt;=   EPS * norm(A, B) / S(I),
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where EPS is the machine precision.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The reciprocal of the condition number of the right eigenvector u
</span><span class="comment">*</span><span class="comment">  and left eigenvector v corresponding to the generalized eigenvalue w
</span><span class="comment">*</span><span class="comment">  is defined as follows. Suppose
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                   (A, B) = ( a   *  ) ( b  *  )  1
</span><span class="comment">*</span><span class="comment">                            ( 0  A22 ),( 0 B22 )  n-1
</span><span class="comment">*</span><span class="comment">                              1  n-1     1 n-1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Then the reciprocal condition number DIF(I) is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          Difl[(a, b), (A22, B22)]  = sigma-min( Zl )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where sigma-min(Zl) denotes the smallest singular value of
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Zl = [ kron(a, In-1) -kron(1, A22) ]
</span><span class="comment">*</span><span class="comment">              [ kron(b, In-1) -kron(1, B22) ].
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Here In-1 is the identity matrix of size n-1 and X' is the conjugate
</span><span class="comment">*</span><span class="comment">  transpose of X. kron(X, Y) is the Kronecker product between the
</span><span class="comment">*</span><span class="comment">  matrices X and Y.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  We approximate the smallest singular value of Zl with an upper
</span><span class="comment">*</span><span class="comment">  bound. This is done by <a name="ZLATDF.173"></a><a href="zlatdf.f.html#ZLATDF.1">ZLATDF</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound for a computed eigenvector VL(i) or
</span><span class="comment">*</span><span class="comment">  VR(i) is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      EPS * norm(A, B) / DIF(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  See ref. [2-3] for more details and further references.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
</span><span class="comment">*</span><span class="comment">     Umea University, S-901 87 Umea, Sweden.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  References
</span><span class="comment">*</span><span class="comment">  ==========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
</span><span class="comment">*</span><span class="comment">      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
</span><span class="comment">*</span><span class="comment">      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
</span><span class="comment">*</span><span class="comment">      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
</span><span class="comment">*</span><span class="comment">      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
</span><span class="comment">*</span><span class="comment">      Estimation: Theory, Algorithms and Software, Report
</span><span class="comment">*</span><span class="comment">      UMINF - 94.04, Department of Computing Science, Umea University,
</span><span class="comment">*</span><span class="comment">      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
</span><span class="comment">*</span><span class="comment">      To appear in Numerical Algorithms, 1996.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
</span><span class="comment">*</span><span class="comment">      for Solving the Generalized Sylvester Equation and Estimating the
</span><span class="comment">*</span><span class="comment">      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
</span><span class="comment">*</span><span class="comment">      Department of Computing Science, Umea University, S-901 87 Umea,
</span><span class="comment">*</span><span class="comment">      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
</span><span class="comment">*</span><span class="comment">      Note 75.
</span><span class="comment">*</span><span class="comment">      To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      INTEGER            IDIFJB
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, SOMCON, WANTBH, WANTDF, WANTS
      INTEGER            I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
      DOUBLE PRECISION   BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
      COMPLEX*16         YHAX, YHBX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      COMPLEX*16         DUMMY( 1 ), DUMMY1( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.226"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      DOUBLE PRECISION   <a name="DLAMCH.227"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLAPY2.227"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>, DZNRM2
      COMPLEX*16         ZDOTC
      EXTERNAL           <a name="LSAME.229"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.229"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLAPY2.229"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>, DZNRM2, ZDOTC
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLABAD.232"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, <a name="XERBLA.232"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, ZGEMV, <a name="ZLACPY.232"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZTGEXC.232"></a><a href="ztgexc.f.html#ZTGEXC.1">ZTGEXC</a>, <a name="ZTGSYL.232"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DCMPLX, MAX
<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.241"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.242"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTDF = <a name="LSAME.243"></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.245"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'S'</span> )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      LQUERY = ( LWORK.EQ.-1 )
<span class="comment">*</span><span class="comment">
</span>      IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.252"></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( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -8
      ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
         INFO = -10
      ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
         INFO = -12