# COMPLEX routines for triangular matrices, generalized problem (i.e., a pair of triangular matrices) matrix

## ctgevc

```USAGE:
m, info, vl, vr = NumRu::Lapack.ctgevc( side, howmny, select, s, p, vl, vr, [:usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )

*  Purpose
*  =======
*
*  CTGEVC computes some or all of the right and/or left eigenvectors of
*  a pair of complex matrices (S,P), where S and P are upper triangular.
*  Matrix pairs of this type are produced by the generalized Schur
*  factorization of a complex matrix pair (A,B):
*
*     A = Q*S*Z**H,  B = Q*P*Z**H
*
*  as computed by CGGHRD + CHGEQZ.
*
*  The right eigenvector x and the left eigenvector y of (S,P)
*  corresponding to an eigenvalue w are defined by:
*
*     S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
*
*  where y**H denotes the conjugate tranpose of y.
*  The eigenvalues are not input to this routine, but are computed
*  directly from the diagonal elements of S and P.
*
*  This routine returns the matrices X and/or Y of right and left
*  eigenvectors of (S,P), or the products Z*X and/or Q*Y,
*  where Z and Q are input matrices.
*  If Q and Z are the unitary factors from the generalized Schur
*  factorization of a matrix pair (A,B), then Z*X and Q*Y
*  are the matrices of right and left eigenvectors of (A,B).
*

*  Arguments
*  =========
*
*  SIDE    (input) CHARACTER*1
*          = 'R': compute right eigenvectors only;
*          = 'L': compute left eigenvectors only;
*          = 'B': compute both right and left eigenvectors.
*
*  HOWMNY  (input) CHARACTER*1
*          = 'A': compute all right and/or left eigenvectors;
*          = 'B': compute all right and/or left eigenvectors,
*                 backtransformed by the matrices in VR and/or VL;
*          = 'S': compute selected right and/or left eigenvectors,
*                 specified by the logical array SELECT.
*
*  SELECT  (input) LOGICAL array, dimension (N)
*          If HOWMNY='S', SELECT specifies the eigenvectors to be
*          computed.  The eigenvector corresponding to the j-th
*          eigenvalue is computed if SELECT(j) = .TRUE..
*          Not referenced if HOWMNY = 'A' or 'B'.
*
*  N       (input) INTEGER
*          The order of the matrices S and P.  N >= 0.
*
*  S       (input) COMPLEX array, dimension (LDS,N)
*          The upper triangular matrix S from a generalized Schur
*          factorization, as computed by CHGEQZ.
*
*  LDS     (input) INTEGER
*          The leading dimension of array S.  LDS >= max(1,N).
*
*  P       (input) COMPLEX array, dimension (LDP,N)
*          The upper triangular matrix P from a generalized Schur
*          factorization, as computed by CHGEQZ.  P must have real
*          diagonal elements.
*
*  LDP     (input) INTEGER
*          The leading dimension of array P.  LDP >= max(1,N).
*
*  VL      (input/output) COMPLEX array, dimension (LDVL,MM)
*          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
*          contain an N-by-N matrix Q (usually the unitary matrix Q
*          of left Schur vectors returned by CHGEQZ).
*          On exit, if SIDE = 'L' or 'B', VL contains:
*          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
*          if HOWMNY = 'B', the matrix Q*Y;
*          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
*                      SELECT, stored consecutively in the columns of
*                      VL, in the same order as their eigenvalues.
*          Not referenced if SIDE = 'R'.
*
*  LDVL    (input) INTEGER
*          The leading dimension of array VL.  LDVL >= 1, and if
*          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
*
*  VR      (input/output) COMPLEX array, dimension (LDVR,MM)
*          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
*          contain an N-by-N matrix Q (usually the unitary matrix Z
*          of right Schur vectors returned by CHGEQZ).
*          On exit, if SIDE = 'R' or 'B', VR contains:
*          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
*          if HOWMNY = 'B', the matrix Z*X;
*          if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
*                      SELECT, stored consecutively in the columns of
*                      VR, in the same order as their eigenvalues.
*          Not referenced if SIDE = 'L'.
*
*  LDVR    (input) INTEGER
*          The leading dimension of the array VR.  LDVR >= 1, and if
*          SIDE = 'R' or 'B', LDVR >= N.
*
*  MM      (input) INTEGER
*          The number of columns in the arrays VL and/or VR. MM >= M.
*
*  M       (output) INTEGER
*          The number of columns in the arrays VL and/or VR actually
*          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
*          is set to N.  Each selected eigenvector occupies one column.
*
*  WORK    (workspace) COMPLEX array, dimension (2*N)
*
*  RWORK   (workspace) REAL array, dimension (2*N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit.
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*

*  =====================================================================
*

```
go to the page top

## ctgex2

