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

dtgevc

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

dtgex2

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

dtgexc

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

dtgsen

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

dtgsja

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

dtgsna

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

dtgsy2

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

dtgsyl

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