dtgsen.f.html

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

</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK. LWORK &gt;=  4*N+16.
</span><span class="comment">*</span><span class="comment">          If IJOB = 1, 2 or 4, LWORK &gt;= MAX(4*N+16, 2*M*(N-M)).
</span><span class="comment">*</span><span class="comment">          If IJOB = 3 or 5, LWORK &gt;= MAX(4*N+16, 4*M*(N-M)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.178"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
</span><span class="comment">*</span><span class="comment">          IF IJOB = 0, IWORK is not referenced.  Otherwise,
</span><span class="comment">*</span><span class="comment">          on exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LIWORK  (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array IWORK. LIWORK &gt;= 1.
</span><span class="comment">*</span><span class="comment">          If IJOB = 1, 2 or 4, LIWORK &gt;=  N+6.
</span><span class="comment">*</span><span class="comment">          If IJOB = 3 or 5, LIWORK &gt;= MAX(2*M*(N-M), N+6).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LIWORK = -1, then a workspace query is assumed; the
</span><span class="comment">*</span><span class="comment">          routine only calculates the optimal size of the IWORK array,
</span><span class="comment">*</span><span class="comment">          returns this value as the first entry of the IWORK array, and
</span><span class="comment">*</span><span class="comment">          no error message related to LIWORK is issued by <a name="XERBLA.192"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</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">            =1: Reordering of (A, B) failed because the transformed
</span><span class="comment">*</span><span class="comment">                matrix pair (A, B) would be too far from generalized
</span><span class="comment">*</span><span class="comment">                Schur form; the problem is very ill-conditioned.
</span><span class="comment">*</span><span class="comment">                (A, B) may have been partially reordered.
</span><span class="comment">*</span><span class="comment">                If requested, 0 is returned in DIF(*), PL and PR.
</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">  <a name="DTGSEN.206"></a><a href="dtgsen.f.html#DTGSEN.1">DTGSEN</a> first collects the selected eigenvalues by computing
</span><span class="comment">*</span><span class="comment">  orthogonal U and W that move them to the top left corner of (A, B).
</span><span class="comment">*</span><span class="comment">  In other words, the selected eigenvalues are the eigenvalues of
</span><span class="comment">*</span><span class="comment">  (A11, B11) in:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                U'*(A, B)*W = (A11 A12) (B11 B12) n1
</span><span class="comment">*</span><span class="comment">                              ( 0  A22),( 0  B22) n2
</span><span class="comment">*</span><span class="comment">                                n1  n2    n1  n2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where N = n1+n2 and U' means the transpose of U. The first n1 columns
</span><span class="comment">*</span><span class="comment">  of U and W span the specified pair of left and right eigenspaces
</span><span class="comment">*</span><span class="comment">  (deflating subspaces) of (A, B).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If (A, B) has been obtained from the generalized real Schur
</span><span class="comment">*</span><span class="comment">  decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the
</span><span class="comment">*</span><span class="comment">  reordered generalized real Schur form of (C, D) is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  and the first n1 columns of Q*U and Z*W span the corresponding
</span><span class="comment">*</span><span class="comment">  deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Note that if the selected eigenvalue is sufficiently ill-conditioned,
</span><span class="comment">*</span><span class="comment">  then its value may differ significantly from its value before
</span><span class="comment">*</span><span class="comment">  reordering.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The reciprocal condition numbers of the left and right eigenspaces
</span><span class="comment">*</span><span class="comment">  spanned by the first n1 columns of U and W (or Q*U and Z*W) may
</span><span class="comment">*</span><span class="comment">  be returned in DIF(1:2), corresponding to Difu and Difl, resp.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The Difu and Difl are defined as:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
</span><span class="comment">*</span><span class="comment">  and
</span><span class="comment">*</span><span class="comment">       Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where sigma-min(Zu) is the smallest singular value of the
</span><span class="comment">*</span><span class="comment">  (2*n1*n2)-by-(2*n1*n2) matrix
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       Zu = [ kron(In2, A11)  -kron(A22', In1) ]
</span><span class="comment">*</span><span class="comment">            [ kron(In2, B11)  -kron(B22', In1) ].
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Here, Inx is the identity matrix of size nx and A22' is the
</span><span class="comment">*</span><span class="comment">  transpose of A22. kron(X, Y) is the Kronecker product between
</span><span class="comment">*</span><span class="comment">  the matrices X and Y.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  When DIF(2) is small, small changes in (A, B) can cause large changes
</span><span class="comment">*</span><span class="comment">  in the deflating subspace. An approximate (asymptotic) bound on the
</span><span class="comment">*</span><span class="comment">  maximum angular error in the computed deflating subspaces is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       EPS * norm((A, B)) / DIF(2),
</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 norm of the projectors on the left and right
</span><span class="comment">*</span><span class="comment">  eigenspaces associated with (A11, B11) may be returned in PL and PR.
</span><span class="comment">*</span><span class="comment">  They are computed as follows. First we compute L and R so that
</span><span class="comment">*</span><span class="comment">  P*(A, B)*Q is block diagonal, where
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       P = ( I -L ) n1           Q = ( I R ) n1
</span><span class="comment">*</span><span class="comment">           ( 0  I ) n2    and        ( 0 I ) n2
</span><span class="comment">*</span><span class="comment">             n1 n2                    n1 n2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  and (L, R) is the solution to the generalized Sylvester equation
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       A11*R - L*A22 = -A12
</span><span class="comment">*</span><span class="comment">       B11*R - L*B22 = -B12
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
</span><span class="comment">*</span><span class="comment">  An approximate (asymptotic) bound on the average absolute error of
</span><span class="comment">*</span><span class="comment">  the selected eigenvalues is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       EPS * norm((A, B)) / PL.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  There are also global error bounds which valid for perturbations up
</span><span class="comment">*</span><span class="comment">  to a certain restriction:  A lower bound (x) on the smallest
</span><span class="comment">*</span><span class="comment">  F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
</span><span class="comment">*</span><span class="comment">  coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
</span><span class="comment">*</span><span class="comment">  (i.e. (A + E, B + F), is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">   x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate bound on x can be computed from DIF(1:2), PL and PR.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If y = ( F-norm(E,F) / x) &lt;= 1, the angles between the perturbed
</span><span class="comment">*</span><span class="comment">  (L', R') and unperturbed (L, R) left and right deflating subspaces
</span><span class="comment">*</span><span class="comment">  associated with the selected cluster in the (1,1)-blocks can be
</span><span class="comment">*</span><span class="comment">  bounded as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">   max-angle(L, L') &lt;= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
</span><span class="comment">*</span><span class="comment">   max-angle(R, R') &lt;= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  See LAPACK User's Guide section 4.11 or the following references
</span><span class="comment">*</span><span class="comment">  for more information.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Note that if the default method for computing the Frobenius-norm-
</span><span class="comment">*</span><span class="comment">  based estimate DIF is not wanted (see <a name="DLATDF.302"></a><a href="dlatdf.f.html#DLATDF.1">DLATDF</a>), then the parameter
</span><span class="comment">*</span><span class="comment">  IDIFJB (see below) should be changed from 3 to 4 (routine <a name="DLATDF.303"></a><a href="dlatdf.f.html#DLATDF.1">DLATDF</a>
</span><span class="comment">*</span><span class="comment">  (IJOB = 2 will be used)). See <a name="DTGSYL.304"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</a> for more details.
</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,
</span><span class="comment">*</span><span class="comment">      Report UMINF - 94.04, Department of Computing Science, Umea
</span><span class="comment">*</span><span class="comment">      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
</span><span class="comment">*</span><span class="comment">      Note 87. 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. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
</span><span class="comment">*</span><span class="comment">      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>      INTEGER            IDIFJB
      PARAMETER          ( IDIFJB = 3 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2,
     $                   WANTP
      INTEGER            I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN,
     $                   MN2, N1, N2
      DOUBLE PRECISION   DSCALE, DSUM, EPS, RDSCAL, SMLNUM
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLACN2.352"></a><a href="dlacn2.f.html#DLACN2.1">DLACN2</a>, <a name="DLACPY.352"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>, <a name="DLAG2.352"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>, <a name="DLASSQ.352"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>, <a name="DTGEXC.352"></a><a href="dtgexc.f.html#DTGEXC.1">DTGEXC</a>, <a name="DTGSYL.352"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</a>,
     $                   <a name="XERBLA.353"></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">     .. External Functions ..