```USAGE:
info, a, b, q, z = NumRu::Lapack.ctgex2( wantq, wantz, a, b, q, ldq, z, ldz, j1, [:usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, INFO )

*  Purpose
*  =======
*
*  CTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
*  in an upper triangular matrix pair (A, B) by an unitary equivalence
*  transformation.
*
*  (A, B) must be in generalized Schur canonical form, that is, A and
*  B are both upper triangular.
*
*  Optionally, the matrices Q and Z of generalized Schur vectors are
*  updated.
*
*         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
*         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
*
*

*  Arguments
*  =========
*
*  WANTQ   (input) LOGICAL
*          .TRUE. : update the left transformation matrix Q;
*          .FALSE.: do not update Q.
*
*  WANTZ   (input) LOGICAL
*          .TRUE. : update the right transformation matrix Z;
*          .FALSE.: do not update Z.
*
*  N       (input) INTEGER
*          The order of the matrices A and B. N >= 0.
*
*  A       (input/output) COMPLEX arrays, dimensions (LDA,N)
*          On entry, the matrix A in the pair (A, B).
*          On exit, the updated matrix A.
*
*  LDA     (input)  INTEGER
*          The leading dimension of the array A. LDA >= max(1,N).
*
*  B       (input/output) COMPLEX arrays, dimensions (LDB,N)
*          On entry, the matrix B in the pair (A, B).
*          On exit, the updated matrix B.
*
*  LDB     (input)  INTEGER
*          The leading dimension of the array B. LDB >= max(1,N).
*
*  Q       (input/output) COMPLEX array, dimension (LDZ,N)
*          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
*          the updated matrix Q.
*          Not referenced if WANTQ = .FALSE..
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q. LDQ >= 1;
*          If WANTQ = .TRUE., LDQ >= N.
*
*  Z       (input/output) COMPLEX array, dimension (LDZ,N)
*          If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
*          the updated matrix Z.
*          Not referenced if WANTZ = .FALSE..
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z. LDZ >= 1;
*          If WANTZ = .TRUE., LDZ >= N.
*
*  J1      (input) INTEGER
*          The index to the first block (A11, B11).
*
*  INFO    (output) INTEGER
*           =0:  Successful exit.
*           =1:  The transformed matrix pair (A, B) would be too far
*                from generalized Schur form; the problem is ill-
*                conditioned.
*
*

*  Further Details
*  ===============
*
*  Based on contributions by
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*     Umea University, S-901 87 Umea, Sweden.
*
*  In the current code both weak and strong stability tests are
*  performed. The user can omit the strong stability test by changing
*  the internal logical parameter WANDS to .FALSE.. See ref. [2] for
*  details.
*
*  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
*      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
*      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
*      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
*  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
*      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
*      Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
*      Department of Computing Science, Umea University, S-901 87 Umea,
*      Sweden, 1994. Also as LAPACK Working Note 87. To appear in
*      Numerical Algorithms, 1996.
*
*  =====================================================================
*

```
go to the page top

## ctgexc

```USAGE:
info, a, b, q, z, ilst = NumRu::Lapack.ctgexc( wantq, wantz, a, b, q, ldq, z, ifst, ilst, [:usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO )

*  Purpose
*  =======
*
*  CTGEXC reorders the generalized Schur decomposition of a complex
*  matrix pair (A,B), using an unitary equivalence transformation
*  (A, B) := Q * (A, B) * Z', so that the diagonal block of (A, B) with
*  row index IFST is moved to row ILST.
*
*  (A, B) must be in generalized Schur canonical form, that is, A and
*  B are both upper triangular.
*
*  Optionally, the matrices Q and Z of generalized Schur vectors are
*  updated.
*
*         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
*         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
*

*  Arguments
*  =========
*
*  WANTQ   (input) LOGICAL
*          .TRUE. : update the left transformation matrix Q;
*          .FALSE.: do not update Q.
*
*  WANTZ   (input) LOGICAL
*          .TRUE. : update the right transformation matrix Z;
*          .FALSE.: do not update Z.
*
*  N       (input) INTEGER
*          The order of the matrices A and B. N >= 0.
*
*  A       (input/output) COMPLEX array, dimension (LDA,N)
*          On entry, the upper triangular matrix A in the pair (A, B).
*          On exit, the updated matrix A.
*
*  LDA     (input)  INTEGER
*          The leading dimension of the array A. LDA >= max(1,N).
*
*  B       (input/output) COMPLEX array, dimension (LDB,N)
*          On entry, the upper triangular matrix B in the pair (A, B).
*          On exit, the updated matrix B.
*
*  LDB     (input)  INTEGER
*          The leading dimension of the array B. LDB >= max(1,N).
*
*  Q       (input/output) COMPLEX array, dimension (LDZ,N)
*          On entry, if WANTQ = .TRUE., the unitary matrix Q.
*          On exit, the updated matrix Q.
*          If WANTQ = .FALSE., Q is not referenced.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q. LDQ >= 1;
*          If WANTQ = .TRUE., LDQ >= N.
*
*  Z       (input/output) COMPLEX array, dimension (LDZ,N)
*          On entry, if WANTZ = .TRUE., the unitary matrix Z.
*          On exit, the updated matrix Z.
*          If WANTZ = .FALSE., Z is not referenced.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z. LDZ >= 1;
*          If WANTZ = .TRUE., LDZ >= N.
*
*  IFST    (input) INTEGER
*  ILST    (input/output) INTEGER
*          Specify the reordering of the diagonal blocks of (A, B).
*          The block with row index IFST is moved to row ILST, by a
*          sequence of swapping between adjacent blocks.
*
*  INFO    (output) INTEGER
*           =0:  Successful exit.
*           <0:  if INFO = -i, the i-th argument had an illegal value.
*           =1:  The transformed matrix pair (A, B) would be too far
*                from generalized Schur form; the problem is ill-
*                conditioned. (A, B) may have been partially reordered,
*                and ILST points to the first row of the current
*                position of the block being moved.
*
*

*  Further Details
*  ===============
*
*  Based on contributions by
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*     Umea University, S-901 87 Umea, Sweden.
*
*  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
*      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
*      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
*      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
*  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
*      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
*      Estimation: Theory, Algorithms and Software, Report
*      UMINF - 94.04, Department of Computing Science, Umea University,
*      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
*      To appear in Numerical Algorithms, 1996.
*
*  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
*      for Solving the Generalized Sylvester Equation and Estimating the
*      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
*      Department of Computing Science, Umea University, S-901 87 Umea,
*      Sweden, December 1993, Revised April 1994, Also as LAPACK working
*      Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
*      1996.
*
*  =====================================================================
*
*     .. Local Scalars ..
INTEGER            HERE
*     ..
*     .. External Subroutines ..
EXTERNAL           CTGEX2, XERBLA
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          MAX
*     ..

```
go to the page top

## ctgsen

