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

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

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

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

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

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

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

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