--- a
+++ b/src/SB03SY.f
@@ -0,0 +1,451 @@
+      SUBROUTINE SB03SY( JOB, TRANA, LYAPUN, N, T, LDT, U, LDU, XA,
+     $                   LDXA, SEPD, THNORM, IWORK, DWORK, LDWORK,
+     $                   INFO )
+C
+C     SLICOT RELEASE 5.0.
+C
+C     Copyright (c) 2002-2009 NICONET e.V.
+C
+C     This program is free software: you can redistribute it and/or
+C     modify it under the terms of the GNU General Public License as
+C     published by the Free Software Foundation, either version 2 of
+C     the License, or (at your option) any later version.
+C
+C     This program is distributed in the hope that it will be useful,
+C     but WITHOUT ANY WARRANTY; without even the implied warranty of
+C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+C     GNU General Public License for more details.
+C
+C     You should have received a copy of the GNU General Public License
+C     along with this program.  If not, see
+C     <http://www.gnu.org/licenses/>.
+C
+C     PURPOSE
+C
+C     To estimate the "separation" between the matrices op(A) and
+C     op(A)',
+C
+C     sepd(op(A),op(A)') = min norm(op(A)'*X*op(A) - X)/norm(X)
+C                        = 1 / norm(inv(Omega))
+C
+C     and/or the 1-norm of Theta, where op(A) = A or A' (A**T), and
+C     Omega and Theta are linear operators associated to the real
+C     discrete-time Lyapunov matrix equation
+C
+C            op(A)'*X*op(A) - X = C,
+C
+C     defined by
+C
+C     Omega(W) = op(A)'*W*op(A) - W,
+C     Theta(W) = inv(Omega(op(W)'*X*op(A) + op(A)'*X*op(W))).
+C
+C     The 1-norm condition estimators are used.
+C
+C     ARGUMENTS
+C
+C     Mode Parameters
+C
+C     JOB     CHARACTER*1
+C             Specifies the computation to be performed, as follows:
+C             = 'S':  Compute the separation only;
+C             = 'T':  Compute the norm of Theta only;
+C             = 'B':  Compute both the separation and the norm of Theta.
+C
+C     TRANA   CHARACTER*1
+C             Specifies the form of op(A) to be used, as follows:
+C             = 'N':  op(A) = A    (No transpose);
+C             = 'T':  op(A) = A**T (Transpose);
+C             = 'C':  op(A) = A**T (Conjugate transpose = Transpose).
+C
+C     LYAPUN  CHARACTER*1
+C             Specifies whether or not the original Lyapunov equations
+C             should be solved, as follows:
+C             = 'O':  Solve the original Lyapunov equations, updating
+C                     the right-hand sides and solutions with the
+C                     matrix U, e.g., X <-- U'*X*U;
+C             = 'R':  Solve reduced Lyapunov equations only, without
+C                     updating the right-hand sides and solutions.
+C
+C     Input/Output Parameters
+C
+C     N       (input) INTEGER
+C             The order of the matrices A and X.  N >= 0.
+C
+C     T       (input) DOUBLE PRECISION array, dimension (LDT,N)
+C             The leading N-by-N upper Hessenberg part of this array
+C             must contain the upper quasi-triangular matrix T in Schur
+C             canonical form from a Schur factorization of A.
+C
+C     LDT     INTEGER
+C             The leading dimension of array T.  LDT >= MAX(1,N).
+C
+C     U       (input) DOUBLE PRECISION array, dimension (LDU,N)
+C             The leading N-by-N part of this array must contain the
+C             orthogonal matrix U from a real Schur factorization of A.
+C             If LYAPUN = 'R', the array U is not referenced.
+C
+C     LDU     INTEGER
+C             The leading dimension of array U.
+C             LDU >= 1,        if LYAPUN = 'R';
+C             LDU >= MAX(1,N), if LYAPUN = 'O'.
+C
+C     XA      (input) DOUBLE PRECISION array, dimension (LDXA,N)
+C             The leading N-by-N part of this array must contain the
+C             matrix product X*op(A), if LYAPUN = 'O', or U'*X*U*op(T),
+C             if LYAPUN = 'R', in the Lyapunov equation.
+C             If JOB = 'S', the array XA is not referenced.
+C
+C     LDXA    INTEGER
+C             The leading dimension of array XA.
+C             LDXA >= 1,        if JOB = 'S';
+C             LDXA >= MAX(1,N), if JOB = 'T' or 'B'.
+C
+C     SEPD    (output) DOUBLE PRECISION
+C             If JOB = 'S' or JOB = 'B', and INFO >= 0, SEPD contains
+C             the estimated quantity sepd(op(A),op(A)').
+C             If JOB = 'T' or N = 0, SEPD is not referenced.
+C
+C     THNORM  (output) DOUBLE PRECISION
+C             If JOB = 'T' or JOB = 'B', and INFO >= 0, THNORM contains
+C             the estimated 1-norm of operator Theta.
+C             If JOB = 'S' or N = 0, THNORM is not referenced.
+C
+C     Workspace
+C
+C     IWORK   INTEGER array, dimension (N*N)
+C
+C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
+C
+C     LDWORK  INTEGER
+C             The length of the array DWORK.
+C             LDWORK >= 0,            if N = 0;
+C             LDWORK >= MAX(3,2*N*N), if N > 0.
+C
+C     Error Indicator
+C
+C     INFO    INTEGER
+C             = 0:  successful exit;
+C             < 0:  if INFO = -i, the i-th argument had an illegal
+C                   value;
+C             = N+1:  if T has (almost) reciprocal eigenvalues;
+C                   perturbed values were used to solve Lyapunov
+C                   equations (but the matrix T is unchanged).
+C
+C     METHOD
+C
+C     SEPD is defined as
+C
+C            sepd( op(A), op(A)' ) = sigma_min( K )
+C
+C     where sigma_min(K) is the smallest singular value of the
+C     N*N-by-N*N matrix
+C
+C        K = kprod( op(A)', op(A)' ) - I(N**2).
+C
+C     I(N**2) is an N*N-by-N*N identity matrix, and kprod denotes the
+C     Kronecker product. The routine estimates sigma_min(K) by the
+C     reciprocal of an estimate of the 1-norm of inverse(K), computed as
+C     suggested in [1]. This involves the solution of several discrete-
+C     time Lyapunov equations, either direct or transposed. The true
+C     reciprocal 1-norm of inverse(K) cannot differ from sigma_min(K) by
+C     more than a factor of N.
+C     The 1-norm of Theta is estimated similarly.
+C
+C     REFERENCES
+C
+C     [1] Higham, N.J.
+C         FORTRAN codes for estimating the one-norm of a real or
+C         complex matrix, with applications to condition estimation.
+C         ACM Trans. Math. Softw., 14, pp. 381-396, 1988.
+C
+C     NUMERICAL ASPECTS
+C                               3
+C     The algorithm requires 0(N ) operations.
+C
+C     FURTHER COMMENTS
+C
+C     When SEPD is zero, the routine returns immediately, with THNORM
+C     (if requested) not set. In this case, the equation is singular.
+C     The option LYAPUN = 'R' may occasionally produce slightly worse
+C     or better estimates, and it is much faster than the option 'O'.
+C
+C     CONTRIBUTOR
+C
+C     V. Sima, Research Institute for Informatics, Bucharest, Romania,
+C     Oct. 1998. Partly based on DDLSVX (and then SB03SD) by P. Petkov,
+C     Tech. University of Sofia, March 1998 (and December 1998).
+C
+C     REVISIONS
+C
+C     February 6, 1999, V. Sima, Katholieke Univ. Leuven, Belgium.
+C     V. Sima, Research Institute for Informatics, Bucharest, Oct. 2004.
+C
+C     KEYWORDS
+C
+C     Lyapunov equation, orthogonal transformation, real Schur form.
+C
+C     ******************************************************************
+C
+C     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE, HALF
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 )
+C     ..
+C     .. Scalar Arguments ..
+      CHARACTER          JOB, LYAPUN, TRANA
+      INTEGER            INFO, LDT, LDU, LDWORK, LDXA, N
+      DOUBLE PRECISION   SEPD, THNORM
+C     ..
+C     .. Array Arguments ..
+      INTEGER            IWORK( * )
+      DOUBLE PRECISION   DWORK( * ), T( LDT, * ), U( LDU, * ),
+     $                   XA( LDXA, * )
+C     ..
+C     .. Local Scalars ..
+      LOGICAL            NOTRNA, UPDATE, WANTS, WANTT
+      CHARACTER          TRANAT, UPLO
+      INTEGER            INFO2, ITMP, KASE, NN
+      DOUBLE PRECISION   BIGNUM, EST, SCALE
+C     ..
+C     .. External Functions ..
+      LOGICAL            LSAME
+      DOUBLE PRECISION   DLAMCH, DLANSY
+      EXTERNAL           DLAMCH, DLANSY, LSAME
+C     ..
+C     .. External Subroutines ..
+      EXTERNAL           DLACON, DLACPY, DSCAL, DSYR2K, MA02ED, MB01RU,
+     $                   SB03MX, XERBLA
+C     ..
+C     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+C     ..
+C     .. Executable Statements ..
+C
+C     Decode and Test input parameters.
+C
+      WANTS  = LSAME( JOB,    'S' )
+      WANTT  = LSAME( JOB,    'T' )
+      NOTRNA = LSAME( TRANA,  'N' )
+      UPDATE = LSAME( LYAPUN, 'O' )
+C
+      NN   = N*N
+      INFO = 0
+      IF( .NOT. ( WANTS .OR. WANTT .OR. LSAME( JOB, 'B' ) ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.( NOTRNA .OR. LSAME( TRANA, 'T' ) .OR.
+     $                            LSAME( TRANA, 'C' ) ) ) THEN
+         INFO = -2
+      ELSE IF( .NOT.( UPDATE .OR. LSAME( LYAPUN, 'R' ) ) ) THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+         INFO = -6
+      ELSE IF( LDU.LT.1 .OR. ( UPDATE .AND. LDU.LT.N ) ) THEN
+         INFO = -8
+      ELSE IF( LDXA.LT.1 .OR. ( .NOT.WANTS .AND. LDXA.LT.N ) ) THEN
+         INFO = -10
+      ELSE IF( LDWORK.LT.0 .OR.
+     $       ( LDWORK.LT.MAX( 3, 2*NN ) .AND. N.GT.0 ) ) THEN
+         INFO = -15
+      END IF
+C
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SB03SY', -INFO )
+         RETURN
+      END IF
+C
+C     Quick return if possible.
+C
+      IF( N.EQ.0 )
+     $   RETURN
+C
+      ITMP = NN + 1
+C
+      IF( NOTRNA ) THEN
+         TRANAT = 'T'
+      ELSE
+         TRANAT = 'N'
+      END IF
+C
+      IF( .NOT.WANTT ) THEN
+C
+C        Estimate sepd(op(A),op(A)').
+C        Workspace:  max(3,2*N*N).
+C
+         KASE = 0
+C
+C        REPEAT
+   10    CONTINUE
+         CALL DLACON( NN, DWORK( ITMP ), DWORK, IWORK, EST, KASE )
+         IF( KASE.NE.0 ) THEN
+C
+C           Select the triangular part of symmetric matrix to be used.
+C
+            IF( DLANSY( '1-norm', 'Upper', N, DWORK, N, DWORK( ITMP ) )
+     $          .GE.
+     $          DLANSY( '1-norm', 'Lower', N, DWORK, N, DWORK( ITMP ) )
+     $        ) THEN
+               UPLO = 'U'
+            ELSE
+               UPLO = 'L'
+            END IF
+C
+            IF( UPDATE ) THEN
+C
+C              Transform the right-hand side: RHS := U'*RHS*U.
+C
+               CALL MB01RU( UPLO, 'Transpose', N, N, ZERO, ONE, DWORK,
+     $                      N, U, LDU, DWORK, N, DWORK( ITMP ), NN,
+     $                      INFO2 )
+               CALL DSCAL( N, HALF, DWORK, N+1 )
+            END IF
+            CALL MA02ED( UPLO, N, DWORK, N )
+C
+            IF( KASE.EQ.1 ) THEN
+C
+C              Solve op(T)'*Y*op(T) - Y = scale*RHS.
+C
+               CALL SB03MX( TRANA, N, T, LDT, DWORK, N, SCALE,
+     $                      DWORK( ITMP ), INFO2 )
+            ELSE
+C
+C              Solve op(T)*W*op(T)' - W = scale*RHS.
+C
+               CALL SB03MX( TRANAT, N, T, LDT, DWORK, N, SCALE,
+     $                      DWORK( ITMP ), INFO2 )
+            END IF
+C
+            IF( INFO2.GT.0 )
+     $         INFO = N + 1
+C
+            IF( UPDATE ) THEN
+C
+C              Transform back to obtain the solution: Z := U*Z*U', with
+C              Z = Y or Z = W.
+C
+               CALL MB01RU( UPLO, 'No transpose', N, N, ZERO, ONE,
+     $                      DWORK, N, U, LDU, DWORK, N, DWORK( ITMP ),
+     $                      NN, INFO2 )
+               CALL DSCAL( N, HALF, DWORK, N+1 )
+C
+C              Fill in the remaining triangle of the symmetric matrix.
+C
+               CALL MA02ED( UPLO, N, DWORK, N )
+            END IF
+C
+            GO TO 10
+         END IF
+C        UNTIL KASE = 0
+C
+         IF( EST.GT.SCALE ) THEN
+            SEPD = SCALE / EST
+         ELSE
+            BIGNUM = ONE / DLAMCH( 'Safe minimum' )
+            IF( SCALE.LT.EST*BIGNUM ) THEN
+               SEPD = SCALE / EST
+            ELSE
+               SEPD = BIGNUM
+            END IF
+         END IF
+C
+C        Return if the equation is singular.
+C
+         IF( SEPD.EQ.ZERO )
+     $      RETURN
+      END IF
+C
+      IF( .NOT.WANTS ) THEN
+C
+C        Estimate norm(Theta).
+C        Workspace:  max(3,2*N*N).
+C
+         KASE = 0
+C
+C        REPEAT
+   20    CONTINUE
+         CALL DLACON( NN, DWORK( ITMP ), DWORK, IWORK, EST, KASE )
+         IF( KASE.NE.0 ) THEN
+C
+C           Select the triangular part of symmetric matrix to be used.
+C
+            IF( DLANSY( '1-norm', 'Upper', N, DWORK, N, DWORK( ITMP ) )
+     $          .GE.
+     $          DLANSY( '1-norm', 'Lower', N, DWORK, N, DWORK( ITMP ) )
+     $        ) THEN
+               UPLO = 'U'
+            ELSE
+               UPLO = 'L'
+            END IF
+C
+C           Fill in the remaining triangle of the symmetric matrix.
+C
+            CALL MA02ED( UPLO, N, DWORK, N )
+C
+C           Compute RHS = op(W)'*X*op(A) + op(A)'*X*op(W).
+C
+            CALL DSYR2K( UPLO, TRANAT, N, N, ONE, DWORK, N, XA, LDXA,
+     $                   ZERO, DWORK( ITMP ), N )
+            CALL DLACPY( UPLO, N, N, DWORK( ITMP ), N, DWORK, N )
+C
+            IF( UPDATE ) THEN
+C
+C              Transform the right-hand side: RHS := U'*RHS*U.
+C
+               CALL MB01RU( UPLO, 'Transpose', N, N, ZERO, ONE, DWORK,
+     $                      N, U, LDU, DWORK, N, DWORK( ITMP ), NN,
+     $                      INFO2 )
+               CALL DSCAL( N, HALF, DWORK, N+1 )
+            END IF
+            CALL MA02ED( UPLO, N, DWORK, N )
+C
+            IF( KASE.EQ.1 ) THEN
+C
+C              Solve op(T)'*Y*op(T) - Y = scale*RHS.
+C
+               CALL SB03MX( TRANA, N, T, LDT, DWORK, N, SCALE,
+     $                      DWORK( ITMP ), INFO2 )
+            ELSE
+C
+C              Solve op(T)*W*op(T)' - W = scale*RHS.
+C
+               CALL SB03MX( TRANAT, N, T, LDT, DWORK, N, SCALE,
+     $                      DWORK( ITMP ), INFO2 )
+            END IF
+C
+            IF( INFO2.GT.0 )
+     $         INFO = N + 1
+C
+            IF( UPDATE ) THEN
+C
+C              Transform back to obtain the solution: Z := U*Z*U', with
+C              Z = Y or Z = W.
+C
+               CALL MB01RU( UPLO, 'No transpose', N, N, ZERO, ONE,
+     $                      DWORK, N, U, LDU, DWORK, N, DWORK( ITMP ),
+     $                      NN, INFO2 )
+               CALL DSCAL( N, HALF, DWORK, N+1 )
+C
+C              Fill in the remaining triangle of the symmetric matrix.
+C
+               CALL MA02ED( UPLO, N, DWORK, N )
+            END IF
+C
+            GO TO 20
+         END IF
+C        UNTIL KASE = 0
+C
+         IF( EST.LT.SCALE ) THEN
+            THNORM = EST / SCALE
+         ELSE
+            BIGNUM = ONE / DLAMCH( 'Safe minimum' )
+            IF( EST.LT.SCALE*BIGNUM ) THEN
+               THNORM = EST / SCALE
+            ELSE
+               THNORM = BIGNUM
+            END IF
+         END IF
+      END IF
+C
+      RETURN
+C *** Last line of SB03SY ***
+      END