```USAGE:
alpha, beta, m, pl, pr, dif, work, iwork, info, a, b, q, z = NumRu::Lapack.ctgsen( ijob, wantq, wantz, select, a, b, q, z, [:lwork => lwork, :liwork => liwork, :usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO )

*  Purpose
*  =======
*
*  CTGSEN reorders the generalized Schur decomposition of a complex
*  matrix pair (A, B) (in terms of an unitary equivalence trans-
*  formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues
*  appears in the leading diagonal blocks of the pair (A,B). The leading
*  columns of Q and Z form unitary bases of the corresponding left and
*  right eigenspaces (deflating subspaces). (A, B) must be in
*  generalized Schur canonical form, that is, A and B are both upper
*  triangular.
*
*  CTGSEN also computes the generalized eigenvalues
*
*           w(j)= ALPHA(j) / BETA(j)
*
*  of the reordered matrix pair (A, B).
*
*  Optionally, the routine computes estimates of reciprocal condition
*  numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
*  (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
*  between the matrix pairs (A11, B11) and (A22,B22) that correspond to
*  the selected cluster and the eigenvalues outside the cluster, resp.,
*  and norms of "projections" onto left and right eigenspaces w.r.t.
*  the selected cluster in the (1,1)-block.
*
*

*  Arguments
*  =========
*
*  IJOB    (input) integer
*          Specifies whether condition numbers are required for the
*          cluster of eigenvalues (PL and PR) or the deflating subspaces
*          (Difu and Difl):
*           =0: Only reorder w.r.t. SELECT. No extras.
*           =1: Reciprocal of norms of "projections" onto left and right
*               eigenspaces w.r.t. the selected cluster (PL and PR).
*           =2: Upper bounds on Difu and Difl. F-norm-based estimate
*               (DIF(1:2)).
*           =3: Estimate of Difu and Difl. 1-norm-based estimate
*               (DIF(1:2)).
*               About 5 times as expensive as IJOB = 2.
*           =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
*               version to get it all.
*           =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
*
*  WANTQ   (input) LOGICAL
*          .TRUE. : update the left transformation matrix Q;
*          .FALSE.: do not update Q.
*
*  WANTZ   (input) LOGICAL
*          .TRUE. : update the right transformation matrix Z;
*          .FALSE.: do not update Z.
*
*  SELECT  (input) LOGICAL array, dimension (N)
*          SELECT specifies the eigenvalues in the selected cluster. To
*          select an eigenvalue w(j), SELECT(j) must be set to
*          .TRUE..
*
*  N       (input) INTEGER
*          The order of the matrices A and B. N >= 0.
*
*  A       (input/output) COMPLEX array, dimension(LDA,N)
*          On entry, the upper triangular matrix A, in generalized
*          Schur canonical form.
*          On exit, A is overwritten by the reordered matrix A.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A. LDA >= max(1,N).
*
*  B       (input/output) COMPLEX array, dimension(LDB,N)
*          On entry, the upper triangular matrix B, in generalized
*          Schur canonical form.
*          On exit, B is overwritten by the reordered matrix B.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B. LDB >= max(1,N).
*
*  ALPHA   (output) COMPLEX array, dimension (N)
*  BETA    (output) COMPLEX array, dimension (N)
*          The diagonal elements of A and B, respectively,
*          when the pair (A,B) has been reduced to generalized Schur
*          form.  ALPHA(i)/BETA(i) i=1,...,N are the generalized
*          eigenvalues.
*
*  Q       (input/output) COMPLEX array, dimension (LDQ,N)
*          On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.
*          On exit, Q has been postmultiplied by the left unitary
*          transformation matrix which reorder (A, B); The leading M
*          columns of Q form orthonormal bases for the specified pair of
*          left eigenspaces (deflating subspaces).
*          If WANTQ = .FALSE., Q is not referenced.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q. LDQ >= 1.
*          If WANTQ = .TRUE., LDQ >= N.
*
*  Z       (input/output) COMPLEX array, dimension (LDZ,N)
*          On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.
*          On exit, Z has been postmultiplied by the left unitary
*          transformation matrix which reorder (A, B); The leading M
*          columns of Z form orthonormal bases for the specified pair of
*          left eigenspaces (deflating subspaces).
*          If WANTZ = .FALSE., Z is not referenced.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z. LDZ >= 1.
*          If WANTZ = .TRUE., LDZ >= N.
*
*  M       (output) INTEGER
*          The dimension of the specified pair of left and right
*          eigenspaces, (deflating subspaces) 0 <= M <= N.
*
*  PL	   (output) REAL
*  PR	   (output) REAL
*          If IJOB = 1, 4 or 5, PL, PR are lower bounds on the
*          reciprocal  of the norm of "projections" onto left and right
*          eigenspace with respect to the selected cluster.
*          0 < PL, PR <= 1.
*          If M = 0 or M = N, PL = PR  = 1.
*          If IJOB = 0, 2 or 3 PL, PR are not referenced.
*
*  DIF     (output) REAL array, dimension (2).
*          If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
*          If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
*          Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
*          estimates of Difu and Difl, computed using reversed
*          communication with CLACN2.
*          If M = 0 or N, DIF(1:2) = F-norm([A, B]).
*          If IJOB = 0 or 1, DIF is not referenced.
*
*  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK. LWORK >=  1
*          If IJOB = 1, 2 or 4, LWORK >=  2*M*(N-M)
*          If IJOB = 3 or 5, LWORK >=  4*M*(N-M)
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by XERBLA.
*
*  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
*          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
*
*  LIWORK  (input) INTEGER
*          The dimension of the array IWORK. LIWORK >= 1.
*          If IJOB = 1, 2 or 4, LIWORK >=  N+2;
*          If IJOB = 3 or 5, LIWORK >= MAX(N+2, 2*M*(N-M));
*
*          If LIWORK = -1, then a workspace query is assumed; the
*          routine only calculates the optimal size of the IWORK array,
*          returns this value as the first entry of the IWORK array, and
*          no error message related to LIWORK is issued by XERBLA.
*
*  INFO    (output) INTEGER
*            =0: Successful exit.
*            <0: If INFO = -i, the i-th argument had an illegal value.
*            =1: Reordering of (A, B) failed because the transformed
*                matrix pair (A, B) would be too far from generalized
*                Schur form; the problem is very ill-conditioned.
*                (A, B) may have been partially reordered.
*                If requested, 0 is returned in DIF(*), PL and PR.
*
*

*  Further Details
*  ===============
*
*  CTGSEN first collects the selected eigenvalues by computing unitary
*  U and W that move them to the top left corner of (A, B). In other
*  words, the selected eigenvalues are the eigenvalues of (A11, B11) in
*
*                U'*(A, B)*W = (A11 A12) (B11 B12) n1
*                              ( 0  A22),( 0  B22) n2
*                                n1  n2    n1  n2
*
*  where N = n1+n2 and U' means the conjugate transpose of U. The first
*  n1 columns of U and W span the specified pair of left and right
*  eigenspaces (deflating subspaces) of (A, B).
*
*  If (A, B) has been obtained from the generalized real Schur
*  decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the
*  reordered generalized Schur form of (C, D) is given by
*
*           (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
*
*  and the first n1 columns of Q*U and Z*W span the corresponding
*  deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
*
*  Note that if the selected eigenvalue is sufficiently ill-conditioned,
*  then its value may differ significantly from its value before
*  reordering.
*
*  The reciprocal condition numbers of the left and right eigenspaces
*  spanned by the first n1 columns of U and W (or Q*U and Z*W) may
*  be returned in DIF(1:2), corresponding to Difu and Difl, resp.
*
*  The Difu and Difl are defined as:
*
*       Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
*  and
*       Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
*
*  where sigma-min(Zu) is the smallest singular value of the
*  (2*n1*n2)-by-(2*n1*n2) matrix
*
*       Zu = [ kron(In2, A11)  -kron(A22', In1) ]
*            [ kron(In2, B11)  -kron(B22', In1) ].
*
*  Here, Inx is the identity matrix of size nx and A22' is the
*  transpose of A22. kron(X, Y) is the Kronecker product between
*  the matrices X and Y.
*
*  When DIF(2) is small, small changes in (A, B) can cause large changes
*  in the deflating subspace. An approximate (asymptotic) bound on the
*  maximum angular error in the computed deflating subspaces is
*
*       EPS * norm((A, B)) / DIF(2),
*
*  where EPS is the machine precision.
*
*  The reciprocal norm of the projectors on the left and right
*  eigenspaces associated with (A11, B11) may be returned in PL and PR.
*  They are computed as follows. First we compute L and R so that
*  P*(A, B)*Q is block diagonal, where
*
*       P = ( I -L ) n1           Q = ( I R ) n1
*           ( 0  I ) n2    and        ( 0 I ) n2
*             n1 n2                    n1 n2
*
*  and (L, R) is the solution to the generalized Sylvester equation
*
*       A11*R - L*A22 = -A12
*       B11*R - L*B22 = -B12
*
*  Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
*  An approximate (asymptotic) bound on the average absolute error of
*  the selected eigenvalues is
*
*       EPS * norm((A, B)) / PL.
*
*  There are also global error bounds which valid for perturbations up
*  to a certain restriction:  A lower bound (x) on the smallest
*  F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
*  coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
*  (i.e. (A + E, B + F), is
*
*   x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
*
*  An approximate bound on x can be computed from DIF(1:2), PL and PR.
*
*  If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
*  (L', R') and unperturbed (L, R) left and right deflating subspaces
*  associated with the selected cluster in the (1,1)-blocks can be
*  bounded as
*
*   max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
*   max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
*
*  See LAPACK User's Guide section 4.11 or the following references
*  for more information.
*
*  Note that if the default method for computing the Frobenius-norm-
*  based estimate DIF is not wanted (see CLATDF), then the parameter
*  IDIFJB (see below) should be changed from 3 to 4 (routine CLATDF
*  (IJOB = 2 will be used)). See CTGSYL for more details.
*
*  Based on contributions by
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*     Umea University, S-901 87 Umea, Sweden.
*
*  References
*  ==========
*
*  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
*      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
*      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
*      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
*  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
*      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
*      Estimation: Theory, Algorithms and Software, Report
*      UMINF - 94.04, Department of Computing Science, Umea University,
*      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
*      To appear in Numerical Algorithms, 1996.
*
*  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
*      for Solving the Generalized Sylvester Equation and Estimating the
*      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
*      Department of Computing Science, Umea University, S-901 87 Umea,
*      Sweden, December 1993, Revised April 1994, Also as LAPACK working
*      Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
*      1996.
*
*  =====================================================================
*

```
go to the page top

