USAGE: m, info, vl, vr = NumRu::Lapack.dtgevc( side, howmny, select, s, p, vl, vr, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO ) * Purpose * ======= * * DTGEVC computes some or all of the right and/or left eigenvectors of * a pair of real matrices (S,P), where S is a quasi-triangular matrix * and P is upper triangular. Matrix pairs of this type are produced by * the generalized Schur factorization of a matrix pair (A,B): * * A = Q*S*Z**T, B = Q*P*Z**T * * as computed by DGGHRD + DHGEQZ. * * 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 blocks 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 orthogonal 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. If w(j) is a real eigenvalue, the corresponding * real eigenvector is computed if SELECT(j) is .TRUE.. * If w(j) and w(j+1) are the real and imaginary parts of a * complex eigenvalue, the corresponding complex eigenvector * is computed if either SELECT(j) or SELECT(j+1) is .TRUE., * and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is * set to .FALSE.. * Not referenced if HOWMNY = 'A' or 'B'. * * N (input) INTEGER * The order of the matrices S and P. N >= 0. * * S (input) DOUBLE PRECISION array, dimension (LDS,N) * The upper quasi-triangular matrix S from a generalized Schur * factorization, as computed by DHGEQZ. * * LDS (input) INTEGER * The leading dimension of array S. LDS >= max(1,N). * * P (input) DOUBLE PRECISION array, dimension (LDP,N) * The upper triangular matrix P from a generalized Schur * factorization, as computed by DHGEQZ. * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks * of S must be in positive diagonal form. * * LDP (input) INTEGER * The leading dimension of array P. LDP >= max(1,N). * * VL (input/output) DOUBLE PRECISION 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 orthogonal matrix Q * of left Schur vectors returned by DHGEQZ). * 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. * * A complex eigenvector corresponding to a complex eigenvalue * is stored in two consecutive columns, the first holding the * real part, and the second the imaginary part. * * Not referenced if SIDE = 'R'. * * LDVL (input) INTEGER * The leading dimension of array VL. LDVL >= 1, and if * SIDE = 'L' or 'B', LDVL >= N. * * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must * contain an N-by-N matrix Z (usually the orthogonal matrix Z * of right Schur vectors returned by DHGEQZ). * * On exit, if SIDE = 'R' or 'B', VR contains: * if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); * if HOWMNY = 'B' or 'b', the matrix Z*X; * if HOWMNY = 'S' or 's', the right eigenvectors of (S,P) * specified by SELECT, stored consecutively in the * columns of VR, in the same order as their * eigenvalues. * * A complex eigenvector corresponding to a complex eigenvalue * is stored in two consecutive columns, the first holding the * real part and the second the imaginary part. * * 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 real eigenvector occupies one * column and each selected complex eigenvector occupies two * columns. * * WORK (workspace) DOUBLE PRECISION array, dimension (6*N) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex * eigenvalue. * * Further Details * =============== * * Allocation of workspace: * ---------- -- --------- * * WORK( j ) = 1-norm of j-th column of A, above the diagonal * WORK( N+j ) = 1-norm of j-th column of B, above the diagonal * WORK( 2*N+1:3*N ) = real part of eigenvector * WORK( 3*N+1:4*N ) = imaginary part of eigenvector * WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector * WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector * * Rowwise vs. columnwise solution methods: * ------- -- ---------- -------- ------- * * Finding a generalized eigenvector consists basically of solving the * singular triangular system * * (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left) * * Consider finding the i-th right eigenvector (assume all eigenvalues * are real). The equation to be solved is: * n i * 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1 * k=j k=j * * where C = (A - w B) (The components v(i+1:n) are 0.) * * The "rowwise" method is: * * (1) v(i) := 1 * for j = i-1,. . .,1: * i * (2) compute s = - sum C(j,k) v(k) and * k=j+1 * * (3) v(j) := s / C(j,j) * * Step 2 is sometimes called the "dot product" step, since it is an * inner product between the j-th row and the portion of the eigenvector * that has been computed so far. * * The "columnwise" method consists basically in doing the sums * for all the rows in parallel. As each v(j) is computed, the * contribution of v(j) times the j-th column of C is added to the * partial sums. Since FORTRAN arrays are stored columnwise, this has * the advantage that at each step, the elements of C that are accessed * are adjacent to one another, whereas with the rowwise method, the * elements accessed at a step are spaced LDS (and LDP) words apart. * * When finding left eigenvectors, the matrix in question is the * transpose of the one in storage, so the rowwise method then * actually accesses columns of A and B at each step, and so is the * preferred method. * * ===================================================================== *go to the page top

