ctgsen.f.html

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

</span><span class="comment">*</span><span class="comment">          estimates of Difu and Difl, computed using reversed
</span><span class="comment">*</span><span class="comment">          communication with <a name="CLACN2.152"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>.
</span><span class="comment">*</span><span class="comment">          If M = 0 or N, DIF(1:2) = F-norm([A, B]).
</span><span class="comment">*</span><span class="comment">          If IJOB = 0 or 1, DIF is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          IF IJOB = 0, WORK is not referenced.  Otherwise,
</span><span class="comment">*</span><span class="comment">          on exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</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;=  1
</span><span class="comment">*</span><span class="comment">          If IJOB = 1, 2 or 4, LWORK &gt;=  2*M*(N-M)
</span><span class="comment">*</span><span class="comment">          If IJOB = 3 or 5, LWORK &gt;=  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.168"></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+2;
</span><span class="comment">*</span><span class="comment">          If IJOB = 3 or 5, LIWORK &gt;= MAX(N+2, 2*M*(N-M));
</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.182"></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">
</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="CTGSEN.197"></a><a href="ctgsen.f.html#CTGSEN.1">CTGSEN</a> first collects the selected eigenvalues by computing unitary
</span><span class="comment">*</span><span class="comment">  U and W that move them to the top left corner of (A, B). In other
</span><span class="comment">*</span><span class="comment">  words, the selected eigenvalues are the eigenvalues of (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 conjugate transpose of U. The first
</span><span class="comment">*</span><span class="comment">  n1 columns of U and W span the specified pair of left and right
</span><span class="comment">*</span><span class="comment">  eigenspaces (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 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="CLATDF.292"></a><a href="clatdf.f.html#CLATDF.1">CLATDF</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="CLATDF.293"></a><a href="clatdf.f.html#CLATDF.1">CLATDF</a>
</span><span class="comment">*</span><span class="comment">  (IJOB = 2 will be used)). See <a name="CTGSYL.294"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</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, 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