## ctgsja

```USAGE:
alpha, beta, ncycle, info, a, b, u, v, q = NumRu::Lapack.ctgsja( jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q, [:usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, NCYCLE, INFO )

*  Purpose
*  =======
*
*  CTGSJA computes the generalized singular value decomposition (GSVD)
*  of two complex upper triangular (or trapezoidal) matrices A and B.
*
*  On entry, it is assumed that matrices A and B have the following
*  forms, which may be obtained by the preprocessing subroutine CGGSVP
*  from a general M-by-N matrix A and P-by-N matrix B:
*
*               N-K-L  K    L
*     A =    K ( 0    A12  A13 ) if M-K-L >= 0;
*            L ( 0     0   A23 )
*        M-K-L ( 0     0    0  )
*
*             N-K-L  K    L
*     A =  K ( 0    A12  A13 ) if M-K-L < 0;
*        M-K ( 0     0   A23 )
*
*             N-K-L  K    L
*     B =  L ( 0     0   B13 )
*        P-L ( 0     0    0  )
*
*  where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
*  upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
*  otherwise A23 is (M-K)-by-L upper trapezoidal.
*
*  On exit,
*
*         U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R ),
*
*  where U, V and Q are unitary matrices, Z' denotes the conjugate
*  transpose of Z, R is a nonsingular upper triangular matrix, and D1
*  and D2 are ``diagonal'' matrices, which are of the following
*  structures:
*
*  If M-K-L >= 0,
*
*                      K  L
*         D1 =     K ( I  0 )
*                  L ( 0  C )
*              M-K-L ( 0  0 )
*
*                     K  L
*         D2 = L   ( 0  S )
*              P-L ( 0  0 )
*
*                 N-K-L  K    L
*    ( 0 R ) = K (  0   R11  R12 ) K
*              L (  0    0   R22 ) L
*
*  where
*
*    C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
*    S = diag( BETA(K+1),  ... , BETA(K+L) ),
*    C**2 + S**2 = I.
*
*    R is stored in A(1:K+L,N-K-L+1:N) on exit.
*
*  If M-K-L < 0,
*
*                 K M-K K+L-M
*      D1 =   K ( I  0    0   )
*           M-K ( 0  C    0   )
*
*                   K M-K K+L-M
*      D2 =   M-K ( 0  S    0   )
*           K+L-M ( 0  0    I   )
*             P-L ( 0  0    0   )
*
*                 N-K-L  K   M-K  K+L-M
* ( 0 R ) =    K ( 0    R11  R12  R13  )
*            M-K ( 0     0   R22  R23  )
*          K+L-M ( 0     0    0   R33  )
*
*  where
*  C = diag( ALPHA(K+1), ... , ALPHA(M) ),
*  S = diag( BETA(K+1),  ... , BETA(M) ),
*  C**2 + S**2 = I.
*
*  R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
*      (  0  R22 R23 )
*  in B(M-K+1:L,N+M-K-L+1:N) on exit.
*
*  The computation of the unitary transformation matrices U, V or Q
*  is optional.  These matrices may either be formed explicitly, or they
*  may be postmultiplied into input matrices U1, V1, or Q1.
*

*  Arguments
*  =========
*
*  JOBU    (input) CHARACTER*1
*          = 'U':  U must contain a unitary matrix U1 on entry, and
*                  the product U1*U is returned;
*          = 'I':  U is initialized to the unit matrix, and the
*                  unitary matrix U is returned;
*          = 'N':  U is not computed.
*
*  JOBV    (input) CHARACTER*1
*          = 'V':  V must contain a unitary matrix V1 on entry, and
*                  the product V1*V is returned;
*          = 'I':  V is initialized to the unit matrix, and the
*                  unitary matrix V is returned;
*          = 'N':  V is not computed.
*
*  JOBQ    (input) CHARACTER*1
*          = 'Q':  Q must contain a unitary matrix Q1 on entry, and
*                  the product Q1*Q is returned;
*          = 'I':  Q is initialized to the unit matrix, and the
*                  unitary matrix Q is returned;
*          = 'N':  Q is not computed.
*
*  M       (input) INTEGER
*          The number of rows of the matrix A.  M >= 0.
*
*  P       (input) INTEGER
*          The number of rows of the matrix B.  P >= 0.
*
*  N       (input) INTEGER
*          The number of columns of the matrices A and B.  N >= 0.
*
*  K       (input) INTEGER
*  L       (input) INTEGER
*          K and L specify the subblocks in the input matrices A and B:
*          A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N)
*          of A and B, whose GSVD is going to be computed by CTGSJA.
*          See Further Details.
*
*  A       (input/output) COMPLEX array, dimension (LDA,N)
*          On entry, the M-by-N matrix A.
*          On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular
*          matrix R or part of R.  See Purpose for details.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A. LDA >= max(1,M).
*
*  B       (input/output) COMPLEX array, dimension (LDB,N)
*          On entry, the P-by-N matrix B.
*          On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
*          a part of R.  See Purpose for details.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B. LDB >= max(1,P).
*
*  TOLA    (input) REAL
*  TOLB    (input) REAL
*          TOLA and TOLB are the convergence criteria for the Jacobi-
*          Kogbetliantz iteration procedure. Generally, they are the
*          same as used in the preprocessing step, say
*              TOLA = MAX(M,N)*norm(A)*MACHEPS,
*              TOLB = MAX(P,N)*norm(B)*MACHEPS.
*
*  ALPHA   (output) REAL array, dimension (N)
*  BETA    (output) REAL array, dimension (N)
*          On exit, ALPHA and BETA contain the generalized singular
*          value pairs of A and B;
*            ALPHA(1:K) = 1,
*            BETA(1:K)  = 0,
*          and if M-K-L >= 0,
*            ALPHA(K+1:K+L) = diag(C),
*            BETA(K+1:K+L)  = diag(S),
*          or if M-K-L < 0,
*            ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
*            BETA(K+1:M) = S, BETA(M+1:K+L) = 1.
*          Furthermore, if K+L < N,
*            ALPHA(K+L+1:N) = 0
*            BETA(K+L+1:N)  = 0.
*
*  U       (input/output) COMPLEX array, dimension (LDU,M)
*          On entry, if JOBU = 'U', U must contain a matrix U1 (usually
*          the unitary matrix returned by CGGSVP).
*          On exit,
*          if JOBU = 'I', U contains the unitary matrix U;
*          if JOBU = 'U', U contains the product U1*U.
*          If JOBU = 'N', U is not referenced.
*
*  LDU     (input) INTEGER
*          The leading dimension of the array U. LDU >= max(1,M) if
*          JOBU = 'U'; LDU >= 1 otherwise.
*
*  V       (input/output) COMPLEX array, dimension (LDV,P)
*          On entry, if JOBV = 'V', V must contain a matrix V1 (usually
*          the unitary matrix returned by CGGSVP).
*          On exit,
*          if JOBV = 'I', V contains the unitary matrix V;
*          if JOBV = 'V', V contains the product V1*V.
*          If JOBV = 'N', V is not referenced.
*
*  LDV     (input) INTEGER
*          The leading dimension of the array V. LDV >= max(1,P) if
*          JOBV = 'V'; LDV >= 1 otherwise.
*
*  Q       (input/output) COMPLEX array, dimension (LDQ,N)
*          On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
*          the unitary matrix returned by CGGSVP).
*          On exit,
*          if JOBQ = 'I', Q contains the unitary matrix Q;
*          if JOBQ = 'Q', Q contains the product Q1*Q.
*          If JOBQ = 'N', Q is not referenced.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q. LDQ >= max(1,N) if
*          JOBQ = 'Q'; LDQ >= 1 otherwise.
*
*  WORK    (workspace) COMPLEX array, dimension (2*N)
*
*  NCYCLE  (output) INTEGER
*          The number of cycles required for convergence.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*          = 1:  the procedure does not converge after MAXIT cycles.
*
*  Internal Parameters
*  ===================
*
*  MAXIT   INTEGER
*          MAXIT specifies the total loops that the iterative procedure
*          may take. If after MAXIT cycles, the routine fails to
*          converge, we return INFO = 1.
*

*  Further Details
*  ===============
*
*  CTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
*  min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
*  matrix B13 to the form:
*
*           U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1,
*
*  where U1, V1 and Q1 are unitary matrix, and Z' is the conjugate
*  transpose of Z.  C1 and S1 are diagonal matrices satisfying
*
*                C1**2 + S1**2 = I,
*
*  and R1 is an L-by-L nonsingular upper triangular matrix.
*
*  =====================================================================
*

```
go to the page top

