COMPLEX routines for triangular matrices, generalized problem (i.e., a pair of triangular matrices) matrix
ctgevc
USAGE:
m, info, vl, vr = NumRu::Lapack.ctgevc( side, howmny, select, s, p, vl, vr, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
* Purpose
* =======
*
* CTGEVC computes some or all of the right and/or left eigenvectors of
* a pair of complex matrices (S,P), where S and P are upper triangular.
* Matrix pairs of this type are produced by the generalized Schur
* factorization of a complex matrix pair (A,B):
*
* A = Q*S*Z**H, B = Q*P*Z**H
*
* as computed by CGGHRD + CHGEQZ.
*
* The right eigenvector x and the left eigenvector y of (S,P)
* corresponding to an eigenvalue w are defined by:
*
* S*x = w*P*x, (y**H)*S = w*(y**H)*P,
*
* where y**H denotes the conjugate tranpose of y.
* The eigenvalues are not input to this routine, but are computed
* directly from the diagonal elements of S and P.
*
* This routine returns the matrices X and/or Y of right and left
* eigenvectors of (S,P), or the products Z*X and/or Q*Y,
* where Z and Q are input matrices.
* If Q and Z are the unitary factors from the generalized Schur
* factorization of a matrix pair (A,B), then Z*X and Q*Y
* are the matrices of right and left eigenvectors of (A,B).
*
* Arguments
* =========
*
* SIDE (input) CHARACTER*1
* = 'R': compute right eigenvectors only;
* = 'L': compute left eigenvectors only;
* = 'B': compute both right and left eigenvectors.
*
* HOWMNY (input) CHARACTER*1
* = 'A': compute all right and/or left eigenvectors;
* = 'B': compute all right and/or left eigenvectors,
* backtransformed by the matrices in VR and/or VL;
* = 'S': compute selected right and/or left eigenvectors,
* specified by the logical array SELECT.
*
* SELECT (input) LOGICAL array, dimension (N)
* If HOWMNY='S', SELECT specifies the eigenvectors to be
* computed. The eigenvector corresponding to the j-th
* eigenvalue is computed if SELECT(j) = .TRUE..
* Not referenced if HOWMNY = 'A' or 'B'.
*
* N (input) INTEGER
* The order of the matrices S and P. N >= 0.
*
* S (input) COMPLEX array, dimension (LDS,N)
* The upper triangular matrix S from a generalized Schur
* factorization, as computed by CHGEQZ.
*
* LDS (input) INTEGER
* The leading dimension of array S. LDS >= max(1,N).
*
* P (input) COMPLEX array, dimension (LDP,N)
* The upper triangular matrix P from a generalized Schur
* factorization, as computed by CHGEQZ. P must have real
* diagonal elements.
*
* LDP (input) INTEGER
* The leading dimension of array P. LDP >= max(1,N).
*
* VL (input/output) COMPLEX array, dimension (LDVL,MM)
* On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
* contain an N-by-N matrix Q (usually the unitary matrix Q
* of left Schur vectors returned by CHGEQZ).
* On exit, if SIDE = 'L' or 'B', VL contains:
* if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
* if HOWMNY = 'B', the matrix Q*Y;
* if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
* SELECT, stored consecutively in the columns of
* VL, in the same order as their eigenvalues.
* Not referenced if SIDE = 'R'.
*
* LDVL (input) INTEGER
* The leading dimension of array VL. LDVL >= 1, and if
* SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
*
* VR (input/output) COMPLEX array, dimension (LDVR,MM)
* On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
* contain an N-by-N matrix Q (usually the unitary matrix Z
* of right Schur vectors returned by CHGEQZ).
* On exit, if SIDE = 'R' or 'B', VR contains:
* if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
* if HOWMNY = 'B', the matrix Z*X;
* if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
* SELECT, stored consecutively in the columns of
* VR, in the same order as their eigenvalues.
* Not referenced if SIDE = 'L'.
*
* LDVR (input) INTEGER
* The leading dimension of the array VR. LDVR >= 1, and if
* SIDE = 'R' or 'B', LDVR >= N.
*
* MM (input) INTEGER
* The number of columns in the arrays VL and/or VR. MM >= M.
*
* M (output) INTEGER
* The number of columns in the arrays VL and/or VR actually
* used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
* is set to N. Each selected eigenvector occupies one column.
*
* WORK (workspace) COMPLEX array, dimension (2*N)
*
* RWORK (workspace) REAL array, dimension (2*N)
*
* INFO (output) INTEGER
* = 0: successful exit.
* < 0: if INFO = -i, the i-th argument had an illegal value.
*
* =====================================================================
*
go to the page top
ctgex2
USAGE:
info, a, b, q, z = NumRu::Lapack.ctgex2( wantq, wantz, a, b, q, ldq, z, ldz, j1, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, INFO )
* Purpose
* =======
*
* CTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
* in an upper triangular matrix pair (A, B) by an unitary equivalence
* transformation.
*
* (A, B) must be in generalized Schur canonical form, that is, A and
* B are both upper triangular.
*
* Optionally, the matrices Q and Z of generalized Schur vectors are
* updated.
*
* Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
* Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
*
*
* Arguments
* =========
*
* WANTQ (input) LOGICAL
* .TRUE. : update the left transformation matrix Q;
* .FALSE.: do not update Q.
*
* WANTZ (input) LOGICAL
* .TRUE. : update the right transformation matrix Z;
* .FALSE.: do not update Z.
*
* N (input) INTEGER
* The order of the matrices A and B. N >= 0.
*
* A (input/output) COMPLEX arrays, dimensions (LDA,N)
* On entry, the matrix A in the pair (A, B).
* On exit, the updated matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* B (input/output) COMPLEX arrays, dimensions (LDB,N)
* On entry, the matrix B in the pair (A, B).
* On exit, the updated matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* Q (input/output) COMPLEX array, dimension (LDZ,N)
* If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
* the updated matrix Q.
* Not referenced if WANTQ = .FALSE..
*
* LDQ (input) INTEGER
* The leading dimension of the array Q. LDQ >= 1;
* If WANTQ = .TRUE., LDQ >= N.
*
* Z (input/output) COMPLEX array, dimension (LDZ,N)
* If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
* the updated matrix Z.
* Not referenced if WANTZ = .FALSE..
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1;
* If WANTZ = .TRUE., LDZ >= N.
*
* J1 (input) INTEGER
* The index to the first block (A11, B11).
*
* INFO (output) INTEGER
* =0: Successful exit.
* =1: The transformed matrix pair (A, B) would be too far
* from generalized Schur form; the problem is ill-
* conditioned.
*
*
* Further Details
* ===============
*
* Based on contributions by
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
* Umea University, S-901 87 Umea, Sweden.
*
* In the current code both weak and strong stability tests are
* performed. The user can omit the strong stability test by changing
* the internal logical parameter WANDS to .FALSE.. See ref. [2] for
* details.
*
* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
* M.S. Moonen et al (eds), Linear Algebra for Large Scale and
* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
* Eigenvalues of a Regular Matrix Pair (A, B) and Condition
* Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
* Department of Computing Science, Umea University, S-901 87 Umea,
* Sweden, 1994. Also as LAPACK Working Note 87. To appear in
* Numerical Algorithms, 1996.
*
* =====================================================================
*
go to the page top
ctgexc
USAGE:
info, a, b, q, z, ilst = NumRu::Lapack.ctgexc( wantq, wantz, a, b, q, ldq, z, ifst, ilst, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO )
* Purpose
* =======
*
* CTGEXC reorders the generalized Schur decomposition of a complex
* matrix pair (A,B), using an unitary equivalence transformation
* (A, B) := Q * (A, B) * Z', so that the diagonal block of (A, B) with
* row index IFST is moved to row ILST.
*
* (A, B) must be in generalized Schur canonical form, that is, A and
* B are both upper triangular.
*
* Optionally, the matrices Q and Z of generalized Schur vectors are
* updated.
*
* Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
* Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
*
* Arguments
* =========
*
* WANTQ (input) LOGICAL
* .TRUE. : update the left transformation matrix Q;
* .FALSE.: do not update Q.
*
* WANTZ (input) LOGICAL
* .TRUE. : update the right transformation matrix Z;
* .FALSE.: do not update Z.
*
* N (input) INTEGER
* The order of the matrices A and B. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the upper triangular matrix A in the pair (A, B).
* On exit, the updated matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* B (input/output) COMPLEX array, dimension (LDB,N)
* On entry, the upper triangular matrix B in the pair (A, B).
* On exit, the updated matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* Q (input/output) COMPLEX array, dimension (LDZ,N)
* On entry, if WANTQ = .TRUE., the unitary matrix Q.
* On exit, the updated matrix Q.
* If WANTQ = .FALSE., Q is not referenced.
*
* LDQ (input) INTEGER
* The leading dimension of the array Q. LDQ >= 1;
* If WANTQ = .TRUE., LDQ >= N.
*
* Z (input/output) COMPLEX array, dimension (LDZ,N)
* On entry, if WANTZ = .TRUE., the unitary matrix Z.
* On exit, the updated matrix Z.
* If WANTZ = .FALSE., Z is not referenced.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1;
* If WANTZ = .TRUE., LDZ >= N.
*
* IFST (input) INTEGER
* ILST (input/output) INTEGER
* Specify the reordering of the diagonal blocks of (A, B).
* The block with row index IFST is moved to row ILST, by a
* sequence of swapping between adjacent blocks.
*
* INFO (output) INTEGER
* =0: Successful exit.
* <0: if INFO = -i, the i-th argument had an illegal value.
* =1: The transformed matrix pair (A, B) would be too far
* from generalized Schur form; the problem is ill-
* conditioned. (A, B) may have been partially reordered,
* and ILST points to the first row of the current
* position of the block being moved.
*
*
* Further Details
* ===============
*
* Based on contributions by
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
* Umea University, S-901 87 Umea, Sweden.
*
* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
* M.S. Moonen et al (eds), Linear Algebra for Large Scale and
* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
* Eigenvalues of a Regular Matrix Pair (A, B) and Condition
* Estimation: Theory, Algorithms and Software, Report
* UMINF - 94.04, Department of Computing Science, Umea University,
* S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
* To appear in Numerical Algorithms, 1996.
*
* [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
* for Solving the Generalized Sylvester Equation and Estimating the
* Separation between Regular Matrix Pairs, Report UMINF - 93.23,
* Department of Computing Science, Umea University, S-901 87 Umea,
* Sweden, December 1993, Revised April 1994, Also as LAPACK working
* Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
* 1996.
*
* =====================================================================
*
* .. Local Scalars ..
INTEGER HERE
* ..
* .. External Subroutines ..
EXTERNAL CTGEX2, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
go to the page top
ctgsen
USAGE:
alpha, beta, m, pl, pr, dif, work, iwork, info, a, b, q, z = NumRu::Lapack.ctgsen( ijob, wantq, wantz, select, a, b, q, z, [:lwork => lwork, :liwork => liwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO )
* Purpose
* =======
*
* CTGSEN reorders the generalized Schur decomposition of a complex
* matrix pair (A, B) (in terms of an unitary equivalence trans-
* formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues
* appears in the leading diagonal blocks of the pair (A,B). The leading
* columns of Q and Z form unitary bases of the corresponding left and
* right eigenspaces (deflating subspaces). (A, B) must be in
* generalized Schur canonical form, that is, A and B are both upper
* triangular.
*
* CTGSEN also computes the generalized eigenvalues
*
* w(j)= ALPHA(j) / BETA(j)
*
* of the reordered matrix pair (A, B).
*
* Optionally, the routine computes estimates of reciprocal condition
* numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
* (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
* between the matrix pairs (A11, B11) and (A22,B22) that correspond to
* the selected cluster and the eigenvalues outside the cluster, resp.,
* and norms of "projections" onto left and right eigenspaces w.r.t.
* the selected cluster in the (1,1)-block.
*
*
* Arguments
* =========
*
* IJOB (input) integer
* Specifies whether condition numbers are required for the
* cluster of eigenvalues (PL and PR) or the deflating subspaces
* (Difu and Difl):
* =0: Only reorder w.r.t. SELECT. No extras.
* =1: Reciprocal of norms of "projections" onto left and right
* eigenspaces w.r.t. the selected cluster (PL and PR).
* =2: Upper bounds on Difu and Difl. F-norm-based estimate
* (DIF(1:2)).
* =3: Estimate of Difu and Difl. 1-norm-based estimate
* (DIF(1:2)).
* About 5 times as expensive as IJOB = 2.
* =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
* version to get it all.
* =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
*
* WANTQ (input) LOGICAL
* .TRUE. : update the left transformation matrix Q;
* .FALSE.: do not update Q.
*
* WANTZ (input) LOGICAL
* .TRUE. : update the right transformation matrix Z;
* .FALSE.: do not update Z.
*
* SELECT (input) LOGICAL array, dimension (N)
* SELECT specifies the eigenvalues in the selected cluster. To
* select an eigenvalue w(j), SELECT(j) must be set to
* .TRUE..
*
* N (input) INTEGER
* The order of the matrices A and B. N >= 0.
*
* A (input/output) COMPLEX array, dimension(LDA,N)
* On entry, the upper triangular matrix A, in generalized
* Schur canonical form.
* On exit, A is overwritten by the reordered matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* B (input/output) COMPLEX array, dimension(LDB,N)
* On entry, the upper triangular matrix B, in generalized
* Schur canonical form.
* On exit, B is overwritten by the reordered matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* ALPHA (output) COMPLEX array, dimension (N)
* BETA (output) COMPLEX array, dimension (N)
* The diagonal elements of A and B, respectively,
* when the pair (A,B) has been reduced to generalized Schur
* form. ALPHA(i)/BETA(i) i=1,...,N are the generalized
* eigenvalues.
*
* Q (input/output) COMPLEX array, dimension (LDQ,N)
* On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.
* On exit, Q has been postmultiplied by the left unitary
* transformation matrix which reorder (A, B); The leading M
* columns of Q form orthonormal bases for the specified pair of
* left eigenspaces (deflating subspaces).
* If WANTQ = .FALSE., Q is not referenced.
*
* LDQ (input) INTEGER
* The leading dimension of the array Q. LDQ >= 1.
* If WANTQ = .TRUE., LDQ >= N.
*
* Z (input/output) COMPLEX array, dimension (LDZ,N)
* On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.
* On exit, Z has been postmultiplied by the left unitary
* transformation matrix which reorder (A, B); The leading M
* columns of Z form orthonormal bases for the specified pair of
* left eigenspaces (deflating subspaces).
* If WANTZ = .FALSE., Z is not referenced.
*
* LDZ (input) INTEGER
* The leading dimension of the array Z. LDZ >= 1.
* If WANTZ = .TRUE., LDZ >= N.
*
* M (output) INTEGER
* The dimension of the specified pair of left and right
* eigenspaces, (deflating subspaces) 0 <= M <= N.
*
* PL (output) REAL
* PR (output) REAL
* If IJOB = 1, 4 or 5, PL, PR are lower bounds on the
* reciprocal of the norm of "projections" onto left and right
* eigenspace with respect to the selected cluster.
* 0 < PL, PR <= 1.
* If M = 0 or M = N, PL = PR = 1.
* If IJOB = 0, 2 or 3 PL, PR are not referenced.
*
* DIF (output) REAL array, dimension (2).
* If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
* If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
* Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
* estimates of Difu and Difl, computed using reversed
* communication with CLACN2.
* If M = 0 or N, DIF(1:2) = F-norm([A, B]).
* If IJOB = 0 or 1, DIF is not referenced.
*
* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK. LWORK >= 1
* If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M)
* If IJOB = 3 or 5, LWORK >= 4*M*(N-M)
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
*
* LIWORK (input) INTEGER
* The dimension of the array IWORK. LIWORK >= 1.
* If IJOB = 1, 2 or 4, LIWORK >= N+2;
* If IJOB = 3 or 5, LIWORK >= MAX(N+2, 2*M*(N-M));
*
* If LIWORK = -1, then a workspace query is assumed; the
* routine only calculates the optimal size of the IWORK array,
* returns this value as the first entry of the IWORK array, and
* no error message related to LIWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* =0: Successful exit.
* <0: If INFO = -i, the i-th argument had an illegal value.
* =1: Reordering of (A, B) failed because the transformed
* matrix pair (A, B) would be too far from generalized
* Schur form; the problem is very ill-conditioned.
* (A, B) may have been partially reordered.
* If requested, 0 is returned in DIF(*), PL and PR.
*
*
* Further Details
* ===============
*
* CTGSEN first collects the selected eigenvalues by computing unitary
* U and W that move them to the top left corner of (A, B). In other
* words, the selected eigenvalues are the eigenvalues of (A11, B11) in
*
* U'*(A, B)*W = (A11 A12) (B11 B12) n1
* ( 0 A22),( 0 B22) n2
* n1 n2 n1 n2
*
* where N = n1+n2 and U' means the conjugate transpose of U. The first
* n1 columns of U and W span the specified pair of left and right
* eigenspaces (deflating subspaces) of (A, B).
*
* If (A, B) has been obtained from the generalized real Schur
* decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the
* reordered generalized Schur form of (C, D) is given by
*
* (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
*
* and the first n1 columns of Q*U and Z*W span the corresponding
* deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
*
* Note that if the selected eigenvalue is sufficiently ill-conditioned,
* then its value may differ significantly from its value before
* reordering.
*
* The reciprocal condition numbers of the left and right eigenspaces
* spanned by the first n1 columns of U and W (or Q*U and Z*W) may
* be returned in DIF(1:2), corresponding to Difu and Difl, resp.
*
* The Difu and Difl are defined as:
*
* Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
* and
* Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
*
* where sigma-min(Zu) is the smallest singular value of the
* (2*n1*n2)-by-(2*n1*n2) matrix
*
* Zu = [ kron(In2, A11) -kron(A22', In1) ]
* [ kron(In2, B11) -kron(B22', In1) ].
*
* Here, Inx is the identity matrix of size nx and A22' is the
* transpose of A22. kron(X, Y) is the Kronecker product between
* the matrices X and Y.
*
* When DIF(2) is small, small changes in (A, B) can cause large changes
* in the deflating subspace. An approximate (asymptotic) bound on the
* maximum angular error in the computed deflating subspaces is
*
* EPS * norm((A, B)) / DIF(2),
*
* where EPS is the machine precision.
*
* The reciprocal norm of the projectors on the left and right
* eigenspaces associated with (A11, B11) may be returned in PL and PR.
* They are computed as follows. First we compute L and R so that
* P*(A, B)*Q is block diagonal, where
*
* P = ( I -L ) n1 Q = ( I R ) n1
* ( 0 I ) n2 and ( 0 I ) n2
* n1 n2 n1 n2
*
* and (L, R) is the solution to the generalized Sylvester equation
*
* A11*R - L*A22 = -A12
* B11*R - L*B22 = -B12
*
* Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
* An approximate (asymptotic) bound on the average absolute error of
* the selected eigenvalues is
*
* EPS * norm((A, B)) / PL.
*
* There are also global error bounds which valid for perturbations up
* to a certain restriction: A lower bound (x) on the smallest
* F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
* coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
* (i.e. (A + E, B + F), is
*
* x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
*
* An approximate bound on x can be computed from DIF(1:2), PL and PR.
*
* If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
* (L', R') and unperturbed (L, R) left and right deflating subspaces
* associated with the selected cluster in the (1,1)-blocks can be
* bounded as
*
* max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
* max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
*
* See LAPACK User's Guide section 4.11 or the following references
* for more information.
*
* Note that if the default method for computing the Frobenius-norm-
* based estimate DIF is not wanted (see CLATDF), then the parameter
* IDIFJB (see below) should be changed from 3 to 4 (routine CLATDF
* (IJOB = 2 will be used)). See CTGSYL for more details.
*
* Based on contributions by
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
* Umea University, S-901 87 Umea, Sweden.
*
* References
* ==========
*
* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
* M.S. Moonen et al (eds), Linear Algebra for Large Scale and
* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
* Eigenvalues of a Regular Matrix Pair (A, B) and Condition
* Estimation: Theory, Algorithms and Software, Report
* UMINF - 94.04, Department of Computing Science, Umea University,
* S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
* To appear in Numerical Algorithms, 1996.
*
* [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
* for Solving the Generalized Sylvester Equation and Estimating the
* Separation between Regular Matrix Pairs, Report UMINF - 93.23,
* Department of Computing Science, Umea University, S-901 87 Umea,
* Sweden, December 1993, Revised April 1994, Also as LAPACK working
* Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
* 1996.
*
* =====================================================================
*
go to the page top
ctgsja
USAGE:
alpha, beta, ncycle, info, a, b, u, v, q = NumRu::Lapack.ctgsja( jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, NCYCLE, INFO )
* Purpose
* =======
*
* CTGSJA computes the generalized singular value decomposition (GSVD)
* of two complex upper triangular (or trapezoidal) matrices A and B.
*
* On entry, it is assumed that matrices A and B have the following
* forms, which may be obtained by the preprocessing subroutine CGGSVP
* from a general M-by-N matrix A and P-by-N matrix B:
*
* N-K-L K L
* A = K ( 0 A12 A13 ) if M-K-L >= 0;
* L ( 0 0 A23 )
* M-K-L ( 0 0 0 )
*
* N-K-L K L
* A = K ( 0 A12 A13 ) if M-K-L < 0;
* M-K ( 0 0 A23 )
*
* N-K-L K L
* B = L ( 0 0 B13 )
* P-L ( 0 0 0 )
*
* where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
* upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
* otherwise A23 is (M-K)-by-L upper trapezoidal.
*
* On exit,
*
* U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
*
* where U, V and Q are unitary matrices, Z' denotes the conjugate
* transpose of Z, R is a nonsingular upper triangular matrix, and D1
* and D2 are ``diagonal'' matrices, which are of the following
* structures:
*
* If M-K-L >= 0,
*
* K L
* D1 = K ( I 0 )
* L ( 0 C )
* M-K-L ( 0 0 )
*
* K L
* D2 = L ( 0 S )
* P-L ( 0 0 )
*
* N-K-L K L
* ( 0 R ) = K ( 0 R11 R12 ) K
* L ( 0 0 R22 ) L
*
* where
*
* C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
* S = diag( BETA(K+1), ... , BETA(K+L) ),
* C**2 + S**2 = I.
*
* R is stored in A(1:K+L,N-K-L+1:N) on exit.
*
* If M-K-L < 0,
*
* K M-K K+L-M
* D1 = K ( I 0 0 )
* M-K ( 0 C 0 )
*
* K M-K K+L-M
* D2 = M-K ( 0 S 0 )
* K+L-M ( 0 0 I )
* P-L ( 0 0 0 )
*
* N-K-L K M-K K+L-M
* ( 0 R ) = K ( 0 R11 R12 R13 )
* M-K ( 0 0 R22 R23 )
* K+L-M ( 0 0 0 R33 )
*
* where
* C = diag( ALPHA(K+1), ... , ALPHA(M) ),
* S = diag( BETA(K+1), ... , BETA(M) ),
* C**2 + S**2 = I.
*
* R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
* ( 0 R22 R23 )
* in B(M-K+1:L,N+M-K-L+1:N) on exit.
*
* The computation of the unitary transformation matrices U, V or Q
* is optional. These matrices may either be formed explicitly, or they
* may be postmultiplied into input matrices U1, V1, or Q1.
*
* Arguments
* =========
*
* JOBU (input) CHARACTER*1
* = 'U': U must contain a unitary matrix U1 on entry, and
* the product U1*U is returned;
* = 'I': U is initialized to the unit matrix, and the
* unitary matrix U is returned;
* = 'N': U is not computed.
*
* JOBV (input) CHARACTER*1
* = 'V': V must contain a unitary matrix V1 on entry, and
* the product V1*V is returned;
* = 'I': V is initialized to the unit matrix, and the
* unitary matrix V is returned;
* = 'N': V is not computed.
*
* JOBQ (input) CHARACTER*1
* = 'Q': Q must contain a unitary matrix Q1 on entry, and
* the product Q1*Q is returned;
* = 'I': Q is initialized to the unit matrix, and the
* unitary matrix Q is returned;
* = 'N': Q is not computed.
*
* M (input) INTEGER
* The number of rows of the matrix A. M >= 0.
*
* P (input) INTEGER
* The number of rows of the matrix B. P >= 0.
*
* N (input) INTEGER
* The number of columns of the matrices A and B. N >= 0.
*
* K (input) INTEGER
* L (input) INTEGER
* K and L specify the subblocks in the input matrices A and B:
* A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N)
* of A and B, whose GSVD is going to be computed by CTGSJA.
* See Further Details.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the M-by-N matrix A.
* On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular
* matrix R or part of R. See Purpose for details.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
*
* B (input/output) COMPLEX array, dimension (LDB,N)
* On entry, the P-by-N matrix B.
* On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
* a part of R. See Purpose for details.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,P).
*
* TOLA (input) REAL
* TOLB (input) REAL
* TOLA and TOLB are the convergence criteria for the Jacobi-
* Kogbetliantz iteration procedure. Generally, they are the
* same as used in the preprocessing step, say
* TOLA = MAX(M,N)*norm(A)*MACHEPS,
* TOLB = MAX(P,N)*norm(B)*MACHEPS.
*
* ALPHA (output) REAL array, dimension (N)
* BETA (output) REAL array, dimension (N)
* On exit, ALPHA and BETA contain the generalized singular
* value pairs of A and B;
* ALPHA(1:K) = 1,
* BETA(1:K) = 0,
* and if M-K-L >= 0,
* ALPHA(K+1:K+L) = diag(C),
* BETA(K+1:K+L) = diag(S),
* or if M-K-L < 0,
* ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
* BETA(K+1:M) = S, BETA(M+1:K+L) = 1.
* Furthermore, if K+L < N,
* ALPHA(K+L+1:N) = 0
* BETA(K+L+1:N) = 0.
*
* U (input/output) COMPLEX array, dimension (LDU,M)
* On entry, if JOBU = 'U', U must contain a matrix U1 (usually
* the unitary matrix returned by CGGSVP).
* On exit,
* if JOBU = 'I', U contains the unitary matrix U;
* if JOBU = 'U', U contains the product U1*U.
* If JOBU = 'N', U is not referenced.
*
* LDU (input) INTEGER
* The leading dimension of the array U. LDU >= max(1,M) if
* JOBU = 'U'; LDU >= 1 otherwise.
*
* V (input/output) COMPLEX array, dimension (LDV,P)
* On entry, if JOBV = 'V', V must contain a matrix V1 (usually
* the unitary matrix returned by CGGSVP).
* On exit,
* if JOBV = 'I', V contains the unitary matrix V;
* if JOBV = 'V', V contains the product V1*V.
* If JOBV = 'N', V is not referenced.
*
* LDV (input) INTEGER
* The leading dimension of the array V. LDV >= max(1,P) if
* JOBV = 'V'; LDV >= 1 otherwise.
*
* Q (input/output) COMPLEX array, dimension (LDQ,N)
* On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
* the unitary matrix returned by CGGSVP).
* On exit,
* if JOBQ = 'I', Q contains the unitary matrix Q;
* if JOBQ = 'Q', Q contains the product Q1*Q.
* If JOBQ = 'N', Q is not referenced.
*
* LDQ (input) INTEGER
* The leading dimension of the array Q. LDQ >= max(1,N) if
* JOBQ = 'Q'; LDQ >= 1 otherwise.
*
* WORK (workspace) COMPLEX array, dimension (2*N)
*
* NCYCLE (output) INTEGER
* The number of cycles required for convergence.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
* = 1: the procedure does not converge after MAXIT cycles.
*
* Internal Parameters
* ===================
*
* MAXIT INTEGER
* MAXIT specifies the total loops that the iterative procedure
* may take. If after MAXIT cycles, the routine fails to
* converge, we return INFO = 1.
*
* Further Details
* ===============
*
* CTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
* min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
* matrix B13 to the form:
*
* U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1,
*
* where U1, V1 and Q1 are unitary matrix, and Z' is the conjugate
* transpose of Z. C1 and S1 are diagonal matrices satisfying
*
* C1**2 + S1**2 = I,
*
* and R1 is an L-by-L nonsingular upper triangular matrix.
*
* =====================================================================
*
go to the page top
ctgsna
USAGE:
s, dif, m, work, info = NumRu::Lapack.ctgsna( job, howmny, select, a, b, vl, vr, [:lwork => lwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO )
* Purpose
* =======
*
* CTGSNA estimates reciprocal condition numbers for specified
* eigenvalues and/or eigenvectors of a matrix pair (A, B).
*
* (A, B) must be in generalized Schur canonical form, that is, A and
* B are both upper triangular.
*
* Arguments
* =========
*
* JOB (input) CHARACTER*1
* Specifies whether condition numbers are required for
* eigenvalues (S) or eigenvectors (DIF):
* = 'E': for eigenvalues only (S);
* = 'V': for eigenvectors only (DIF);
* = 'B': for both eigenvalues and eigenvectors (S and DIF).
*
* HOWMNY (input) CHARACTER*1
* = 'A': compute condition numbers for all eigenpairs;
* = 'S': compute condition numbers for selected eigenpairs
* specified by the array SELECT.
*
* SELECT (input) LOGICAL array, dimension (N)
* If HOWMNY = 'S', SELECT specifies the eigenpairs for which
* condition numbers are required. To select condition numbers
* for the corresponding j-th eigenvalue and/or eigenvector,
* SELECT(j) must be set to .TRUE..
* If HOWMNY = 'A', SELECT is not referenced.
*
* N (input) INTEGER
* The order of the square matrix pair (A, B). N >= 0.
*
* A (input) COMPLEX array, dimension (LDA,N)
* The upper triangular matrix A in the pair (A,B).
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* B (input) COMPLEX array, dimension (LDB,N)
* The upper triangular matrix B in the pair (A, B).
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* VL (input) COMPLEX array, dimension (LDVL,M)
* IF JOB = 'E' or 'B', VL must contain left eigenvectors of
* (A, B), corresponding to the eigenpairs specified by HOWMNY
* and SELECT. The eigenvectors must be stored in consecutive
* columns of VL, as returned by CTGEVC.
* If JOB = 'V', VL is not referenced.
*
* LDVL (input) INTEGER
* The leading dimension of the array VL. LDVL >= 1; and
* If JOB = 'E' or 'B', LDVL >= N.
*
* VR (input) COMPLEX array, dimension (LDVR,M)
* IF JOB = 'E' or 'B', VR must contain right eigenvectors of
* (A, B), corresponding to the eigenpairs specified by HOWMNY
* and SELECT. The eigenvectors must be stored in consecutive
* columns of VR, as returned by CTGEVC.
* If JOB = 'V', VR is not referenced.
*
* LDVR (input) INTEGER
* The leading dimension of the array VR. LDVR >= 1;
* If JOB = 'E' or 'B', LDVR >= N.
*
* S (output) REAL array, dimension (MM)
* If JOB = 'E' or 'B', the reciprocal condition numbers of the
* selected eigenvalues, stored in consecutive elements of the
* array.
* If JOB = 'V', S is not referenced.
*
* DIF (output) REAL array, dimension (MM)
* If JOB = 'V' or 'B', the estimated reciprocal condition
* numbers of the selected eigenvectors, stored in consecutive
* elements of the array.
* If the eigenvalues cannot be reordered to compute DIF(j),
* DIF(j) is set to 0; this can only occur when the true value
* would be very small anyway.
* For each eigenvalue/vector specified by SELECT, DIF stores
* a Frobenius norm-based estimate of Difl.
* If JOB = 'E', DIF is not referenced.
*
* MM (input) INTEGER
* The number of elements in the arrays S and DIF. MM >= M.
*
* M (output) INTEGER
* The number of elements of the arrays S and DIF used to store
* the specified condition numbers; for each selected eigenvalue
* one element is used. If HOWMNY = 'A', M is set to N.
*
* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK. LWORK >= max(1,N).
* If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
*
* IWORK (workspace) INTEGER array, dimension (N+2)
* If JOB = 'E', IWORK is not referenced.
*
* INFO (output) INTEGER
* = 0: Successful exit
* < 0: If INFO = -i, the i-th argument had an illegal value
*
* Further Details
* ===============
*
* The reciprocal of the condition number of the i-th generalized
* eigenvalue w = (a, b) is defined as
*
* S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v))
*
* where u and v are the right and left eigenvectors of (A, B)
* corresponding to w; |z| denotes the absolute value of the complex
* number, and norm(u) denotes the 2-norm of the vector u. The pair
* (a, b) corresponds to an eigenvalue w = a/b (= v'Au/v'Bu) of the
* matrix pair (A, B). If both a and b equal zero, then (A,B) is
* singular and S(I) = -1 is returned.
*
* An approximate error bound on the chordal distance between the i-th
* computed generalized eigenvalue w and the corresponding exact
* eigenvalue lambda is
*
* chord(w, lambda) <= EPS * norm(A, B) / S(I),
*
* where EPS is the machine precision.
*
* The reciprocal of the condition number of the right eigenvector u
* and left eigenvector v corresponding to the generalized eigenvalue w
* is defined as follows. Suppose
*
* (A, B) = ( a * ) ( b * ) 1
* ( 0 A22 ),( 0 B22 ) n-1
* 1 n-1 1 n-1
*
* Then the reciprocal condition number DIF(I) is
*
* Difl[(a, b), (A22, B22)] = sigma-min( Zl )
*
* where sigma-min(Zl) denotes the smallest singular value of
*
* Zl = [ kron(a, In-1) -kron(1, A22) ]
* [ kron(b, In-1) -kron(1, B22) ].
*
* Here In-1 is the identity matrix of size n-1 and X' is the conjugate
* transpose of X. kron(X, Y) is the Kronecker product between the
* matrices X and Y.
*
* We approximate the smallest singular value of Zl with an upper
* bound. This is done by CLATDF.
*
* An approximate error bound for a computed eigenvector VL(i) or
* VR(i) is given by
*
* EPS * norm(A, B) / DIF(i).
*
* See ref. [2-3] for more details and further references.
*
* Based on contributions by
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
* Umea University, S-901 87 Umea, Sweden.
*
* References
* ==========
*
* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
* M.S. Moonen et al (eds), Linear Algebra for Large Scale and
* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
* Eigenvalues of a Regular Matrix Pair (A, B) and Condition
* Estimation: Theory, Algorithms and Software, Report
* UMINF - 94.04, Department of Computing Science, Umea University,
* S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
* To appear in Numerical Algorithms, 1996.
*
* [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
* for Solving the Generalized Sylvester Equation and Estimating the
* Separation between Regular Matrix Pairs, Report UMINF - 93.23,
* Department of Computing Science, Umea University, S-901 87 Umea,
* Sweden, December 1993, Revised April 1994, Also as LAPACK Working
* Note 75.
* To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
*
* =====================================================================
*
go to the page top
ctgsy2
USAGE:
scale, info, c, f, rdsum, rdscal = NumRu::Lapack.ctgsy2( trans, ijob, a, b, c, d, e, f, rdsum, rdscal, [:usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, INFO )
* Purpose
* =======
*
* CTGSY2 solves the generalized Sylvester equation
*
* A * R - L * B = scale * C (1)
* D * R - L * E = scale * F
*
* using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
* (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
* N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
* (i.e., (A,D) and (B,E) in generalized Schur form).
*
* The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
* scaling factor chosen to avoid overflow.
*
* In matrix notation solving equation (1) corresponds to solve
* Zx = scale * b, where Z is defined as
*
* Z = [ kron(In, A) -kron(B', Im) ] (2)
* [ kron(In, D) -kron(E', Im) ],
*
* Ik is the identity matrix of size k and X' is the transpose of X.
* kron(X, Y) is the Kronecker product between the matrices X and Y.
*
* If TRANS = 'C', y in the conjugate transposed system Z'y = scale*b
* is solved for, which is equivalent to solve for R and L in
*
* A' * R + D' * L = scale * C (3)
* R * B' + L * E' = scale * -F
*
* This case is used to compute an estimate of Dif[(A, D), (B, E)] =
* = sigma_min(Z) using reverse communicaton with CLACON.
*
* CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL
* of an upper bound on the separation between to matrix pairs. Then
* the input (A, D), (B, E) are sub-pencils of two matrix pairs in
* CTGSYL.
*
* Arguments
* =========
*
* TRANS (input) CHARACTER*1
* = 'N', solve the generalized Sylvester equation (1).
* = 'T': solve the 'transposed' system (3).
*
* IJOB (input) INTEGER
* Specifies what kind of functionality to be performed.
* =0: solve (1) only.
* =1: A contribution from this subsystem to a Frobenius
* norm-based estimate of the separation between two matrix
* pairs is computed. (look ahead strategy is used).
* =2: A contribution from this subsystem to a Frobenius
* norm-based estimate of the separation between two matrix
* pairs is computed. (SGECON on sub-systems is used.)
* Not referenced if TRANS = 'T'.
*
* M (input) INTEGER
* On entry, M specifies the order of A and D, and the row
* dimension of C, F, R and L.
*
* N (input) INTEGER
* On entry, N specifies the order of B and E, and the column
* dimension of C, F, R and L.
*
* A (input) COMPLEX array, dimension (LDA, M)
* On entry, A contains an upper triangular matrix.
*
* LDA (input) INTEGER
* The leading dimension of the matrix A. LDA >= max(1, M).
*
* B (input) COMPLEX array, dimension (LDB, N)
* On entry, B contains an upper triangular matrix.
*
* LDB (input) INTEGER
* The leading dimension of the matrix B. LDB >= max(1, N).
*
* C (input/output) COMPLEX array, dimension (LDC, N)
* On entry, C contains the right-hand-side of the first matrix
* equation in (1).
* On exit, if IJOB = 0, C has been overwritten by the solution
* R.
*
* LDC (input) INTEGER
* The leading dimension of the matrix C. LDC >= max(1, M).
*
* D (input) COMPLEX array, dimension (LDD, M)
* On entry, D contains an upper triangular matrix.
*
* LDD (input) INTEGER
* The leading dimension of the matrix D. LDD >= max(1, M).
*
* E (input) COMPLEX array, dimension (LDE, N)
* On entry, E contains an upper triangular matrix.
*
* LDE (input) INTEGER
* The leading dimension of the matrix E. LDE >= max(1, N).
*
* F (input/output) COMPLEX array, dimension (LDF, N)
* On entry, F contains the right-hand-side of the second matrix
* equation in (1).
* On exit, if IJOB = 0, F has been overwritten by the solution
* L.
*
* LDF (input) INTEGER
* The leading dimension of the matrix F. LDF >= max(1, M).
*
* SCALE (output) REAL
* On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
* R and L (C and F on entry) will hold the solutions to a
* slightly perturbed system but the input matrices A, B, D and
* E have not been changed. If SCALE = 0, R and L will hold the
* solutions to the homogeneous system with C = F = 0.
* Normally, SCALE = 1.
*
* RDSUM (input/output) REAL
* On entry, the sum of squares of computed contributions to
* the Dif-estimate under computation by CTGSYL, where the
* scaling factor RDSCAL (see below) has been factored out.
* On exit, the corresponding sum of squares updated with the
* contributions from the current sub-system.
* If TRANS = 'T' RDSUM is not touched.
* NOTE: RDSUM only makes sense when CTGSY2 is called by
* CTGSYL.
*
* RDSCAL (input/output) REAL
* On entry, scaling factor used to prevent overflow in RDSUM.
* On exit, RDSCAL is updated w.r.t. the current contributions
* in RDSUM.
* If TRANS = 'T', RDSCAL is not touched.
* NOTE: RDSCAL only makes sense when CTGSY2 is called by
* CTGSYL.
*
* INFO (output) INTEGER
* On exit, if INFO is set to
* =0: Successful exit
* <0: If INFO = -i, input argument number i is illegal.
* >0: The matrix pairs (A, D) and (B, E) have common or very
* close eigenvalues.
*
* Further Details
* ===============
*
* Based on contributions by
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
* Umea University, S-901 87 Umea, Sweden.
*
* =====================================================================
*
go to the page top
ctgsyl
USAGE:
scale, dif, work, info, c, f = NumRu::Lapack.ctgsyl( trans, ijob, a, b, c, d, e, f, [:lwork => lwork, :usage => usage, :help => help])
FORTRAN MANUAL
SUBROUTINE CTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO )
* Purpose
* =======
*
* CTGSYL solves the generalized Sylvester equation:
*
* A * R - L * B = scale * C (1)
* D * R - L * E = scale * F
*
* where R and L are unknown m-by-n matrices, (A, D), (B, E) and
* (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
* respectively, with complex entries. A, B, D and E are upper
* triangular (i.e., (A,D) and (B,E) in generalized Schur form).
*
* The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
* is an output scaling factor chosen to avoid overflow.
*
* In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
* is defined as
*
* Z = [ kron(In, A) -kron(B', Im) ] (2)
* [ kron(In, D) -kron(E', Im) ],
*
* Here Ix is the identity matrix of size x and X' is the conjugate
* transpose of X. Kron(X, Y) is the Kronecker product between the
* matrices X and Y.
*
* If TRANS = 'C', y in the conjugate transposed system Z'*y = scale*b
* is solved for, which is equivalent to solve for R and L in
*
* A' * R + D' * L = scale * C (3)
* R * B' + L * E' = scale * -F
*
* This case (TRANS = 'C') is used to compute an one-norm-based estimate
* of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
* and (B,E), using CLACON.
*
* If IJOB >= 1, CTGSYL computes a Frobenius norm-based estimate of
* Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
* reciprocal of the smallest singular value of Z.
*
* This is a level-3 BLAS algorithm.
*
* Arguments
* =========
*
* TRANS (input) CHARACTER*1
* = 'N': solve the generalized sylvester equation (1).
* = 'C': solve the "conjugate transposed" system (3).
*
* IJOB (input) INTEGER
* Specifies what kind of functionality to be performed.
* =0: solve (1) only.
* =1: The functionality of 0 and 3.
* =2: The functionality of 0 and 4.
* =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
* (look ahead strategy is used).
* =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
* (CGECON on sub-systems is used).
* Not referenced if TRANS = 'C'.
*
* M (input) INTEGER
* The order of the matrices A and D, and the row dimension of
* the matrices C, F, R and L.
*
* N (input) INTEGER
* The order of the matrices B and E, and the column dimension
* of the matrices C, F, R and L.
*
* A (input) COMPLEX array, dimension (LDA, M)
* The upper triangular matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1, M).
*
* B (input) COMPLEX array, dimension (LDB, N)
* The upper triangular matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1, N).
*
* C (input/output) COMPLEX array, dimension (LDC, N)
* On entry, C contains the right-hand-side of the first matrix
* equation in (1) or (3).
* On exit, if IJOB = 0, 1 or 2, C has been overwritten by
* the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
* the solution achieved during the computation of the
* Dif-estimate.
*
* LDC (input) INTEGER
* The leading dimension of the array C. LDC >= max(1, M).
*
* D (input) COMPLEX array, dimension (LDD, M)
* The upper triangular matrix D.
*
* LDD (input) INTEGER
* The leading dimension of the array D. LDD >= max(1, M).
*
* E (input) COMPLEX array, dimension (LDE, N)
* The upper triangular matrix E.
*
* LDE (input) INTEGER
* The leading dimension of the array E. LDE >= max(1, N).
*
* F (input/output) COMPLEX array, dimension (LDF, N)
* On entry, F contains the right-hand-side of the second matrix
* equation in (1) or (3).
* On exit, if IJOB = 0, 1 or 2, F has been overwritten by
* the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
* the solution achieved during the computation of the
* Dif-estimate.
*
* LDF (input) INTEGER
* The leading dimension of the array F. LDF >= max(1, M).
*
* DIF (output) REAL
* On exit DIF is the reciprocal of a lower bound of the
* reciprocal of the Dif-function, i.e. DIF is an upper bound of
* Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2).
* IF IJOB = 0 or TRANS = 'C', DIF is not referenced.
*
* SCALE (output) REAL
* On exit SCALE is the scaling factor in (1) or (3).
* If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
* to a slightly perturbed system but the input matrices A, B,
* D and E have not been changed. If SCALE = 0, R and L will
* hold the solutions to the homogenious system with C = F = 0.
*
* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK. LWORK > = 1.
* If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* IWORK (workspace) INTEGER array, dimension (M+N+2)
*
* INFO (output) INTEGER
* =0: successful exit
* <0: If INFO = -i, the i-th argument had an illegal value.
* >0: (A, D) and (B, E) have common or very close
* eigenvalues.
*
* Further Details
* ===============
*
* Based on contributions by
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
* Umea University, S-901 87 Umea, Sweden.
*
* [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
* for Solving the Generalized Sylvester Equation and Estimating the
* Separation between Regular Matrix Pairs, Report UMINF - 93.23,
* Department of Computing Science, Umea University, S-901 87 Umea,
* Sweden, December 1993, Revised April 1994, Also as LAPACK Working
* Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
* No 1, 1996.
*
* [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
* Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
* Appl., 15(4):1045-1060, 1994.
*
* [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
* Condition Estimators for Solving the Generalized Sylvester
* Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
* July 1989, pp 745-751.
*
* =====================================================================
* Replaced various illegal calls to CCOPY by calls to CLASET.
* Sven Hammarling, 1/5/02.
*
go to the page top
back to matrix types
back to data types