USAGE: info, a, b, q, z = NumRu::Lapack.dtgex2( wantq, wantz, a, b, q, z, j1, n1, n2, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, N1, N2, WORK, LWORK, INFO ) * Purpose * ======= * * DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22) * of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair * (A, B) by an orthogonal equivalence transformation. * * (A, B) must be in generalized real Schur canonical form (as returned * by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 * diagonal blocks. B is 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) DOUBLE PRECISION array, 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) DOUBLE PRECISION array, 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) DOUBLE PRECISION array, dimension (LDQ,N) * On entry, if WANTQ = .TRUE., the orthogonal 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) DOUBLE PRECISION array, dimension (LDZ,N) * On entry, if WANTZ =.TRUE., the orthogonal 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). 1 <= J1 <= N. * * N1 (input) INTEGER * The order of the first block (A11, B11). N1 = 0, 1 or 2. * * N2 (input) INTEGER * The order of the second block (A22, B22). N2 = 0, 1 or 2. * * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)). * * LWORK (input) INTEGER * The dimension of the array WORK. * LWORK >= MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 ) * * INFO (output) INTEGER * =0: Successful exit * >0: If INFO = 1, the transformed matrix (A, B) would be * too far from generalized Schur form; the blocks are * not swapped and (A, B) and (Q, Z) are unchanged. * The problem of swapping is too ill-conditioned. * <0: If INFO = -16: LWORK is too small. Appropriate value * for LWORK is returned in WORK(1). * * 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. * * ===================================================================== * Replaced various illegal calls to DCOPY by calls to DLASET, or by DO * loops. Sven Hammarling, 1/5/02. *go to the page top

USAGE: work, info, a, b, q, z, ifst, ilst = NumRu::Lapack.dtgexc( wantq, wantz, a, b, q, z, ifst, ilst, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO ) * Purpose * ======= * * DTGEXC reorders the generalized real Schur decomposition of a real * matrix pair (A,B) using an orthogonal 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 real Schur canonical form (as returned * by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 * diagonal blocks. B is 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) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the matrix A in generalized real Schur canonical * form. * On exit, the updated matrix A, again in generalized * real Schur canonical form. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,N) * On entry, the matrix B in generalized real Schur canonical * form (A,B). * On exit, the updated matrix B, again in generalized * real Schur canonical form (A,B). * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) * On entry, if WANTQ = .TRUE., the orthogonal 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) DOUBLE PRECISION array, dimension (LDZ,N) * On entry, if WANTZ = .TRUE., the orthogonal 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/output) 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. * On exit, if IFST pointed on entry to the second row of * a 2-by-2 block, it is changed to point to the first row; * ILST always points to the first row of the block in its * final position (which may differ from its input value by * +1 or -1). 1 <= IFST, ILST <= N. * * WORK (workspace/output) DOUBLE PRECISION 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 when N <= 1, otherwise LWORK >= 4*N + 16. * * 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. * * 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. * * ===================================================================== *go to the page top