## ctgsna

```USAGE:
s, dif, m, work, info = NumRu::Lapack.ctgsna( job, howmny, select, a, b, vl, vr, [:lwork => lwork, :usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO )

*  Purpose
*  =======
*
*  CTGSNA estimates reciprocal condition numbers for specified
*  eigenvalues and/or eigenvectors of a matrix pair (A, B).
*
*  (A, B) must be in generalized Schur canonical form, that is, A and
*  B are both upper triangular.
*

*  Arguments
*  =========
*
*  JOB     (input) CHARACTER*1
*          Specifies whether condition numbers are required for
*          eigenvalues (S) or eigenvectors (DIF):
*          = 'E': for eigenvalues only (S);
*          = 'V': for eigenvectors only (DIF);
*          = 'B': for both eigenvalues and eigenvectors (S and DIF).
*
*  HOWMNY  (input) CHARACTER*1
*          = 'A': compute condition numbers for all eigenpairs;
*          = 'S': compute condition numbers for selected eigenpairs
*                 specified by the array SELECT.
*
*  SELECT  (input) LOGICAL array, dimension (N)
*          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
*          condition numbers are required. To select condition numbers
*          for the corresponding j-th eigenvalue and/or eigenvector,
*          SELECT(j) must be set to .TRUE..
*          If HOWMNY = 'A', SELECT is not referenced.
*
*  N       (input) INTEGER
*          The order of the square matrix pair (A, B). N >= 0.
*
*  A       (input) COMPLEX array, dimension (LDA,N)
*          The upper triangular matrix A in the pair (A,B).
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A. LDA >= max(1,N).
*
*  B       (input) COMPLEX array, dimension (LDB,N)
*          The upper triangular matrix B in the pair (A, B).
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B. LDB >= max(1,N).
*
*  VL      (input) COMPLEX array, dimension (LDVL,M)
*          IF JOB = 'E' or 'B', VL must contain left eigenvectors of
*          (A, B), corresponding to the eigenpairs specified by HOWMNY
*          and SELECT.  The eigenvectors must be stored in consecutive
*          columns of VL, as returned by CTGEVC.
*          If JOB = 'V', VL is not referenced.
*
*  LDVL    (input) INTEGER
*          The leading dimension of the array VL. LDVL >= 1; and
*          If JOB = 'E' or 'B', LDVL >= N.
*
*  VR      (input) COMPLEX array, dimension (LDVR,M)
*          IF JOB = 'E' or 'B', VR must contain right eigenvectors of
*          (A, B), corresponding to the eigenpairs specified by HOWMNY
*          and SELECT.  The eigenvectors must be stored in consecutive
*          columns of VR, as returned by CTGEVC.
*          If JOB = 'V', VR is not referenced.
*
*  LDVR    (input) INTEGER
*          The leading dimension of the array VR. LDVR >= 1;
*          If JOB = 'E' or 'B', LDVR >= N.
*
*  S       (output) REAL array, dimension (MM)
*          If JOB = 'E' or 'B', the reciprocal condition numbers of the
*          selected eigenvalues, stored in consecutive elements of the
*          array.
*          If JOB = 'V', S is not referenced.
*
*  DIF     (output) REAL array, dimension (MM)
*          If JOB = 'V' or 'B', the estimated reciprocal condition
*          numbers of the selected eigenvectors, stored in consecutive
*          elements of the array.
*          If the eigenvalues cannot be reordered to compute DIF(j),
*          DIF(j) is set to 0; this can only occur when the true value
*          would be very small anyway.
*          For each eigenvalue/vector specified by SELECT, DIF stores
*          a Frobenius norm-based estimate of Difl.
*          If JOB = 'E', DIF is not referenced.
*
*  MM      (input) INTEGER
*          The number of elements in the arrays S and DIF. MM >= M.
*
*  M       (output) INTEGER
*          The number of elements of the arrays S and DIF used to store
*          the specified condition numbers; for each selected eigenvalue
*          one element is used. If HOWMNY = 'A', M is set to N.
*
*  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK  (input) INTEGER
*          The dimension of the array WORK. LWORK >= max(1,N).
*          If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
*
*  IWORK   (workspace) INTEGER array, dimension (N+2)
*          If JOB = 'E', IWORK is not referenced.
*
*  INFO    (output) INTEGER
*          = 0: Successful exit
*          < 0: If INFO = -i, the i-th argument had an illegal value
*

*  Further Details
*  ===============
*
*  The reciprocal of the condition number of the i-th generalized
*  eigenvalue w = (a, b) is defined as
*
*          S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v))
*
*  where u and v are the right and left eigenvectors of (A, B)
*  corresponding to w; |z| denotes the absolute value of the complex
*  number, and norm(u) denotes the 2-norm of the vector u. The pair
*  (a, b) corresponds to an eigenvalue w = a/b (= v'Au/v'Bu) of the
*  matrix pair (A, B). If both a and b equal zero, then (A,B) is
*  singular and S(I) = -1 is returned.
*
*  An approximate error bound on the chordal distance between the i-th
*  computed generalized eigenvalue w and the corresponding exact
*  eigenvalue lambda is
*
*          chord(w, lambda) <=   EPS * norm(A, B) / S(I),
*
*  where EPS is the machine precision.
*
*  The reciprocal of the condition number of the right eigenvector u
*  and left eigenvector v corresponding to the generalized eigenvalue w
*  is defined as follows. Suppose
*
*                   (A, B) = ( a   *  ) ( b  *  )  1
*                            ( 0  A22 ),( 0 B22 )  n-1
*                              1  n-1     1 n-1
*
*  Then the reciprocal condition number DIF(I) is
*
*          Difl[(a, b), (A22, B22)]  = sigma-min( Zl )
*
*  where sigma-min(Zl) denotes the smallest singular value of
*
*         Zl = [ kron(a, In-1) -kron(1, A22) ]
*              [ kron(b, In-1) -kron(1, B22) ].
*
*  Here In-1 is the identity matrix of size n-1 and X' is the conjugate
*  transpose of X. kron(X, Y) is the Kronecker product between the
*  matrices X and Y.
*
*  We approximate the smallest singular value of Zl with an upper
*  bound. This is done by CLATDF.
*
*  An approximate error bound for a computed eigenvector VL(i) or
*  VR(i) is given by
*
*                      EPS * norm(A, B) / DIF(i).
*
*  See ref. [2-3] for more details and further references.
*
*  Based on contributions by
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*     Umea University, S-901 87 Umea, Sweden.
*
*  References
*  ==========
*
*  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
*      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
*      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
*      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
*  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
*      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
*      Estimation: Theory, Algorithms and Software, Report
*      UMINF - 94.04, Department of Computing Science, Umea University,
*      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
*      To appear in Numerical Algorithms, 1996.
*
*  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
*      for Solving the Generalized Sylvester Equation and Estimating the
*      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
*      Department of Computing Science, Umea University, S-901 87 Umea,
*      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
*      Note 75.
*      To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
*
*  =====================================================================
*

```
go to the page top

