From 879f74caf4ab9128e50a0f61ab0d4f06156beadf Mon Sep 17 00:00:00 2001
From: Olivier Stasse <ostasse@laas.fr>
Date: Fri, 23 Sep 2016 19:24:01 +0200
Subject: [PATCH] Fix and first running test with the Fortran solution.
---
src/Fortran/CMakeLists.txt | 5 +-
src/Fortran/{dgges.f => gdgges.f} | 18 +-
src/Fortran/gdtgex2.f | 697 ++++++++++++++
src/Fortran/gdtgexc.f | 544 +++++++++++
src/Fortran/gdtgsen.f | 866 ++++++++++++++++++
.../OptimalControllerSolver.cpp | 4 +-
tests/CMakeLists.txt | 2 +
7 files changed, 2124 insertions(+), 12 deletions(-)
rename src/Fortran/{dgges.f => gdgges.f} (97%)
create mode 100644 src/Fortran/gdtgex2.f
create mode 100644 src/Fortran/gdtgexc.f
create mode 100644 src/Fortran/gdtgsen.f
diff --git a/src/Fortran/CMakeLists.txt b/src/Fortran/CMakeLists.txt
index db17ae78..7daa78ad 100644
--- a/src/Fortran/CMakeLists.txt
+++ b/src/Fortran/CMakeLists.txt
@@ -4,10 +4,13 @@ SET(SOURCES_FORTRAN
dgesvd.f
dgesvx.f
dgetrf.f
- dgges.f
+ gdgges.f
dlamch.f
dlapy2.f
dgetrf2.f
+ gdtgsen.f
+ gdtgexc.f
+ gdtgex2.f
gdgemm.f
gsgemm.f
)
diff --git a/src/Fortran/dgges.f b/src/Fortran/gdgges.f
similarity index 97%
rename from src/Fortran/dgges.f
rename to src/Fortran/gdgges.f
index 76d6d399..1a9a2f6d 100644
--- a/src/Fortran/dgges.f
+++ b/src/Fortran/gdgges.f
@@ -1,4 +1,4 @@
-*> \brief <b> DGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
+*> \brief <b> GDGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
*
* =========== DOCUMENTATION ===========
*
@@ -6,7 +6,7 @@
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download DGGES + dependencies
+*> Download GDGGES + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgges.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgges.f">
@@ -18,7 +18,7 @@
* Definition:
* ===========
*
-* SUBROUTINE DGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
+* SUBROUTINE GDGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
* SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
* LDVSR, WORK, LWORK, BWORK, INFO )
*
@@ -43,7 +43,7 @@
*>
*> \verbatim
*>
-*> DGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
+*> GDGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
*> the generalized eigenvalues, the generalized real Schur form (S,T),
*> optionally, the left and/or right matrices of Schur vectors (VSL and
*> VSR). This gives the generalized Schur factorization
@@ -264,8 +264,8 @@
*> eigenvalues in the Generalized Schur form no
*> longer satisfy SELCTG=.TRUE. This could also
*> be caused due to scaling.
-*> =N+3: reordering failed in DTGSEN.
-*> \endverbatim
+*> =N+3: reordering failed in GDTGSEN.
+*> \endverbatnim
*
* Authors:
* ========
@@ -280,7 +280,7 @@
*> \ingroup doubleGEeigen
*
* =====================================================================
- SUBROUTINE DGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
+ SUBROUTINE GDGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,LDB,
$ SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
$ LDVSR, WORK, LWORK, BWORK, INFO )
*
@@ -325,7 +325,7 @@
* ..
* .. External Subroutines ..
EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
- $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
+ $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR,GDTGSEN,
$ XERBLA
* ..
* .. External Functions ..
@@ -557,7 +557,7 @@
BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
10 CONTINUE
*
- CALL DTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
+ CALL GDTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,ALPHAR,
$ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
$ PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
$ IERR )
diff --git a/src/Fortran/gdtgex2.f b/src/Fortran/gdtgex2.f
new file mode 100644
index 00000000..51cf797e
--- /dev/null
+++ b/src/Fortran/gdtgex2.f
@@ -0,0 +1,697 @@
+*> \brief \b GDTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download GDTGEX2 + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgex2.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgex2.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgex2.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+* LDZ, J1, N1, N2, WORK, LWORK, INFO )
+*
+* .. Scalar Arguments ..
+* LOGICAL WANTQ, WANTZ
+* INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
+* $ WORK( * ), Z( LDZ, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> GDTGEX2 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)**T = Q(out) * A(out) * Z(out)**T
+*> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
+*>
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] WANTQ
+*> \verbatim
+*> WANTQ is LOGICAL
+*> .TRUE. : update the left transformation matrix Q;
+*> .FALSE.: do not update Q.
+*> \endverbatim
+*>
+*> \param[in] WANTZ
+*> \verbatim
+*> WANTZ is LOGICAL
+*> .TRUE. : update the right transformation matrix Z;
+*> .FALSE.: do not update Z.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A and B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is DOUBLE PRECISION array, dimensions (LDA,N)
+*> On entry, the matrix A in the pair (A, B).
+*> On exit, the updated matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is DOUBLE PRECISION array, dimensions (LDB,N)
+*> On entry, the matrix B in the pair (A, B).
+*> On exit, the updated matrix B.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] Q
+*> \verbatim
+*> Q is 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..
+*> \endverbatim
+*>
+*> \param[in] LDQ
+*> \verbatim
+*> LDQ is INTEGER
+*> The leading dimension of the array Q. LDQ >= 1.
+*> If WANTQ = .TRUE., LDQ >= N.
+*> \endverbatim
+*>
+*> \param[in,out] Z
+*> \verbatim
+*> Z is 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..
+*> \endverbatim
+*>
+*> \param[in] LDZ
+*> \verbatim
+*> LDZ is INTEGER
+*> The leading dimension of the array Z. LDZ >= 1.
+*> If WANTZ = .TRUE., LDZ >= N.
+*> \endverbatim
+*>
+*> \param[in] J1
+*> \verbatim
+*> J1 is INTEGER
+*> The index to the first block (A11, B11). 1 <= J1 <= N.
+*> \endverbatim
+*>
+*> \param[in] N1
+*> \verbatim
+*> N1 is INTEGER
+*> The order of the first block (A11, B11). N1 = 0, 1 or 2.
+*> \endverbatim
+*>
+*> \param[in] N2
+*> \verbatim
+*> N2 is INTEGER
+*> The order of the second block (A22, B22). N2 = 0, 1 or 2.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)).
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> LWORK >= MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 )
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is 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).
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2015
+*
+*> \ingroup doubleGEauxiliary
+*
+*> \par Further Details:
+* =====================
+*>
+*> 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.
+*
+*> \par Contributors:
+* ==================
+*>
+*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
+*> Umea University, S-901 87 Umea, Sweden.
+*
+*> \par References:
+* ================
+*>
+*> \verbatim
+*>
+*> [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.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, J1, N1, N2, WORK, LWORK, INFO )
+*
+* -- LAPACK auxiliary routine (version 3.6.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2015
+*
+* .. Scalar Arguments ..
+ LOGICAL WANTQ, WANTZ
+ INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
+ $ WORK( * ), Z( LDZ, * )
+* ..
+*
+* =====================================================================
+* Replaced various illegal calls to DCOPY by calls to DLASET, or by DO
+* loops. Sven Hammarling, 1/5/02.
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+ DOUBLE PRECISION TWENTY
+ PARAMETER ( TWENTY = 2.0D+01 )
+ INTEGER LDST
+ PARAMETER ( LDST = 4 )
+ LOGICAL WANDS
+ PARAMETER ( WANDS = .TRUE. )
+* ..
+* .. Local Scalars ..
+ LOGICAL DTRONG, WEAK
+ INTEGER I, IDUM, LINFO, M
+ DOUBLE PRECISION BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,
+ $ F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS
+* ..
+* .. Local Arrays ..
+ INTEGER IWORK( LDST )
+ DOUBLE PRECISION AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
+ $ IRCOP( LDST, LDST ), LI( LDST, LDST ),
+ $ LICOP( LDST, LDST ), S( LDST, LDST ),
+ $ SCPY( LDST, LDST ), T( LDST, LDST ),
+ $ TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
+* ..
+* .. External Functions ..
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL DLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL GDGEMM, DGEQR2, DGERQ2, DLACPY, DLAGV2, DLARTG,
+ $ DLASET, DLASSQ, DORG2R, DORGR2, DORM2R, DORMR2,
+ $ DROT, DSCAL, DTGSY2
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, SQRT
+* ..
+* .. Executable Statements ..
+*
+ INFO = 0
+*
+* Quick return if possible
+*
+ IF( N.LE.1 .OR. N1.LE.0 .OR. N2.LE.0 )
+ $ RETURN
+ IF( N1.GT.N .OR. ( J1+N1 ).GT.N )
+ $ RETURN
+ M = N1 + N2
+ IF( LWORK.LT.MAX( 1, N*M, M*M*2 ) ) THEN
+ INFO = -16
+ WORK( 1 ) = MAX( 1, N*M, M*M*2 )
+ RETURN
+ END IF
+*
+ WEAK = .FALSE.
+ DTRONG = .FALSE.
+*
+* Make a local copy of selected block
+*
+ CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, LI, LDST )
+ CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, IR, LDST )
+ CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
+ CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
+*
+* Compute threshold for testing acceptance of swapping.
+*
+ EPS = DLAMCH( 'P' )
+ SMLNUM = DLAMCH( 'S' ) / EPS
+ DSCALE = ZERO
+ DSUM = ONE
+ CALL DLACPY( 'Full', M, M, S, LDST, WORK, M )
+ CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )
+ CALL DLACPY( 'Full', M, M, T, LDST, WORK, M )
+ CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )
+ DNORM = DSCALE*SQRT( DSUM )
+*
+* THRES has been changed from
+* THRESH = MAX( TEN*EPS*SA, SMLNUM )
+* to
+* THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
+* on 04/01/10.
+* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
+* Jim Demmel and Guillaume Revy. See forum post 1783.
+*
+ THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM )
+*
+ IF( M.EQ.2 ) THEN
+*
+* CASE 1: Swap 1-by-1 and 1-by-1 blocks.
+*
+* Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
+* using Givens rotations and perform the swap tentatively.
+*
+ F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
+ G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
+ SB = ABS( T( 2, 2 ) )
+ SA = ABS( S( 2, 2 ) )
+ CALL DLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM )
+ IR( 2, 1 ) = -IR( 1, 2 )
+ IR( 2, 2 ) = IR( 1, 1 )
+ CALL DROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, IR( 1, 1 ),
+ $ IR( 2, 1 ) )
+ CALL DROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, IR( 1, 1 ),
+ $ IR( 2, 1 ) )
+ IF( SA.GE.SB ) THEN
+ CALL DLARTG( S( 1, 1 ), S( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
+ $ DDUM )
+ ELSE
+ CALL DLARTG( T( 1, 1 ), T( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
+ $ DDUM )
+ END IF
+ CALL DROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, LI( 1, 1 ),
+ $ LI( 2, 1 ) )
+ CALL DROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, LI( 1, 1 ),
+ $ LI( 2, 1 ) )
+ LI( 2, 2 ) = LI( 1, 1 )
+ LI( 1, 2 ) = -LI( 2, 1 )
+*
+* Weak stability test:
+* |S21| + |T21| <= O(EPS * F-norm((S, T)))
+*
+ WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
+ WEAK = WS.LE.THRESH
+ IF( .NOT.WEAK )
+ $ GO TO 70
+*
+ IF( WANDS ) THEN
+*
+* Strong stability test:
+* F-norm((A-QL**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A,B)))
+*
+ CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
+ $ M )
+ CALL GDGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, S, LDST,ZERO,
+ $ WORK, M )
+ CALL GDGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST,ONE,
+ $ WORK( M*M+1 ), M )
+ DSCALE = ZERO
+ DSUM = ONE
+ CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
+*
+ CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
+ $ M )
+ CALL GDGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, T, LDST,ZERO,
+ $ WORK, M )
+ CALL GDGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST,ONE,
+ $ WORK( M*M+1 ), M )
+ CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
+ SS = DSCALE*SQRT( DSUM )
+ DTRONG = SS.LE.THRESH
+ IF( .NOT.DTRONG )
+ $ GO TO 70
+ END IF
+*
+* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
+* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
+*
+ CALL DROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, IR( 1, 1 ),
+ $ IR( 2, 1 ) )
+ CALL DROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, IR( 1, 1 ),
+ $ IR( 2, 1 ) )
+ CALL DROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA,
+ $ LI( 1, 1 ), LI( 2, 1 ) )
+ CALL DROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB,
+ $ LI( 1, 1 ), LI( 2, 1 ) )
+*
+* Set N1-by-N2 (2,1) - blocks to ZERO.
+*
+ A( J1+1, J1 ) = ZERO
+ B( J1+1, J1 ) = ZERO
+*
+* Accumulate transformations into Q and Z if requested.
+*
+ IF( WANTZ )
+ $ CALL DROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, IR( 1, 1 ),
+ $ IR( 2, 1 ) )
+ IF( WANTQ )
+ $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, LI( 1, 1 ),
+ $ LI( 2, 1 ) )
+*
+* Exit with INFO = 0 if swap was successfully performed.
+*
+ RETURN
+*
+ ELSE
+*
+* CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
+* and 2-by-2 blocks.
+*
+* Solve the generalized Sylvester equation
+* S11 * R - L * S22 = SCALE * S12
+* T11 * R - L * T22 = SCALE * T12
+* for R and L. Solutions in LI and IR.
+*
+ CALL DLACPY( 'Full', N1, N2, T( 1, N1+1 ), LDST, LI, LDST )
+ CALL DLACPY( 'Full', N1, N2, S( 1, N1+1 ), LDST,
+ $ IR( N2+1, N1+1 ), LDST )
+ CALL DTGSY2( 'N', 0, N1, N2, S, LDST, S( N1+1, N1+1 ), LDST,
+ $ IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ),
+ $ LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM,
+ $ LINFO )
+*
+* Compute orthogonal matrix QL:
+*
+* QL**T * LI = [ TL ]
+* [ 0 ]
+* where
+* LI = [ -L ]
+* [ SCALE * identity(N2) ]
+*
+ DO 10 I = 1, N2
+ CALL DSCAL( N1, -ONE, LI( 1, I ), 1 )
+ LI( N1+I, I ) = SCALE
+ 10 CONTINUE
+ CALL DGEQR2( M, N2, LI, LDST, TAUL, WORK, LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+ CALL DORG2R( M, M, N2, LI, LDST, TAUL, WORK, LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+*
+* Compute orthogonal matrix RQ:
+*
+* IR * RQ**T = [ 0 TR],
+*
+* where IR = [ SCALE * identity(N1), R ]
+*
+ DO 20 I = 1, N1
+ IR( N2+I, I ) = SCALE
+ 20 CONTINUE
+ CALL DGERQ2( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+ CALL DORGR2( M, M, N1, IR, LDST, TAUR, WORK, LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+*
+* Perform the swapping tentatively:
+*
+ CALL GDGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, S, LDST, ZERO,
+ $ WORK, M )
+ CALL GDGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO,S,
+ $ LDST )
+ CALL GDGEMM( 'T', 'N', M, M, M, ONE, LI, LDST, T, LDST, ZERO,
+ $ WORK, M )
+ CALL GDGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO,T,
+ $ LDST )
+ CALL DLACPY( 'F', M, M, S, LDST, SCPY, LDST )
+ CALL DLACPY( 'F', M, M, T, LDST, TCPY, LDST )
+ CALL DLACPY( 'F', M, M, IR, LDST, IRCOP, LDST )
+ CALL DLACPY( 'F', M, M, LI, LDST, LICOP, LDST )
+*
+* Triangularize the B-part by an RQ factorization.
+* Apply transformation (from left) to A-part, giving S.
+*
+ CALL DGERQ2( M, M, T, LDST, TAUR, WORK, LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+ CALL DORMR2( 'R', 'T', M, M, M, T, LDST, TAUR, S, LDST, WORK,
+ $ LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+ CALL DORMR2( 'L', 'N', M, M, M, T, LDST, TAUR, IR, LDST, WORK,
+ $ LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+*
+* Compute F-norm(S21) in BRQA21. (T21 is 0.)
+*
+ DSCALE = ZERO
+ DSUM = ONE
+ DO 30 I = 1, N2
+ CALL DLASSQ( N1, S( N2+1, I ), 1, DSCALE, DSUM )
+ 30 CONTINUE
+ BRQA21 = DSCALE*SQRT( DSUM )
+*
+* Triangularize the B-part by a QR factorization.
+* Apply transformation (from right) to A-part, giving S.
+*
+ CALL DGEQR2( M, M, TCPY, LDST, TAUL, WORK, LINFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+ CALL DORM2R( 'L', 'T', M, M, M, TCPY, LDST, TAUL, SCPY, LDST,
+ $ WORK, INFO )
+ CALL DORM2R( 'R', 'N', M, M, M, TCPY, LDST, TAUL, LICOP, LDST,
+ $ WORK, INFO )
+ IF( LINFO.NE.0 )
+ $ GO TO 70
+*
+* Compute F-norm(S21) in BQRA21. (T21 is 0.)
+*
+ DSCALE = ZERO
+ DSUM = ONE
+ DO 40 I = 1, N2
+ CALL DLASSQ( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM )
+ 40 CONTINUE
+ BQRA21 = DSCALE*SQRT( DSUM )
+*
+* Decide which method to use.
+* Weak stability test:
+* F-norm(S21) <= O(EPS * F-norm((S, T)))
+*
+ IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN
+ CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST )
+ CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST )
+ CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )
+ CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST )
+ ELSE IF( BRQA21.GE.THRESH ) THEN
+ GO TO 70
+ END IF
+*
+* Set lower triangle of B-part to zero
+*
+ CALL DLASET( 'Lower', M-1, M-1, ZERO, ZERO, T(2,1), LDST )
+*
+ IF( WANDS ) THEN
+*
+* Strong stability test:
+* F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,B)))
+*
+ CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
+ $ M )
+ CALL GDGEMM( 'N', 'N', M, M, M,ONE, LI, LDST, S, LDST, ZERO,
+ $ WORK, M )
+ CALL GDGEMM( 'N', 'N', M, M, M,-ONE, WORK, M, IR, LDST, ONE,
+ $ WORK( M*M+1 ), M )
+ DSCALE = ZERO
+ DSUM = ONE
+ CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
+*
+ CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
+ $ M )
+ CALL GDGEMM( 'N', 'N', M, M, M,ONE, LI, LDST, T, LDST, ZERO,
+ $ WORK, M )
+ CALL GDGEMM( 'N', 'N', M, M, M,-ONE, WORK, M, IR, LDST, ONE,
+ $ WORK( M*M+1 ), M )
+ CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
+ SS = DSCALE*SQRT( DSUM )
+ DTRONG = ( SS.LE.THRESH )
+ IF( .NOT.DTRONG )
+ $ GO TO 70
+*
+ END IF
+*
+* If the swap is accepted ("weakly" and "strongly"), apply the
+* transformations and set N1-by-N2 (2,1)-block to zero.
+*
+ CALL DLASET( 'Full', N1, N2, ZERO, ZERO, S(N2+1,1), LDST )
+*
+* copy back M-by-M diagonal block starting at index J1 of (A, B)
+*
+ CALL DLACPY( 'F', M, M, S, LDST, A( J1, J1 ), LDA )
+ CALL DLACPY( 'F', M, M, T, LDST, B( J1, J1 ), LDB )
+ CALL DLASET( 'Full', LDST, LDST, ZERO, ZERO, T, LDST )
+*
+* Standardize existing 2-by-2 blocks.
+*
+ CALL DLASET( 'Full', M, M, ZERO, ZERO, WORK, M )
+ WORK( 1 ) = ONE
+ T( 1, 1 ) = ONE
+ IDUM = LWORK - M*M - 2
+ IF( N2.GT.1 ) THEN
+ CALL DLAGV2( A( J1, J1 ), LDA, B( J1, J1 ), LDB, AR, AI, BE,
+ $ WORK( 1 ), WORK( 2 ), T( 1, 1 ), T( 2, 1 ) )
+ WORK( M+1 ) = -WORK( 2 )
+ WORK( M+2 ) = WORK( 1 )
+ T( N2, N2 ) = T( 1, 1 )
+ T( 1, 2 ) = -T( 2, 1 )
+ END IF
+ WORK( M*M ) = ONE
+ T( M, M ) = ONE
+*
+ IF( N1.GT.1 ) THEN
+ CALL DLAGV2( A( J1+N2, J1+N2 ), LDA, B( J1+N2, J1+N2 ), LDB,
+ $ TAUR, TAUL, WORK( M*M+1 ), WORK( N2*M+N2+1 ),
+ $ WORK( N2*M+N2+2 ), T( N2+1, N2+1 ),
+ $ T( M, M-1 ) )
+ WORK( M*M ) = WORK( N2*M+N2+1 )
+ WORK( M*M-1 ) = -WORK( N2*M+N2+2 )
+ T( M, M ) = T( N2+1, N2+1 )
+ T( M-1, M ) = -T( M, M-1 )
+ END IF
+ CALL GDGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, A( J1, J1+N2),
+ $ LDA, ZERO, WORK( M*M+1 ), N2 )
+ CALL DLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, A( J1, J1+N2 ),
+ $ LDA )
+ CALL GDGEMM( 'T', 'N', N2, N1, N2, ONE, WORK, M, B( J1, J1+N2),
+ $ LDB, ZERO, WORK( M*M+1 ), N2 )
+ CALL DLACPY( 'Full', N2, N1, WORK( M*M+1 ), N2, B( J1, J1+N2 ),
+ $ LDB )
+ CALL GDGEMM( 'N', 'N', M, M, M, ONE, LI, LDST, WORK, M, ZERO,
+ $ WORK( M*M+1 ), M )
+ CALL DLACPY( 'Full', M, M, WORK( M*M+1 ), M, LI, LDST )
+ CALL GDGEMM( 'N', 'N', N2, N1, N1, ONE, A( J1, J1+N2 ), LDA,
+ $ T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
+ CALL DLACPY( 'Full', N2, N1, WORK, N2, A( J1, J1+N2 ), LDA )
+ CALL GDGEMM( 'N', 'N', N2, N1, N1, ONE, B( J1, J1+N2 ), LDB,
+ $ T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
+ CALL DLACPY( 'Full', N2, N1, WORK, N2, B( J1, J1+N2 ), LDB )
+ CALL GDGEMM( 'T', 'N', M, M, M, ONE, IR, LDST, T, LDST, ZERO,
+ $ WORK, M )
+ CALL DLACPY( 'Full', M, M, WORK, M, IR, LDST )
+*
+* Accumulate transformations into Q and Z if requested.
+*
+ IF( WANTQ ) THEN
+ CALL GDGEMM( 'N', 'N', N, M, M, ONE, Q( 1, J1 ), LDQ, LI,
+ $ LDST, ZERO, WORK, N )
+ CALL DLACPY( 'Full', N, M, WORK, N, Q( 1, J1 ), LDQ )
+*
+ END IF
+*
+ IF( WANTZ ) THEN
+ CALL GDGEMM( 'N', 'N', N, M, M, ONE, Z( 1, J1 ), LDZ, IR,
+ $ LDST, ZERO, WORK, N )
+ CALL DLACPY( 'Full', N, M, WORK, N, Z( 1, J1 ), LDZ )
+*
+ END IF
+*
+* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
+* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
+*
+ I = J1 + M
+ IF( I.LE.N ) THEN
+ CALL GDGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
+ $ A( J1, I ), LDA, ZERO, WORK, M )
+ CALL DLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA )
+ CALL GDGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST,
+ $ B( J1, I ), LDB, ZERO, WORK, M )
+ CALL DLACPY( 'Full', M, N-I+1, WORK, M, B( J1, I ), LDB )
+ END IF
+ I = J1 - 1
+ IF( I.GT.0 ) THEN
+ CALL GDGEMM( 'N', 'N', I, M, M, ONE, A( 1, J1 ), LDA, IR,
+ $ LDST, ZERO, WORK, I )
+ CALL DLACPY( 'Full', I, M, WORK, I, A( 1, J1 ), LDA )
+ CALL GDGEMM( 'N', 'N', I, M, M, ONE, B( 1, J1 ), LDB, IR,
+ $ LDST, ZERO, WORK, I )
+ CALL DLACPY( 'Full', I, M, WORK, I, B( 1, J1 ), LDB )
+ END IF
+*
+* Exit with INFO = 0 if swap was successfully performed.
+*
+ RETURN
+*
+ END IF
+*
+* Exit with INFO = 1 if swap was rejected.
+*
+ 70 CONTINUE
+*
+ INFO = 1
+ RETURN
+*
+* End of GDTGEX2
+*
+ END
diff --git a/src/Fortran/gdtgexc.f b/src/Fortran/gdtgexc.f
new file mode 100644
index 00000000..793ee90e
--- /dev/null
+++ b/src/Fortran/gdtgexc.f
@@ -0,0 +1,544 @@
+*> \brief \b GDTGEXC
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download GDTGEXC + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgexc.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgexc.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgexc.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE GDTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+* LDZ, IFST, ILST, WORK, LWORK, INFO )
+*
+* .. Scalar Arguments ..
+* LOGICAL WANTQ, WANTZ
+* INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
+* $ WORK( * ), Z( LDZ, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> GDTGEXC reorders the generalized real Schur decomposition of a real
+*> matrix pair (A,B) using an orthogonal equivalence transformation
+*>
+*> (A, B) = Q * (A, B) * Z**T,
+*>
+*> 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)**T = Q(out) * A(out) * Z(out)**T
+*> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
+*>
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] WANTQ
+*> \verbatim
+*> WANTQ is LOGICAL
+*> .TRUE. : update the left transformation matrix Q;
+*> .FALSE.: do not update Q.
+*> \endverbatim
+*>
+*> \param[in] WANTZ
+*> \verbatim
+*> WANTZ is LOGICAL
+*> .TRUE. : update the right transformation matrix Z;
+*> .FALSE.: do not update Z.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A and B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is 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.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is 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).
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] Q
+*> \verbatim
+*> Q is 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.
+*> \endverbatim
+*>
+*> \param[in] LDQ
+*> \verbatim
+*> LDQ is INTEGER
+*> The leading dimension of the array Q. LDQ >= 1.
+*> If WANTQ = .TRUE., LDQ >= N.
+*> \endverbatim
+*>
+*> \param[in,out] Z
+*> \verbatim
+*> Z is 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.
+*> \endverbatim
+*>
+*> \param[in] LDZ
+*> \verbatim
+*> LDZ is INTEGER
+*> The leading dimension of the array Z. LDZ >= 1.
+*> If WANTZ = .TRUE., LDZ >= N.
+*> \endverbatim
+*>
+*> \param[in,out] IFST
+*> \verbatim
+*> IFST is INTEGER
+*> \endverbatim
+*>
+*> \param[in,out] ILST
+*> \verbatim
+*> ILST is 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.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is 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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is 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.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup doubleGEcomputational
+*
+*> \par Contributors:
+* ==================
+*>
+*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
+*> Umea University, S-901 87 Umea, Sweden.
+*
+*> \par References:
+* ================
+*>
+*> \verbatim
+*>
+*> [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.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE GDTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, IFST, ILST, WORK, LWORK, INFO )
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ LOGICAL WANTQ, WANTZ
+ INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
+ $ WORK( * ), Z( LDZ, * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO
+ PARAMETER ( ZERO = 0.0D+0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL LQUERY
+ INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
+* ..
+* .. External Subroutines ..
+ EXTERNAL GDTGEX2, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX
+* ..
+* .. Executable Statements ..
+*
+* Decode and test input arguments.
+*
+ INFO = 0
+ LQUERY = ( LWORK.EQ.-1 )
+ IF( N.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -5
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
+ INFO = -9
+ ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
+ INFO = -11
+ ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
+ INFO = -12
+ ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
+ INFO = -13
+ END IF
+*
+ IF( INFO.EQ.0 ) THEN
+ IF( N.LE.1 ) THEN
+ LWMIN = 1
+ ELSE
+ LWMIN = 4*N + 16
+ END IF
+ WORK(1) = LWMIN
+*
+ IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
+ INFO = -15
+ END IF
+ END IF
+*
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'GDTGEXC', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.LE.1 )
+ $ RETURN
+*
+* Determine the first row of the specified block and find out
+* if it is 1-by-1 or 2-by-2.
+*
+ IF( IFST.GT.1 ) THEN
+ IF( A( IFST, IFST-1 ).NE.ZERO )
+ $ IFST = IFST - 1
+ END IF
+ NBF = 1
+ IF( IFST.LT.N ) THEN
+ IF( A( IFST+1, IFST ).NE.ZERO )
+ $ NBF = 2
+ END IF
+*
+* Determine the first row of the final block
+* and find out if it is 1-by-1 or 2-by-2.
+*
+ IF( ILST.GT.1 ) THEN
+ IF( A( ILST, ILST-1 ).NE.ZERO )
+ $ ILST = ILST - 1
+ END IF
+ NBL = 1
+ IF( ILST.LT.N ) THEN
+ IF( A( ILST+1, ILST ).NE.ZERO )
+ $ NBL = 2
+ END IF
+ IF( IFST.EQ.ILST )
+ $ RETURN
+*
+ IF( IFST.LT.ILST ) THEN
+*
+* Update ILST.
+*
+ IF( NBF.EQ.2 .AND. NBL.EQ.1 )
+ $ ILST = ILST - 1
+ IF( NBF.EQ.1 .AND. NBL.EQ.2 )
+ $ ILST = ILST + 1
+*
+ HERE = IFST
+*
+ 10 CONTINUE
+*
+* Swap with next one below.
+*
+ IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
+*
+* Current block either 1-by-1 or 2-by-2.
+*
+ NBNEXT = 1
+ IF( HERE+NBF+1.LE.N ) THEN
+ IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
+ $ NBNEXT = 2
+ END IF
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE + NBNEXT
+*
+* Test if 2-by-2 block breaks into two 1-by-1 blocks.
+*
+ IF( NBF.EQ.2 ) THEN
+ IF( A( HERE+1, HERE ).EQ.ZERO )
+ $ NBF = 3
+ END IF
+*
+ ELSE
+*
+* Current block consists of two 1-by-1 blocks, each of which
+* must be swapped individually.
+*
+ NBNEXT = 1
+ IF( HERE+3.LE.N ) THEN
+ IF( A( HERE+3, HERE+2 ).NE.ZERO )
+ $ NBNEXT = 2
+ END IF
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ IF( NBNEXT.EQ.1 ) THEN
+*
+* Swap two 1-by-1 blocks.
+*
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, HERE, 1, 1, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE + 1
+*
+ ELSE
+*
+* Recompute NBNEXT in case of 2-by-2 split.
+*
+ IF( A( HERE+2, HERE+1 ).EQ.ZERO )
+ $ NBNEXT = 1
+ IF( NBNEXT.EQ.2 ) THEN
+*
+* 2-by-2 block did not split.
+*
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
+ $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
+ $ INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE + 2
+ ELSE
+*
+* 2-by-2 block did split.
+*
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
+ $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE + 1
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
+ $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE + 1
+ END IF
+*
+ END IF
+ END IF
+ IF( HERE.LT.ILST )
+ $ GO TO 10
+ ELSE
+ HERE = IFST
+*
+ 20 CONTINUE
+*
+* Swap with next one below.
+*
+ IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
+*
+* Current block either 1-by-1 or 2-by-2.
+*
+ NBNEXT = 1
+ IF( HERE.GE.3 ) THEN
+ IF( A( HERE-1, HERE-2 ).NE.ZERO )
+ $ NBNEXT = 2
+ END IF
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
+ $ INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE - NBNEXT
+*
+* Test if 2-by-2 block breaks into two 1-by-1 blocks.
+*
+ IF( NBF.EQ.2 ) THEN
+ IF( A( HERE+1, HERE ).EQ.ZERO )
+ $ NBF = 3
+ END IF
+*
+ ELSE
+*
+* Current block consists of two 1-by-1 blocks, each of which
+* must be swapped individually.
+*
+ NBNEXT = 1
+ IF( HERE.GE.3 ) THEN
+ IF( A( HERE-1, HERE-2 ).NE.ZERO )
+ $ NBNEXT = 2
+ END IF
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
+ $ INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ IF( NBNEXT.EQ.1 ) THEN
+*
+* Swap two 1-by-1 blocks.
+*
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+ $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE - 1
+ ELSE
+*
+* Recompute NBNEXT in case of 2-by-2 split.
+*
+ IF( A( HERE, HERE-1 ).EQ.ZERO )
+ $ NBNEXT = 1
+ IF( NBNEXT.EQ.2 ) THEN
+*
+* 2-by-2 block did not split.
+*
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
+ $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE - 2
+ ELSE
+*
+* 2-by-2 block did split.
+*
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
+ $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE - 1
+ CALL GDTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
+ $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
+ IF( INFO.NE.0 ) THEN
+ ILST = HERE
+ RETURN
+ END IF
+ HERE = HERE - 1
+ END IF
+ END IF
+ END IF
+ IF( HERE.GT.ILST )
+ $ GO TO 20
+ END IF
+ ILST = HERE
+ WORK( 1 ) = LWMIN
+ RETURN
+*
+* End of GDTGEXC
+*
+ END
diff --git a/src/Fortran/gdtgsen.f b/src/Fortran/gdtgsen.f
new file mode 100644
index 00000000..2109a058
--- /dev/null
+++ b/src/Fortran/gdtgsen.f
@@ -0,0 +1,866 @@
+*> \brief \b GDTGSEN
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download GDTGSEN + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsen.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsen.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsen.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE GDTGSEN( 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 )
+*
+* .. Scalar Arguments ..
+* LOGICAL WANTQ, WANTZ
+* INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
+* $ M, N
+* DOUBLE PRECISION PL, PR
+* ..
+* .. Array Arguments ..
+* LOGICAL SELECT( * )
+* INTEGER IWORK( * )
+* DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+* $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ),
+* $ WORK( * ), Z( LDZ, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> GDTGSEN reorders the generalized real Schur decomposition of a real
+*> matrix pair (A, B) (in terms of an orthonormal equivalence trans-
+*> formation Q**T * (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.
+*>
+*> GDTGSEN also computes the generalized eigenvalues
+*>
+*> w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
+*>
+*> of the reordered matrix pair (A, B).
+*>
+*> Optionally, GDTGSEN 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.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] IJOB
+*> \verbatim
+*> IJOB is 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)
+*> \endverbatim
+*>
+*> \param[in] WANTQ
+*> \verbatim
+*> WANTQ is LOGICAL
+*> .TRUE. : update the left transformation matrix Q;
+*> .FALSE.: do not update Q.
+*> \endverbatim
+*>
+*> \param[in] WANTZ
+*> \verbatim
+*> WANTZ is LOGICAL
+*> .TRUE. : update the right transformation matrix Z;
+*> .FALSE.: do not update Z.
+*> \endverbatim
+*>
+*> \param[in] SELECT
+*> \verbatim
+*> SELECT is 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.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A and B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is 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.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is 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.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHAR
+*> \verbatim
+*> ALPHAR is DOUBLE PRECISION array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] ALPHAI
+*> \verbatim
+*> ALPHAI is DOUBLE PRECISION array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*> BETA is 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.
+*> \endverbatim
+*>
+*> \param[in,out] Q
+*> \verbatim
+*> Q is 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.
+*> \endverbatim
+*>
+*> \param[in] LDQ
+*> \verbatim
+*> LDQ is INTEGER
+*> The leading dimension of the array Q. LDQ >= 1;
+*> and if WANTQ = .TRUE., LDQ >= N.
+*> \endverbatim
+*>
+*> \param[in,out] Z
+*> \verbatim
+*> Z is 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.
+*> \endverbatim
+*>
+*> \param[in] LDZ
+*> \verbatim
+*> LDZ is INTEGER
+*> The leading dimension of the array Z. LDZ >= 1;
+*> If WANTZ = .TRUE., LDZ >= N.
+*> \endverbatim
+*>
+*> \param[out] M
+*> \verbatim
+*> M is INTEGER
+*> The dimension of the specified pair of left and right eigen-
+*> spaces (deflating subspaces). 0 <= M <= N.
+*> \endverbatim
+*>
+*> \param[out] PL
+*> \verbatim
+*> PL is DOUBLE PRECISION
+*> \endverbatim
+
+*> \param[out] PR
+*> \verbatim
+*> PR is 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.
+*> \endverbatim
+*>
+*> \param[out] DIF
+*> \verbatim
+*> DIF is 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.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array,
+*> dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is 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.
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
+*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
+*> \endverbatim
+*>
+*> \param[in] LIWORK
+*> \verbatim
+*> LIWORK is 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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is 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.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date June 2016
+*
+*> \ingroup doubleOTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> GDTGSEN 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**T*(A, B)*W = (A11 A12) (B11 B12) n1
+*> ( 0 A22),( 0 B22) n2
+*> n1 n2 n1 n2
+*>
+*> where N = n1+n2 and U**T 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**T, then the
+*> reordered generalized real Schur form of (C, D) is given by
+*>
+*> (C, D) = (Q*U)*(U**T*(A, B)*W)*(Z*W)**T,
+*>
+*> 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**T, In1) ]
+*> [ kron(In2, B11) -kron(B22**T, In1) ].
+*>
+*> Here, Inx is the identity matrix of size nx and A22**T 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.
+*> \endverbatim
+*
+*> \par Contributors:
+* ==================
+*>
+*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
+*> Umea University, S-901 87 Umea, Sweden.
+*
+*> \par References:
+* ================
+*>
+*> \verbatim
+*>
+*> [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.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE GDTGSEN( 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 )
+*
+* -- LAPACK computational routine (version 3.6.1) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* June 2016
+*
+* .. Scalar Arguments ..
+ LOGICAL WANTQ, WANTZ
+ INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
+ $ M, N
+ DOUBLE PRECISION PL, PR
+* ..
+* .. Array Arguments ..
+ LOGICAL SELECT( * )
+ INTEGER IWORK( * )
+ DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+ $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ),
+ $ WORK( * ), Z( LDZ, * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ INTEGER IDIFJB
+ PARAMETER ( IDIFJB = 3 )
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2,
+ $ WANTP
+ INTEGER I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN,
+ $ MN2, N1, N2
+ DOUBLE PRECISION DSCALE, DSUM, EPS, RDSCAL, SMLNUM
+* ..
+* .. Local Arrays ..
+ INTEGER ISAVE( 3 )
+* ..
+* .. External Subroutines ..
+ EXTERNAL DLACN2, DLACPY, DLAG2, DLASSQ, GDTGEXC, DTGSYL,
+ $ XERBLA
+* ..
+* .. External Functions ..
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL DLAMCH
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX, SIGN, SQRT
+* ..
+* .. Executable Statements ..
+*
+* Decode and test the input parameters
+*
+ INFO = 0
+ LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
+*
+ IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -5
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
+ INFO = -14
+ ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
+ INFO = -16
+ END IF
+*
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'GDTGSEN', -INFO )
+ RETURN
+ END IF
+*
+* Get machine constants
+*
+ EPS = DLAMCH( 'P' )
+ SMLNUM = DLAMCH( 'S' ) / EPS
+ IERR = 0
+*
+ WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
+ WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
+ WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
+ WANTD = WANTD1 .OR. WANTD2
+*
+* Set M to the dimension of the specified pair of deflating
+* subspaces.
+*
+ M = 0
+ PAIR = .FALSE.
+ IF( .NOT.LQUERY .OR. IJOB.NE.0 ) THEN
+ DO 10 K = 1, N
+ IF( PAIR ) THEN
+ PAIR = .FALSE.
+ ELSE
+ IF( K.LT.N ) THEN
+ IF( A( K+1, K ).EQ.ZERO ) THEN
+ IF( SELECT( K ) )
+ $ M = M + 1
+ ELSE
+ PAIR = .TRUE.
+ IF( SELECT( K ) .OR. SELECT( K+1 ) )
+ $ M = M + 2
+ END IF
+ ELSE
+ IF( SELECT( N ) )
+ $ M = M + 1
+ END IF
+ END IF
+ 10 CONTINUE
+ END IF
+*
+ IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
+ LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) )
+ LIWMIN = MAX( 1, N+6 )
+ ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
+ LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) )
+ LIWMIN = MAX( 1, 2*M*( N-M ), N+6 )
+ ELSE
+ LWMIN = MAX( 1, 4*N+16 )
+ LIWMIN = 1
+ END IF
+*
+ WORK( 1 ) = LWMIN
+ IWORK( 1 ) = LIWMIN
+*
+ IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
+ INFO = -22
+ ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
+ INFO = -24
+ END IF
+*
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'GDTGSEN', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF( M.EQ.N .OR. M.EQ.0 ) THEN
+ IF( WANTP ) THEN
+ PL = ONE
+ PR = ONE
+ END IF
+ IF( WANTD ) THEN
+ DSCALE = ZERO
+ DSUM = ONE
+ DO 20 I = 1, N
+ CALL DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
+ CALL DLASSQ( N, B( 1, I ), 1, DSCALE, DSUM )
+ 20 CONTINUE
+ DIF( 1 ) = DSCALE*SQRT( DSUM )
+ DIF( 2 ) = DIF( 1 )
+ END IF
+ GO TO 60
+ END IF
+*
+* Collect the selected blocks at the top-left corner of (A, B).
+*
+ KS = 0
+ PAIR = .FALSE.
+ DO 30 K = 1, N
+ IF( PAIR ) THEN
+ PAIR = .FALSE.
+ ELSE
+*
+ SWAP = SELECT( K )
+ IF( K.LT.N ) THEN
+ IF( A( K+1, K ).NE.ZERO ) THEN
+ PAIR = .TRUE.
+ SWAP = SWAP .OR. SELECT( K+1 )
+ END IF
+ END IF
+*
+ IF( SWAP ) THEN
+ KS = KS + 1
+*
+* Swap the K-th block to position KS.
+* Perform the reordering of diagonal blocks in (A, B)
+* by orthogonal transformation matrices and update
+* Q and Z accordingly (if requested):
+*
+ KK = K
+ IF( K.NE.KS )
+ $ CALL GDTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
+ $ Z, LDZ, KK, KS, WORK, LWORK, IERR )
+*
+ IF( IERR.GT.0 ) THEN
+*
+* Swap is rejected: exit.
+*
+ INFO = 1
+ IF( WANTP ) THEN
+ PL = ZERO
+ PR = ZERO
+ END IF
+ IF( WANTD ) THEN
+ DIF( 1 ) = ZERO
+ DIF( 2 ) = ZERO
+ END IF
+ GO TO 60
+ END IF
+*
+ IF( PAIR )
+ $ KS = KS + 1
+ END IF
+ END IF
+ 30 CONTINUE
+ IF( WANTP ) THEN
+*
+* Solve generalized Sylvester equation for R and L
+* and compute PL and PR.
+*
+ N1 = M
+ N2 = N - M
+ I = N1 + 1
+ IJB = 0
+ CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
+ CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
+ $ N1 )
+ CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
+ $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
+ $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
+ $ LWORK-2*N1*N2, IWORK, IERR )
+*
+* Estimate the reciprocal of norms of "projections" onto left
+* and right eigenspaces.
+*
+ RDSCAL = ZERO
+ DSUM = ONE
+ CALL DLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM )
+ PL = RDSCAL*SQRT( DSUM )
+ IF( PL.EQ.ZERO ) THEN
+ PL = ONE
+ ELSE
+ PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
+ END IF
+ RDSCAL = ZERO
+ DSUM = ONE
+ CALL DLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
+ PR = RDSCAL*SQRT( DSUM )
+ IF( PR.EQ.ZERO ) THEN
+ PR = ONE
+ ELSE
+ PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
+ END IF
+ END IF
+*
+ IF( WANTD ) THEN
+*
+* Compute estimates of Difu and Difl.
+*
+ IF( WANTD1 ) THEN
+ N1 = M
+ N2 = N - M
+ I = N1 + 1
+ IJB = IDIFJB
+*
+* Frobenius norm-based Difu-estimate.
+*
+ CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
+ $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
+ $ N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ),
+ $ LWORK-2*N1*N2, IWORK, IERR )
+*
+* Frobenius norm-based Difl-estimate.
+*
+ CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
+ $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
+ $ N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ),
+ $ LWORK-2*N1*N2, IWORK, IERR )
+ ELSE
+*
+*
+* Compute 1-norm-based estimates of Difu and Difl using
+* reversed communication with DLACN2. In each step a
+* generalized Sylvester equation or a transposed variant
+* is solved.
+*
+ KASE = 0
+ N1 = M
+ N2 = N - M
+ I = N1 + 1
+ IJB = 0
+ MN2 = 2*N1*N2
+*
+* 1-norm-based estimate of Difu.
+*
+ 40 CONTINUE
+ CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 1 ),
+ $ KASE, ISAVE )
+ IF( KASE.NE.0 ) THEN
+ IF( KASE.EQ.1 ) THEN
+*
+* Solve generalized Sylvester equation.
+*
+ CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
+ $ WORK, N1, B, LDB, B( I, I ), LDB,
+ $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
+ $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
+ $ IERR )
+ ELSE
+*
+* Solve the transposed variant.
+*
+ CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA,
+ $ WORK, N1, B, LDB, B( I, I ), LDB,
+ $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
+ $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
+ $ IERR )
+ END IF
+ GO TO 40
+ END IF
+ DIF( 1 ) = DSCALE / DIF( 1 )
+*
+* 1-norm-based estimate of Difl.
+*
+ 50 CONTINUE
+ CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 2 ),
+ $ KASE, ISAVE )
+ IF( KASE.NE.0 ) THEN
+ IF( KASE.EQ.1 ) THEN
+*
+* Solve generalized Sylvester equation.
+*
+ CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
+ $ WORK, N2, B( I, I ), LDB, B, LDB,
+ $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
+ $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
+ $ IERR )
+ ELSE
+*
+* Solve the transposed variant.
+*
+ CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA,
+ $ WORK, N2, B( I, I ), LDB, B, LDB,
+ $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
+ $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
+ $ IERR )
+ END IF
+ GO TO 50
+ END IF
+ DIF( 2 ) = DSCALE / DIF( 2 )
+*
+ END IF
+ END IF
+*
+ 60 CONTINUE
+*
+* Compute generalized eigenvalues of reordered pair (A, B) and
+* normalize the generalized Schur form.
+*
+ PAIR = .FALSE.
+ DO 80 K = 1, N
+ IF( PAIR ) THEN
+ PAIR = .FALSE.
+ ELSE
+*
+ IF( K.LT.N ) THEN
+ IF( A( K+1, K ).NE.ZERO ) THEN
+ PAIR = .TRUE.
+ END IF
+ END IF
+*
+ IF( PAIR ) THEN
+*
+* Compute the eigenvalue(s) at position K.
+*
+ WORK( 1 ) = A( K, K )
+ WORK( 2 ) = A( K+1, K )
+ WORK( 3 ) = A( K, K+1 )
+ WORK( 4 ) = A( K+1, K+1 )
+ WORK( 5 ) = B( K, K )
+ WORK( 6 ) = B( K+1, K )
+ WORK( 7 ) = B( K, K+1 )
+ WORK( 8 ) = B( K+1, K+1 )
+ CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ),
+ $ BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ),
+ $ ALPHAI( K ) )
+ ALPHAI( K+1 ) = -ALPHAI( K )
+*
+ ELSE
+*
+ IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN
+*
+* If B(K,K) is negative, make it positive
+*
+ DO 70 I = 1, N
+ A( K, I ) = -A( K, I )
+ B( K, I ) = -B( K, I )
+ IF( WANTQ ) Q( I, K ) = -Q( I, K )
+ 70 CONTINUE
+ END IF
+*
+ ALPHAR( K ) = A( K, K )
+ ALPHAI( K ) = ZERO
+ BETA( K ) = B( K, K )
+*
+ END IF
+ END IF
+ 80 CONTINUE
+*
+ WORK( 1 ) = LWMIN
+ IWORK( 1 ) = LIWMIN
+*
+ RETURN
+*
+* End of GDTGSEN
+*
+ END
diff --git a/src/PreviewControl/OptimalControllerSolver.cpp b/src/PreviewControl/OptimalControllerSolver.cpp
index 2e22aa4b..d4da3c8a 100644
--- a/src/PreviewControl/OptimalControllerSolver.cpp
+++ b/src/PreviewControl/OptimalControllerSolver.cpp
@@ -44,7 +44,7 @@ typedef int integer ;
extern "C" {
extern doublereal dlapy2_(doublereal *, doublereal *);
extern double dlamch_ (char *);
-extern /* Subroutine */ int dgges_(char *, char *, char *, L_fp, integer *
+extern /* Subroutine */ int gdgges_(char *, char *, char *, L_fp, integer *
, doublereal *, integer *, doublereal *, integer *, integer *,
doublereal *, doublereal *, doublereal *, doublereal *, integer *,
doublereal *, integer *, doublereal *, integer *, logical *,
@@ -163,7 +163,7 @@ bool OptimalControllerSolver::GeneralizedSchur(MAL_MATRIX( &A,double),
bwork[i] =0;
A = MAL_RET_TRANSPOSE(A);
B = MAL_RET_TRANSPOSE(B);
- dgges_ (lV, lV,
+ gdgges_ (lV, lV,
lS,
(logical (*)(...))sb02ox,
&n,
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 630403bb..c9959a86 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -75,6 +75,8 @@ ADD_EXECUTABLE(TestRiccatiEquation
)
# Need this dependency for the use of the fortran functions (it seems)
PKG_CONFIG_USE_DEPENDENCY(TestRiccatiEquation jrl-mal)
+TARGET_LINK_LIBRARIES(TestRiccatiEquation sublapack)
+
# Add test on the ricatti equation
ADD_TEST(TestRiccatiEquation TestRiccatiEquation)
--
GitLab