USAGE: alphar, alphai, beta, m, pl, pr, dif, work, iwork, info, a, b, q, z = NumRu::Lapack.dtgsen( ijob, wantq, wantz, select, a, b, q, z, [:lwork => lwork, :liwork => liwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO ) * Purpose * ======= * * DTGSEN reorders the generalized real Schur decomposition of a real * matrix pair (A, B) (in terms of an orthonormal equivalence trans- * formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues * appears in the leading diagonal blocks of the upper quasi-triangular * matrix A and the upper triangular B. The leading columns of Q and * Z form orthonormal bases of the corresponding left and right eigen- * spaces (deflating subspaces). (A, B) must be in generalized real * Schur canonical form (as returned by DGGES), i.e. A is block upper * triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper * triangular. * * DTGSEN also computes the generalized eigenvalues * * w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) * * of the reordered matrix pair (A, B). * * Optionally, DTGSEN computes the 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 a real eigenvalue w(j), SELECT(j) must be set to * .TRUE.. To select a complex conjugate pair of eigenvalues * w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, * either SELECT(j) or SELECT(j+1) or both must be set to * .TRUE.; a complex conjugate pair of eigenvalues must be * either both included in the cluster or both excluded. * * N (input) INTEGER * The order of the matrices A and B. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension(LDA,N) * On entry, the upper quasi-triangular matrix A, with (A, B) in * generalized real 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) DOUBLE PRECISION array, dimension(LDB,N) * On entry, the upper triangular matrix B, with (A, B) in * generalized real 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). * * ALPHAR (output) DOUBLE PRECISION array, dimension (N) * ALPHAI (output) DOUBLE PRECISION array, dimension (N) * BETA (output) DOUBLE PRECISION array, dimension (N) * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will * be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i * and BETA(j),j=1,...,N are the diagonals of the complex Schur * form (S,T) that would result if the 2-by-2 diagonal blocks of * the real generalized Schur form of (A,B) were further reduced * to triangular form using complex unitary transformations. * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if * positive, then the j-th and (j+1)-st eigenvalues are a * complex conjugate pair, with ALPHAI(j+1) negative. * * Q (input/output) DOUBLE PRECISION 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 orthogonal * 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; * and if WANTQ = .TRUE., LDQ >= N. * * Z (input/output) DOUBLE PRECISION 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 orthogonal * 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 eigen- * spaces (deflating subspaces). 0 <= M <= N. * * PL (output) DOUBLE PRECISION * PR (output) DOUBLE PRECISION * If IJOB = 1, 4 or 5, PL, PR are lower bounds on the * reciprocal of the norm of "projections" onto left and right * eigenspaces 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 and PR are not referenced. * * DIF (output) DOUBLE PRECISION 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. * If M = 0 or N, DIF(1:2) = F-norm([A, B]). * If IJOB = 0 or 1, DIF is not referenced. * * WORK (workspace/output) DOUBLE PRECISION 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 >= 4*N+16. * If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)). * If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 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+6. * If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6). * * 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 * =============== * * DTGSEN first collects the selected eigenvalues by computing * orthogonal 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 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 real 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 DLATDF), then the parameter * IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF * (IJOB = 2 will be used)). See DTGSYL 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.dtgsja( jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGSJA( 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 * ======= * * DTGSJA computes the generalized singular value decomposition (GSVD) * of two real 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 DGGSVP * 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 orthogonal matrices, Z' denotes the 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 orthogonal 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 an orthogonal matrix U1 on entry, and * the product U1*U is returned; * = 'I': U is initialized to the unit matrix, and the * orthogonal matrix U is returned; * = 'N': U is not computed. * * JOBV (input) CHARACTER*1 * = 'V': V must contain an orthogonal matrix V1 on entry, and * the product V1*V is returned; * = 'I': V is initialized to the unit matrix, and the * orthogonal matrix V is returned; * = 'N': V is not computed. * * JOBQ (input) CHARACTER*1 * = 'Q': Q must contain an orthogonal matrix Q1 on entry, and * the product Q1*Q is returned; * = 'I': Q is initialized to the unit matrix, and the * orthogonal 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 DTGSJA. * See Further Details. * * A (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION * TOLB (input) DOUBLE PRECISION * 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)*MAZHEPS, * TOLB = max(P,N)*norm(B)*MAZHEPS. * * ALPHA (output) DOUBLE PRECISION array, dimension (N) * BETA (output) DOUBLE PRECISION 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 and * BETA(K+L+1:N) = 0. * * U (input/output) DOUBLE PRECISION array, dimension (LDU,M) * On entry, if JOBU = 'U', U must contain a matrix U1 (usually * the orthogonal matrix returned by DGGSVP). * On exit, * if JOBU = 'I', U contains the orthogonal 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) DOUBLE PRECISION array, dimension (LDV,P) * On entry, if JOBV = 'V', V must contain a matrix V1 (usually * the orthogonal matrix returned by DGGSVP). * On exit, * if JOBV = 'I', V contains the orthogonal 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) DOUBLE PRECISION array, dimension (LDQ,N) * On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually * the orthogonal matrix returned by DGGSVP). * On exit, * if JOBQ = 'I', Q contains the orthogonal 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) DOUBLE PRECISION 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 * =============== * * DTGSJA 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 orthogonal matrix, and Z' is the 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.dtgsna( job, howmny, select, a, b, vl, vr, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO ) * Purpose * ======= * * DTGSNA estimates reciprocal condition numbers for specified * eigenvalues and/or eigenvectors of a matrix pair (A, B) in * generalized real Schur canonical form (or of any matrix pair * (Q*A*Z', Q*B*Z') with orthogonal matrices Q and Z, where * Z' denotes the transpose of Z. * * (A, B) must be in generalized real Schur form (as returned by DGGES), * i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal * blocks. B is 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 eigenpair corresponding to a real eigenvalue w(j), * SELECT(j) must be set to .TRUE.. To select condition numbers * corresponding to a complex conjugate pair of eigenvalues w(j) * and w(j+1), either SELECT(j) or SELECT(j+1) or both, 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) DOUBLE PRECISION array, dimension (LDA,N) * The upper quasi-triangular matrix A in the pair (A,B). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * B (input) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DTGEVC. * If JOB = 'V', VL is not referenced. * * LDVL (input) INTEGER * The leading dimension of the array VL. LDVL >= 1. * If JOB = 'E' or 'B', LDVL >= N. * * VR (input) DOUBLE PRECISION 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 ov VR, as returned by DTGEVC. * 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) DOUBLE PRECISION array, dimension (MM) * If JOB = 'E' or 'B', the reciprocal condition numbers of the * selected eigenvalues, stored in consecutive elements of the * array. For a complex conjugate pair of eigenvalues two * consecutive elements of S are set to the same value. Thus * S(j), DIF(j), and the j-th columns of VL and VR all * correspond to the same eigenpair (but not in general the * j-th eigenpair, unless all eigenpairs are selected). * If JOB = 'V', S is not referenced. * * DIF (output) DOUBLE PRECISION array, dimension (MM) * If JOB = 'V' or 'B', the estimated reciprocal condition * numbers of the selected eigenvectors, stored in consecutive * elements of the array. For a complex eigenvector two * consecutive elements of DIF are set to the same value. 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. * 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 real * eigenvalue one element is used, and for each selected complex * conjugate pair of eigenvalues, two elements are used. * If HOWMNY = 'A', M is set to N. * * WORK (workspace/output) DOUBLE PRECISION 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 >= 2*N*(N+2)+16. * * 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 (N + 6) * 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 a generalized eigenvalue * w = (a, b) is defined as * * S(w) = (|u'Av|**2 + |u'Bv|**2)**(1/2) / (norm(u)*norm(v)) * * where u and v are the left and right 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 (= u'Av/u'Bv) * 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 DIF(i) of right eigenvector u * and left eigenvector v corresponding to the generalized eigenvalue w * is defined as follows: * * a) If the i-th eigenvalue w = (a,b) is real * * Suppose U and V are orthogonal transformations such that * * U'*(A, B)*V = (S, T) = ( a * ) ( b * ) 1 * ( 0 S22 ),( 0 T22 ) n-1 * 1 n-1 1 n-1 * * Then the reciprocal condition number DIF(i) is * * Difl((a, b), (S22, T22)) = sigma-min( Zl ), * * where sigma-min(Zl) denotes the smallest singular value of the * 2(n-1)-by-2(n-1) matrix * * Zl = [ kron(a, In-1) -kron(1, S22) ] * [ kron(b, In-1) -kron(1, T22) ] . * * Here In-1 is the identity matrix of size n-1. kron(X, Y) is the * Kronecker product between the matrices X and Y. * * Note that if the default method for computing DIF(i) is wanted * (see DLATDF), then the parameter DIFDRI (see below) should be * changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). * See DTGSYL for more details. * * b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair, * * Suppose U and V are orthogonal transformations such that * * U'*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2 * ( 0 S22 ),( 0 T22) n-2 * 2 n-2 2 n-2 * * and (S11, T11) corresponds to the complex conjugate eigenvalue * pair (w, conjg(w)). There exist unitary matrices U1 and V1 such * that * * U1'*S11*V1 = ( s11 s12 ) and U1'*T11*V1 = ( t11 t12 ) * ( 0 s22 ) ( 0 t22 ) * * where the generalized eigenvalues w = s11/t11 and * conjg(w) = s22/t22. * * Then the reciprocal condition number DIF(i) is bounded by * * min( d1, max( 1, |real(s11)/real(s22)| )*d2 ) * * where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where * Z1 is the complex 2-by-2 matrix * * Z1 = [ s11 -s22 ] * [ t11 -t22 ], * * This is done by computing (using real arithmetic) the * roots of the characteristical polynomial det(Z1' * Z1 - lambda I), * where Z1' denotes the conjugate transpose of Z1 and det(X) denotes * the determinant of X. * * and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an * upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2) * * Z2 = [ kron(S11', In-2) -kron(I2, S22) ] * [ kron(T11', In-2) -kron(I2, T22) ] * * Note that if the default method for computing DIF is wanted (see * DLATDF), then the parameter DIFDRI (see below) should be changed * from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL * for more details. * * For each eigenvalue/vector specified by SELECT, DIF stores a * Frobenius norm-based estimate of Difl. * * An approximate error bound for the i-th 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, pq, info, c, f, rdsum, rdscal = NumRu::Lapack.dtgsy2( trans, ijob, a, b, c, d, e, f, rdsum, rdscal, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO ) * Purpose * ======= * * DTGSY2 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, with real entries. (A, D) and (B, E) * must be in generalized Schur canonical form, i.e. A, B are upper * quasi triangular and D, E are upper triangular. 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 * Z*x = 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. * In the process of solving (1), we solve a number of such systems * where Dim(In), Dim(In) = 1 or 2. * * If TRANS = 'T', solve the transposed system Z'*y = scale*b for y, * 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 DLACON. * * DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL * of an upper bound on the separation between to matrix pairs. Then * the input (A, D), (B, E) are sub-pencils of the matrix pair in * DTGSYL. See DTGSYL for details. * * 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. (DGECON 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) DOUBLE PRECISION array, dimension (LDA, M) * On entry, A contains an upper quasi triangular matrix. * * LDA (input) INTEGER * The leading dimension of the matrix A. LDA >= max(1, M). * * B (input) DOUBLE PRECISION array, dimension (LDB, N) * On entry, B contains an upper quasi triangular matrix. * * LDB (input) INTEGER * The leading dimension of the matrix B. LDB >= max(1, N). * * C (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION * 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) DOUBLE PRECISION * On entry, the sum of squares of computed contributions to * the Dif-estimate under computation by DTGSYL, 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 DTGSY2 is called by DTGSYL. * * RDSCAL (input/output) DOUBLE PRECISION * 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 DTGSY2 is called by * DTGSYL. * * IWORK (workspace) INTEGER array, dimension (M+N+2) * * PQ (output) INTEGER * On exit, the number of subsystems (of size 2-by-2, 4-by-4 and * 8-by-8) solved by this routine. * * INFO (output) INTEGER * On exit, if INFO is set to * =0: Successful exit * <0: If INFO = -i, the i-th argument had an illegal value. * >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. * * ===================================================================== * Replaced various illegal calls to DCOPY by calls to DLASET. * Sven Hammarling, 27/5/02. *go to the page top

USAGE: scale, dif, work, info, c, f = NumRu::Lapack.dtgsyl( trans, ijob, a, b, c, d, e, f, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO ) * Purpose * ======= * * DTGSYL 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 real entries. (A, D) and (B, E) must be in * generalized (real) Schur canonical form, i.e. A, B are upper quasi * triangular and D, E are upper triangular. * * 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 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 = 'T', DTGSYL solves the transposed system Z'*y = scale*b, * 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 = 'T') 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 DLACON. * * If IJOB >= 1, DTGSYL 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. See [1-2] for more * information. * * This is a level 3 BLAS algorithm. * * 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: 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 IJOB = 1 is used). * =4: Only an estimate of Dif[(A,D), (B,E)] is computed. * ( DGECON on sub-systems is used ). * Not referenced if TRANS = 'T'. * * 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) DOUBLE PRECISION array, dimension (LDA, M) * The upper quasi triangular matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1, M). * * B (input) DOUBLE PRECISION array, dimension (LDB, N) * The upper quasi triangular matrix B. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1, N). * * C (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION * 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 = 'T', DIF is not touched. * * SCALE (output) DOUBLE PRECISION * 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, C and F hold the * solutions R and L, respectively, to the homogeneous system * with C = F = 0. Normally, SCALE = 1. * * WORK (workspace/output) DOUBLE PRECISION 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+6) * * 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 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 DCOPY by calls to DLASET. * Sven Hammarling, 1/5/02. *go to the page top

back to matrix types

back to data types