## ctgsy2

```USAGE:
scale, info, c, f, rdsum, rdscal = NumRu::Lapack.ctgsy2( trans, ijob, a, b, c, d, e, f, rdsum, rdscal, [:usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, INFO )

*  Purpose
*  =======
*
*  CTGSY2 solves the generalized Sylvester equation
*
*              A * R - L * B = scale *   C               (1)
*              D * R - L * E = scale * F
*
*  using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
*  (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
*  N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
*  (i.e., (A,D) and (B,E) in generalized Schur form).
*
*  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
*  scaling factor chosen to avoid overflow.
*
*  In matrix notation solving equation (1) corresponds to solve
*  Zx = scale * b, where Z is defined as
*
*         Z = [ kron(In, A)  -kron(B', Im) ]             (2)
*             [ kron(In, D)  -kron(E', Im) ],
*
*  Ik is the identity matrix of size k and X' is the transpose of X.
*  kron(X, Y) is the Kronecker product between the matrices X and Y.
*
*  If TRANS = 'C', y in the conjugate transposed system Z'y = scale*b
*  is solved for, which is equivalent to solve for R and L in
*
*              A' * R  + D' * L   = scale *  C           (3)
*              R  * B' + L  * E'  = scale * -F
*
*  This case is used to compute an estimate of Dif[(A, D), (B, E)] =
*  = sigma_min(Z) using reverse communicaton with CLACON.
*
*  CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL
*  of an upper bound on the separation between to matrix pairs. Then
*  the input (A, D), (B, E) are sub-pencils of two matrix pairs in
*  CTGSYL.
*

*  Arguments
*  =========
*
*  TRANS   (input) CHARACTER*1
*          = 'N', solve the generalized Sylvester equation (1).
*          = 'T': solve the 'transposed' system (3).
*
*  IJOB    (input) INTEGER
*          Specifies what kind of functionality to be performed.
*          =0: solve (1) only.
*          =1: A contribution from this subsystem to a Frobenius
*              norm-based estimate of the separation between two matrix
*              pairs is computed. (look ahead strategy is used).
*          =2: A contribution from this subsystem to a Frobenius
*              norm-based estimate of the separation between two matrix
*              pairs is computed. (SGECON on sub-systems is used.)
*          Not referenced if TRANS = 'T'.
*
*  M       (input) INTEGER
*          On entry, M specifies the order of A and D, and the row
*          dimension of C, F, R and L.
*
*  N       (input) INTEGER
*          On entry, N specifies the order of B and E, and the column
*          dimension of C, F, R and L.
*
*  A       (input) COMPLEX array, dimension (LDA, M)
*          On entry, A contains an upper triangular matrix.
*
*  LDA     (input) INTEGER
*          The leading dimension of the matrix A. LDA >= max(1, M).
*
*  B       (input) COMPLEX array, dimension (LDB, N)
*          On entry, B contains an upper triangular matrix.
*
*  LDB     (input) INTEGER
*          The leading dimension of the matrix B. LDB >= max(1, N).
*
*  C       (input/output) COMPLEX array, dimension (LDC, N)
*          On entry, C contains the right-hand-side of the first matrix
*          equation in (1).
*          On exit, if IJOB = 0, C has been overwritten by the solution
*          R.
*
*  LDC     (input) INTEGER
*          The leading dimension of the matrix C. LDC >= max(1, M).
*
*  D       (input) COMPLEX array, dimension (LDD, M)
*          On entry, D contains an upper triangular matrix.
*
*  LDD     (input) INTEGER
*          The leading dimension of the matrix D. LDD >= max(1, M).
*
*  E       (input) COMPLEX array, dimension (LDE, N)
*          On entry, E contains an upper triangular matrix.
*
*  LDE     (input) INTEGER
*          The leading dimension of the matrix E. LDE >= max(1, N).
*
*  F       (input/output) COMPLEX array, dimension (LDF, N)
*          On entry, F contains the right-hand-side of the second matrix
*          equation in (1).
*          On exit, if IJOB = 0, F has been overwritten by the solution
*          L.
*
*  LDF     (input) INTEGER
*          The leading dimension of the matrix F. LDF >= max(1, M).
*
*  SCALE   (output) REAL
*          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
*          R and L (C and F on entry) will hold the solutions to a
*          slightly perturbed system but the input matrices A, B, D and
*          E have not been changed. If SCALE = 0, R and L will hold the
*          solutions to the homogeneous system with C = F = 0.
*          Normally, SCALE = 1.
*
*  RDSUM   (input/output) REAL
*          On entry, the sum of squares of computed contributions to
*          the Dif-estimate under computation by CTGSYL, where the
*          scaling factor RDSCAL (see below) has been factored out.
*          On exit, the corresponding sum of squares updated with the
*          contributions from the current sub-system.
*          If TRANS = 'T' RDSUM is not touched.
*          NOTE: RDSUM only makes sense when CTGSY2 is called by
*          CTGSYL.
*
*  RDSCAL  (input/output) REAL
*          On entry, scaling factor used to prevent overflow in RDSUM.
*          On exit, RDSCAL is updated w.r.t. the current contributions
*          in RDSUM.
*          If TRANS = 'T', RDSCAL is not touched.
*          NOTE: RDSCAL only makes sense when CTGSY2 is called by
*          CTGSYL.
*
*  INFO    (output) INTEGER
*          On exit, if INFO is set to
*            =0: Successful exit
*            <0: If INFO = -i, input argument number i is illegal.
*            >0: The matrix pairs (A, D) and (B, E) have common or very
*                close eigenvalues.
*

*  Further Details
*  ===============
*
*  Based on contributions by
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*     Umea University, S-901 87 Umea, Sweden.
*
*  =====================================================================
*

```
go to the page top

## ctgsyl

```USAGE:
scale, dif, work, info, c, f = NumRu::Lapack.ctgsyl( trans, ijob, a, b, c, d, e, f, [:lwork => lwork, :usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE CTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO )

*  Purpose
*  =======
*
*  CTGSYL solves the generalized Sylvester equation:
*
*              A * R - L * B = scale * C            (1)
*              D * R - L * E = scale * F
*
*  where R and L are unknown m-by-n matrices, (A, D), (B, E) and
*  (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
*  respectively, with complex entries. A, B, D and E are upper
*  triangular (i.e., (A,D) and (B,E) in generalized Schur form).
*
*  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
*  is an output scaling factor chosen to avoid overflow.
*
*  In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
*  is defined as
*
*         Z = [ kron(In, A)  -kron(B', Im) ]        (2)
*             [ kron(In, D)  -kron(E', Im) ],
*
*  Here Ix is the identity matrix of size x and X' is the conjugate
*  transpose of X. Kron(X, Y) is the Kronecker product between the
*  matrices X and Y.
*
*  If TRANS = 'C', y in the conjugate transposed system Z'*y = scale*b
*  is solved for, which is equivalent to solve for R and L in
*
*              A' * R + D' * L = scale * C           (3)
*              R * B' + L * E' = scale * -F
*
*  This case (TRANS = 'C') is used to compute an one-norm-based estimate
*  of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
*  and (B,E), using CLACON.
*
*  If IJOB >= 1, CTGSYL computes a Frobenius norm-based estimate of
*  Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
*  reciprocal of the smallest singular value of Z.
*
*  This is a level-3 BLAS algorithm.
*

*  Arguments
*  =========
*
*  TRANS   (input) CHARACTER*1
*          = 'N': solve the generalized sylvester equation (1).
*          = 'C': solve the "conjugate transposed" system (3).
*
*  IJOB    (input) INTEGER
*          Specifies what kind of functionality to be performed.
*          =0: solve (1) only.
*          =1: The functionality of 0 and 3.
*          =2: The functionality of 0 and 4.
*          =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
*              (look ahead strategy is used).
*          =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
*              (CGECON on sub-systems is used).
*          Not referenced if TRANS = 'C'.
*
*  M       (input) INTEGER
*          The order of the matrices A and D, and the row dimension of
*          the matrices C, F, R and L.
*
*  N       (input) INTEGER
*          The order of the matrices B and E, and the column dimension
*          of the matrices C, F, R and L.
*
*  A       (input) COMPLEX array, dimension (LDA, M)
*          The upper triangular matrix A.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A. LDA >= max(1, M).
*
*  B       (input) COMPLEX array, dimension (LDB, N)
*          The upper triangular matrix B.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B. LDB >= max(1, N).
*
*  C       (input/output) COMPLEX array, dimension (LDC, N)
*          On entry, C contains the right-hand-side of the first matrix
*          equation in (1) or (3).
*          On exit, if IJOB = 0, 1 or 2, C has been overwritten by
*          the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
*          the solution achieved during the computation of the
*          Dif-estimate.
*
*  LDC     (input) INTEGER
*          The leading dimension of the array C. LDC >= max(1, M).
*
*  D       (input) COMPLEX array, dimension (LDD, M)
*          The upper triangular matrix D.
*
*  LDD     (input) INTEGER
*          The leading dimension of the array D. LDD >= max(1, M).
*
*  E       (input) COMPLEX array, dimension (LDE, N)
*          The upper triangular matrix E.
*
*  LDE     (input) INTEGER
*          The leading dimension of the array E. LDE >= max(1, N).
*
*  F       (input/output) COMPLEX array, dimension (LDF, N)
*          On entry, F contains the right-hand-side of the second matrix
*          equation in (1) or (3).
*          On exit, if IJOB = 0, 1 or 2, F has been overwritten by
*          the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
*          the solution achieved during the computation of the
*          Dif-estimate.
*
*  LDF     (input) INTEGER
*          The leading dimension of the array F. LDF >= max(1, M).
*
*  DIF     (output) REAL
*          On exit DIF is the reciprocal of a lower bound of the
*          reciprocal of the Dif-function, i.e. DIF is an upper bound of
*          Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2).
*          IF IJOB = 0 or TRANS = 'C', DIF is not referenced.
*
*  SCALE   (output) REAL
*          On exit SCALE is the scaling factor in (1) or (3).
*          If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
*          to a slightly perturbed system but the input matrices A, B,
*          D and E have not been changed. If SCALE = 0, R and L will
*          hold the solutions to the homogenious system with C = F = 0.
*
*  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK. LWORK > = 1.
*          If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by XERBLA.
*
*  IWORK   (workspace) INTEGER array, dimension (M+N+2)
*
*  INFO    (output) INTEGER
*            =0: successful exit
*            <0: If INFO = -i, the i-th argument had an illegal value.
*            >0: (A, D) and (B, E) have common or very close
*                eigenvalues.
*

*  Further Details
*  ===============
*
*  Based on contributions by
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*     Umea University, S-901 87 Umea, Sweden.
*
*  [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
*      for Solving the Generalized Sylvester Equation and Estimating the
*      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
*      Department of Computing Science, Umea University, S-901 87 Umea,
*      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
*      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
*      No 1, 1996.
*
*  [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
*      Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
*      Appl., 15(4):1045-1060, 1994.
*
*  [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
*      Condition Estimators for Solving the Generalized Sylvester
*      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
*      July 1989, pp 745-751.
*
*  =====================================================================
*  Replaced various illegal calls to CCOPY by calls to CLASET.
*  Sven Hammarling, 1/5/02.
*

```
go to the page top
back to matrix types
back to data types