--- a/extra/control-devel/devel/makefile_modred.m	Wed Feb 22 15:13:45 2012 +0000
+++ b/extra/control-devel/devel/makefile_modred.m	Wed Feb 22 16:37:14 2012 +0000
@@ -13,36 +13,17 @@
 cd (srcdir);
 
 mkoctfile slab09hd.cc \
-          AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \
-          AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \
-          SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \
-          SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \
-          SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \
-          SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \
+          slicotlibrary.a \
           "$(mkoctfile -p LAPACK_LIBS)" \
           "$(mkoctfile -p BLAS_LIBS)"
 
 mkoctfile slab09id.cc \
-          AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \
-          AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \
-          MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \
-          AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \
-          MA02BD.f MB03OY.f MB03QX.f MB01PD.f select.f \
-          MB01YD.f MB04NY.f MB01ZD.f SB03OT.f MB04OX.f \
-          MB04OY.f MB03QD.f SB03OY.f MB03QY.f MB01QD.f \
-          SB03OR.f SB03OV.f SB04PX.f \
+          slicotlibrary.a \
           "$(mkoctfile -p LAPACK_LIBS)" \
           "$(mkoctfile -p BLAS_LIBS)"
 
 mkoctfile slab09jd.cc \
-          AB09JD.f TB01ID.f TB01KD.f AB07ND.f AB09JV.f \
-          AB09JW.f AB09CX.f AG07BD.f AB08MD.f AB04MD.f \
-          TB01LD.f delctg.f SB04PY.f AB09AX.f AB08NX.f \
-          MB01SD.f AB09JX.f MA02AD.f TB01WD.f MB03OY.f \
-          MB03PY.f MA02DD.f MB03UD.f MB03QX.f select.f \
-          SB04PX.f SB03OU.f MB03QD.f MB03QY.f SB03OT.f \
-          MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \
-          MB04OY.f SB03OV.f \
+          slicotlibrary.a \
           "$(mkoctfile -p LAPACK_LIBS)" \
           "$(mkoctfile -p BLAS_LIBS)"
 

--- a/extra/control-devel/src/AB04MD.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,345 +0,0 @@
-      SUBROUTINE AB04MD( TYPE, N, M, P, ALPHA, BETA, A, LDA, B, LDB, C,
-     $                   LDC, D, LDD, 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 perform a transformation on the parameters (A,B,C,D) of a
-C     system, which is equivalent to a bilinear transformation of the
-C     corresponding transfer function matrix.
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     TYPE    CHARACTER*1
-C             Indicates the type of the original system and the
-C             transformation to be performed as follows:
-C             = 'D':  discrete-time   -> continuous-time;
-C             = 'C':  continuous-time -> discrete-time.
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The order of the state matrix A.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of system outputs.  P >= 0.
-C
-C     ALPHA,  (input) DOUBLE PRECISION
-C     BETA    Parameters specifying the bilinear transformation.
-C             Recommended values for stable systems: ALPHA = 1,
-C             BETA = 1.  ALPHA <> 0, BETA <> 0.
-C
-C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C             On entry, the leading N-by-N part of this array must
-C             contain the state matrix A of the original system.
-C             On exit, the leading N-by-N part of this array contains
-C                              _
-C             the state matrix A of the transformed system.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N).
-C
-C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C             On entry, the leading N-by-M part of this array must
-C             contain the input matrix B of the original system.
-C             On exit, the leading N-by-M part of this array contains
-C                              _
-C             the input matrix B of the transformed system.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N).
-C
-C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C             On entry, the leading P-by-N part of this array must
-C             contain the output matrix C of the original system.
-C             On exit, the leading P-by-N part of this array contains
-C                               _
-C             the output matrix C of the transformed system.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.  LDC >= MAX(1,P).
-C
-C     D       (input/output) DOUBLE PRECISION array, dimension (LDD,M)
-C             On entry, the leading P-by-M part of this array must
-C             contain the input/output matrix D for the original system.
-C             On exit, the leading P-by-M part of this array contains
-C                                     _
-C             the input/output matrix D of the transformed system.
-C
-C     LDD     INTEGER
-C             The leading dimension of array D.  LDD >= MAX(1,P).
-C
-C     Workspace
-C
-C     IWORK   INTEGER array, dimension (N)
-C
-C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
-C             On exit, if INFO = 0, DWORK(1) returns the optimal value
-C             of LDWORK.
-C
-C     LDWORK  INTEGER
-C             The length of the array DWORK.  LDWORK >= MAX(1,N).
-C             For optimum performance LDWORK >= MAX(1,N*NB), where NB
-C             is the optimal blocksize.
-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             = 1:  if the matrix (ALPHA*I + A) is exactly singular;
-C             = 2:  if the matrix  (BETA*I - A) is exactly singular.
-C
-C     METHOD
-C
-C     The parameters of the discrete-time system are transformed into
-C     the parameters of the continuous-time system (TYPE = 'D'), or
-C     vice-versa (TYPE = 'C') by the transformation:
-C
-C     1.  Discrete -> continuous
-C         _                     -1
-C         A = beta*(alpha*I + A)  * (A - alpha*I)
-C         _                                     -1
-C         B = sqrt(2*alpha*beta) * (alpha*I + A)  * B
-C         _                                         -1
-C         C = sqrt(2*alpha*beta) * C * (alpha*I + A)
-C         _                        -1
-C         D = D - C * (alpha*I + A)  * B
-C
-C     which is equivalent to the bilinear transformation
-C
-C                       z - alpha
-C         z -> s = beta ---------  .
-C                       z + alpha
-C
-C     of one transfer matrix onto the other.
-C
-C     2.  Continuous -> discrete
-C         _                     -1
-C         A = alpha*(beta*I - A)  * (beta*I + A)
-C         _                                    -1
-C         B = sqrt(2*alpha*beta) * (beta*I - A)  * B
-C         _                                        -1
-C         C = sqrt(2*alpha*beta) * C * (beta*I - A)
-C         _                       -1
-C         D = D + C * (beta*I - A)  * B
-C
-C     which is equivalent to the bilinear transformation
-C
-C                      beta + s
-C       s -> z = alpha -------- .
-C                      beta - s
-C
-C     of one transfer matrix onto the other.
-C
-C     REFERENCES
-C
-C     [1] Al-Saggaf, U.M. and Franklin, G.F.
-C         Model reduction via balanced realizations: a extension and
-C         frequency weighting techniques.
-C         IEEE Trans. Autom. Contr., AC-33, pp. 687-692, 1988.
-C
-C     NUMERICAL ASPECTS
-C                                                      3
-C     The time taken is approximately proportional to N .
-C     The accuracy depends mainly on the condition number of the matrix
-C     to be inverted.
-C
-C     CONTRIBUTORS
-C
-C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, and
-C                  A. Varga, German Aerospace Research Establishment,
-C                  Oberpfaffenhofen, Germany, Nov. 1996.
-C     Supersedes Release 2.0 routine AB04AD by W. van der Linden, and
-C     A.J. Geurts, Technische Hogeschool Eindhoven, Holland.
-C
-C     REVISIONS
-C
-C     -
-C
-C     KEYWORDS
-C
-C     Bilinear transformation, continuous-time system, discrete-time
-C     system, state-space model.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ZERO, ONE, TWO
-      PARAMETER         ( ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0 )
-C     .. Scalar Arguments ..
-      CHARACTER         TYPE
-      INTEGER           INFO, LDA, LDB, LDC, LDD, LDWORK, M, N, P
-      DOUBLE PRECISION  ALPHA, BETA
-C     .. Array Arguments ..
-      INTEGER           IWORK(*)
-      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), DWORK(*)
-C     .. Local Scalars ..
-      LOGICAL           LTYPE
-      INTEGER           I, IP
-      DOUBLE PRECISION  AB2, PALPHA, PBETA, SQRAB2
-C     .. External Functions ..
-      LOGICAL           LSAME
-      EXTERNAL          LSAME
-C     .. External Subroutines ..
-      EXTERNAL          DGEMM, DGETRF, DGETRS, DGETRI, DLASCL, DSCAL,
-     $                  DSWAP, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         ABS, MAX, SIGN, SQRT
-C     .. Executable Statements ..
-C
-      INFO = 0
-      LTYPE = LSAME( TYPE, 'D' )
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT.LTYPE .AND. .NOT.LSAME( TYPE, 'C' ) ) THEN
-         INFO = -1
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -3
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -4
-      ELSE IF( ALPHA.EQ.ZERO ) THEN
-         INFO = -5
-      ELSE IF( BETA.EQ.ZERO ) THEN
-         INFO = -6
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -8
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -10
-      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
-         INFO = -12
-      ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
-         INFO = -14
-      ELSE IF( LDWORK.LT.MAX( 1, N ) ) THEN
-         INFO = -17
-      END IF
-C
-      IF ( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB04MD', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF ( MAX( N, M, P ).EQ.0 )
-     $   RETURN
-C
-C     (Note: Comments in the code beginning "Workspace:" describe the
-C     minimal amount of real workspace needed at that point in the
-C     code, as well as the preferred amount for good performance.
-C     NB refers to the optimal block size for the immediately
-C     following subroutine, as returned by ILAENV.)
-C
-      IF (LTYPE) THEN
-C
-C        Discrete-time to continuous-time with (ALPHA, BETA).
-C
-         PALPHA = ALPHA
-         PBETA = BETA
-      ELSE
-C
-C        Continuous-time to discrete-time with (ALPHA, BETA) is
-C        equivalent with discrete-time to continuous-time with
-C        (-BETA, -ALPHA), if B and C change the sign.
-C
-         PALPHA = -BETA
-         PBETA = -ALPHA
-      END IF
-C
-      AB2 = PALPHA*PBETA*TWO
-      SQRAB2 = SIGN( SQRT( ABS( AB2 ) ), PALPHA )
-C                          -1
-C     Compute (alpha*I + A)  .
-C
-      DO 10 I = 1, N
-         A(I,I)  =  A(I,I) + PALPHA
-   10 CONTINUE
-C
-      CALL DGETRF( N, N, A, LDA, IWORK, INFO )
-C
-      IF (INFO.NE.0) THEN
-C
-C        Error return.
-C
-         IF (LTYPE) THEN
-            INFO = 1
-         ELSE
-            INFO = 2
-         END IF
-         RETURN
-      END IF
-C                         -1
-C     Compute  (alpha*I+A)  *B.
-C
-      CALL DGETRS( 'No transpose', N, M, A, LDA, IWORK, B, LDB, INFO )
-C                               -1
-C     Compute  D - C*(alpha*I+A)  *B.
-C
-      CALL DGEMM( 'No transpose', 'No transpose', P, M, N, -ONE, C,
-     $            LDC, B, LDB, ONE, D, LDD )
-C
-C     Scale B by  sqrt(2*alpha*beta).
-C
-      CALL DLASCL( 'General', 0, 0, ONE, SQRAB2, N, M, B, LDB, INFO )
-C                                                -1
-C     Compute  sqrt(2*alpha*beta)*C*(alpha*I + A)  .
-C
-      CALL DTRSM( 'Right', 'Upper', 'No transpose', 'Non-unit', P, N,
-     $            SQRAB2, A, LDA, C, LDC )
-C
-      CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', P, N, ONE,
-     $            A, LDA, C, LDC )
-C
-C     Apply column interchanges to the solution matrix.
-C
-      DO 20 I = N-1, 1, -1
-         IP = IWORK(I)
-         IF ( IP.NE.I )
-     $      CALL DSWAP( P, C(1,I), 1, C(1,IP), 1 )
-  20  CONTINUE
-C                               -1
-C     Compute beta*(alpha*I + A)  *(A - alpha*I) as
-C                                        -1
-C     beta*I - 2*alpha*beta*(alpha*I + A)  .
-C
-C     Workspace: need N;  prefer N*NB.
-C
-      CALL DGETRI( N, A, LDA, IWORK, DWORK, LDWORK, INFO )
-C
-      DO 30 I = 1, N
-         CALL DSCAL(N, -AB2, A(1,I), 1)
-         A(I,I) = A(I,I) + PBETA
-   30 CONTINUE
-C
-      RETURN
-C *** Last line of AB04MD ***
-      END

--- a/extra/control-devel/src/AB05PD.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,385 +0,0 @@
-      SUBROUTINE AB05PD( OVER, N1, M, P, N2, ALPHA, A1, LDA1, B1, LDB1,
-     $                   C1, LDC1, D1, LDD1, A2, LDA2, B2, LDB2, C2,
-     $                   LDC2, D2, LDD2, N, A, LDA, B, LDB, C, LDC, D,
-     $                   LDD, 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 compute the state-space model G = (A,B,C,D) corresponding to
-C     the sum G = G1 + alpha*G2, where G1 = (A1,B1,C1,D1) and
-C     G2 = (A2,B2,C2,D2).  G, G1, and G2 are the transfer-function
-C     matrices of the corresponding state-space models.
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     OVER    CHARACTER*1
-C             Indicates whether the user wishes to overlap pairs of
-C             arrays, as follows:
-C             = 'N':  Do not overlap;
-C             = 'O':  Overlap pairs of arrays: A1 and A, B1 and B,
-C                     C1 and C, and D1 and D, i.e. the same name is
-C                     effectively used for each pair (for all pairs)
-C                     in the routine call.  In this case, setting
-C                     LDA1 = LDA, LDB1 = LDB, LDC1 = LDC, and LDD1 = LDD
-C                     will give maximum efficiency.
-C
-C     Input/Output Parameters
-C
-C     N1      (input) INTEGER
-C             The number of state variables in the first system, i.e.
-C             the order of the matrix A1, the number of rows of B1 and
-C             the number of columns of C1.  N1 >= 0.
-C
-C     M       (input) INTEGER
-C             The number of input variables of the two systems, i.e. the
-C             number of columns of matrices B1, D1, B2 and D2.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of output variables of the two systems, i.e.
-C             the number of rows of matrices C1, D1, C2 and D2.  P >= 0.
-C
-C     N2      (input) INTEGER
-C             The number of state variables in the second system, i.e.
-C             the order of the matrix A2, the number of rows of B2 and
-C             the number of columns of C2.  N2 >= 0.
-C
-C     ALPHA   (input) DOUBLE PRECISION
-C             The coefficient multiplying G2.
-C
-C     A1      (input) DOUBLE PRECISION array, dimension (LDA1,N1)
-C             The leading N1-by-N1 part of this array must contain the
-C             state transition matrix A1 for the first system.
-C
-C     LDA1    INTEGER
-C             The leading dimension of array A1.  LDA1 >= MAX(1,N1).
-C
-C     B1      (input) DOUBLE PRECISION array, dimension (LDB1,M)
-C             The leading N1-by-M part of this array must contain the
-C             input/state matrix B1 for the first system.
-C
-C     LDB1    INTEGER
-C             The leading dimension of array B1.  LDB1 >= MAX(1,N1).
-C
-C     C1      (input) DOUBLE PRECISION array, dimension (LDC1,N1)
-C             The leading P-by-N1 part of this array must contain the
-C             state/output matrix C1 for the first system.
-C
-C     LDC1    INTEGER
-C             The leading dimension of array C1.
-C             LDC1 >= MAX(1,P) if N1 > 0.
-C             LDC1 >= 1 if N1 = 0.
-C
-C     D1      (input) DOUBLE PRECISION array, dimension (LDD1,M)
-C             The leading P-by-M part of this array must contain the
-C             input/output matrix D1 for the first system.
-C
-C     LDD1    INTEGER
-C             The leading dimension of array D1.  LDD1 >= MAX(1,P).
-C
-C     A2      (input) DOUBLE PRECISION array, dimension (LDA2,N2)
-C             The leading N2-by-N2 part of this array must contain the
-C             state transition matrix A2 for the second system.
-C
-C     LDA2    INTEGER
-C             The leading dimension of array A2.  LDA2 >= MAX(1,N2).
-C
-C     B2      (input) DOUBLE PRECISION array, dimension (LDB2,M)
-C             The leading N2-by-M part of this array must contain the
-C             input/state matrix B2 for the second system.
-C
-C     LDB2    INTEGER
-C             The leading dimension of array B2.  LDB2 >= MAX(1,N2).
-C
-C     C2      (input) DOUBLE PRECISION array, dimension (LDC2,N2)
-C             The leading P-by-N2 part of this array must contain the
-C             state/output matrix C2 for the second system.
-C
-C     LDC2    INTEGER
-C             The leading dimension of array C2.
-C             LDC2 >= MAX(1,P) if N2 > 0.
-C             LDC2 >= 1 if N2 = 0.
-C
-C     D2      (input) DOUBLE PRECISION array, dimension (LDD2,M)
-C             The leading P-by-M part of this array must contain the
-C             input/output matrix D2 for the second system.
-C
-C     LDD2    INTEGER
-C             The leading dimension of array D2.  LDD2 >= MAX(1,P).
-C
-C     N       (output) INTEGER
-C             The number of state variables (N1 + N2) in the resulting
-C             system, i.e. the order of the matrix A, the number of rows
-C             of B and the number of columns of C.
-C
-C     A       (output) DOUBLE PRECISION array, dimension (LDA,N1+N2)
-C             The leading N-by-N part of this array contains the state
-C             transition matrix A for the resulting system.
-C             The array A can overlap A1 if OVER = 'O'.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N1+N2).
-C
-C     B       (output) DOUBLE PRECISION array, dimension (LDB,M)
-C             The leading N-by-M part of this array contains the
-C             input/state matrix B for the resulting system.
-C             The array B can overlap B1 if OVER = 'O'.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N1+N2).
-C
-C     C       (output) DOUBLE PRECISION array, dimension (LDC,N1+N2)
-C             The leading P-by-N part of this array contains the
-C             state/output matrix C for the resulting system.
-C             The array C can overlap C1 if OVER = 'O'.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.
-C             LDC >= MAX(1,P) if N1+N2 > 0.
-C             LDC >= 1 if N1+N2 = 0.
-C
-C     D       (output) DOUBLE PRECISION array, dimension (LDD,M)
-C             The leading P-by-M part of this array contains the
-C             input/output matrix D for the resulting system.
-C             The array D can overlap D1 if OVER = 'O'.
-C
-C     LDD     INTEGER
-C             The leading dimension of array D.  LDD >= MAX(1,P).
-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
-C     METHOD
-C
-C     The matrices of the resulting systems are determined as:
-C
-C           ( A1   0  )             ( B1 )
-C       A = (         ) ,       B = (    ) ,
-C           ( 0    A2 )             ( B2 )
-C
-C       C = ( C1  alpha*C2 ) ,  D = D1 + alpha*D2 .
-C
-C     REFERENCES
-C
-C     None
-C
-C     NUMERICAL ASPECTS
-C
-C     None
-C
-C     CONTRIBUTORS
-C
-C     A. Varga, German Aerospace Research Establishment,
-C     Oberpfaffenhofen, Germany, and V. Sima, Katholieke Univ. Leuven,
-C     Belgium, Nov. 1996.
-C
-C     REVISIONS
-C
-C     V. Sima, Research Institute for Informatics, Bucharest, July 2003,
-C     Feb. 2004.
-C
-C     KEYWORDS
-C
-C     Multivariable system, state-space model, state-space
-C     representation.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ZERO, ONE
-      PARAMETER         ( ZERO=0.0D0, ONE = 1.0D0 )
-C     .. Scalar Arguments ..
-      CHARACTER         OVER
-      INTEGER           INFO, LDA, LDA1, LDA2, LDB, LDB1, LDB2, LDC,
-     $                  LDC1, LDC2, LDD, LDD1, LDD2, M, N, N1, N2, P
-      DOUBLE PRECISION  ALPHA
-C     .. Array Arguments ..
-      DOUBLE PRECISION  A(LDA,*), A1(LDA1,*), A2(LDA2,*), B(LDB,*),
-     $                  B1(LDB1,*), B2(LDB2,*), C(LDC,*), C1(LDC1,*),
-     $                  C2(LDC2,*), D(LDD,*), D1(LDD1,*), D2(LDD2,*)
-C     .. Local Scalars ..
-      LOGICAL           LOVER
-      INTEGER           I, J, N1P1
-C     .. External Functions ..
-      LOGICAL           LSAME
-      EXTERNAL          LSAME
-C     .. External Subroutines ..
-      EXTERNAL          DAXPY, DLACPY, DLASCL, DLASET, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         MAX, MIN
-C     .. Executable Statements ..
-C
-      LOVER = LSAME( OVER, 'O' )
-      N = N1 + N2
-      INFO = 0
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT.LOVER .AND. .NOT.LSAME( OVER, 'N' ) ) THEN
-         INFO = -1
-      ELSE IF( N1.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -3
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -4
-      ELSE IF( N2.LT.0 ) THEN
-         INFO = -5
-      ELSE IF( LDA1.LT.MAX( 1, N1 ) ) THEN
-         INFO = -8
-      ELSE IF( LDB1.LT.MAX( 1, N1 ) ) THEN
-         INFO = -10
-      ELSE IF( ( N1.GT.0 .AND. LDC1.LT.MAX( 1, P ) ) .OR.
-     $         ( N1.EQ.0 .AND. LDC1.LT.1 ) ) THEN
-         INFO = -12
-      ELSE IF( LDD1.LT.MAX( 1, P ) ) THEN
-         INFO = -14
-      ELSE IF( LDA2.LT.MAX( 1, N2 ) ) THEN
-         INFO = -16
-      ELSE IF( LDB2.LT.MAX( 1, N2 ) ) THEN
-         INFO = -18
-      ELSE IF( ( N2.GT.0 .AND. LDC2.LT.MAX( 1, P ) ) .OR.
-     $         ( N2.EQ.0 .AND. LDC2.LT.1 ) ) THEN
-         INFO = -20
-      ELSE IF( LDD2.LT.MAX( 1, P ) ) THEN
-         INFO = -22
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -25
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -27
-      ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, P ) ) .OR.
-     $         ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN
-         INFO = -29
-      ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
-         INFO = -31
-      END IF
-C
-      IF ( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB05PD', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF ( MAX( N, MIN( M, P ) ).EQ.0 )
-     $   RETURN
-C
-      N1P1 = N1 + 1
-C
-C                       ( A1   0  )
-C     Construct     A = (         ) .
-C                       ( 0    A2 )
-C
-      IF ( LOVER .AND. LDA1.LE.LDA ) THEN
-         IF ( LDA1.LT.LDA ) THEN
-C
-            DO 20 J = N1, 1, -1
-               DO 10 I = N1, 1, -1
-                  A(I,J) = A1(I,J)
-   10          CONTINUE
-   20       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', N1, N1, A1, LDA1, A, LDA )
-      END IF
-C
-      IF ( N2.GT.0 ) THEN
-         CALL DLASET( 'F', N1, N2, ZERO, ZERO, A(1,N1P1), LDA )
-         CALL DLASET( 'F', N2, N1, ZERO, ZERO, A(N1P1,1), LDA )
-         CALL DLACPY( 'F', N2, N2, A2, LDA2, A(N1P1,N1P1), LDA )
-      END IF
-C
-C                        ( B1 )
-C     Construct      B = (    ) .
-C                        ( B2 )
-C
-      IF ( LOVER .AND. LDB1.LE.LDB ) THEN
-         IF ( LDB1.LT.LDB ) THEN
-C
-            DO 40 J = M, 1, -1
-               DO 30 I = N1, 1, -1
-                  B(I,J) = B1(I,J)
-   30          CONTINUE
-   40       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', N1, M, B1, LDB1, B, LDB )
-      END IF
-C
-      IF ( N2.GT.0 )
-     $   CALL DLACPY( 'F', N2, M, B2, LDB2, B(N1P1,1), LDB )
-C
-C     Construct      C = ( C1 alpha*C2 ) .
-C
-      IF ( LOVER .AND. LDC1.LE.LDC ) THEN
-         IF ( LDC1.LT.LDC ) THEN
-C
-            DO 60 J = N1, 1, -1
-               DO 50 I = P, 1, -1
-                  C(I,J) = C1(I,J)
-   50          CONTINUE
-   60       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', P, N1, C1, LDC1, C, LDC )
-      END IF
-C
-      IF ( N2.GT.0 ) THEN
-         CALL DLACPY( 'F', P, N2, C2, LDC2, C(1,N1P1), LDC )
-         IF ( ALPHA.NE.ONE )
-     $      CALL DLASCL( 'G', 0, 0, ONE, ALPHA, P, N2, C(1,N1P1), LDC,
-     $                   INFO )
-      END IF
-C
-C     Construct       D = D1 + alpha*D2 .
-C
-      IF ( LOVER .AND. LDD1.LE.LDD ) THEN
-         IF ( LDD1.LT.LDD ) THEN
-C
-            DO 80 J = M, 1, -1
-               DO 70 I = P, 1, -1
-                  D(I,J) = D1(I,J)
-   70          CONTINUE
-   80       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', P, M, D1, LDD1, D, LDD )
-      END IF
-C
-      DO 90 J = 1, M
-         CALL DAXPY( P, ALPHA, D2(1,J), 1, D(1,J), 1 )
-   90 CONTINUE
-C
-      RETURN
-C *** Last line of AB05PD ***
-      END

--- a/extra/control-devel/src/AB05QD.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,419 +0,0 @@
-      SUBROUTINE AB05QD( OVER, N1, M1, P1, N2, M2, P2, A1, LDA1, B1,
-     $                   LDB1, C1, LDC1, D1, LDD1, A2, LDA2, B2, LDB2,
-     $                   C2, LDC2, D2, LDD2, N, M, P, A, LDA, B, LDB,
-     $                   C, LDC, D, LDD, 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 append two systems G1 and G2 in state-space form together.
-C     If G1 = (A1,B1,C1,D1) and G2 = (A2,B2,C2,D2) are the state-space
-C     models of the given two systems having the transfer-function
-C     matrices G1 and G2, respectively, this subroutine constructs the
-C     state-space model G = (A,B,C,D) which corresponds to the
-C     transfer-function matrix
-C
-C                           ( G1 0  )
-C                       G = (       )
-C                           ( 0  G2 )
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     OVER    CHARACTER*1
-C             Indicates whether the user wishes to overlap pairs of
-C             arrays, as follows:
-C             = 'N':  Do not overlap;
-C             = 'O':  Overlap pairs of arrays: A1 and A, B1 and B,
-C                     C1 and C, and D1 and D, i.e. the same name is
-C                     effectively used for each pair (for all pairs)
-C                     in the routine call.  In this case, setting
-C                     LDA1 = LDA, LDB1 = LDB, LDC1 = LDC, and LDD1 = LDD
-C                     will give maximum efficiency.
-C
-C     Input/Output Parameters
-C
-C     N1      (input) INTEGER
-C             The number of state variables in the first system, i.e.
-C             the order of the matrix A1, the number of rows of B1 and
-C             the number of columns of C1.  N1 >= 0.
-C
-C     M1      (input) INTEGER
-C             The number of input variables in the first system, i.e.
-C             the number of columns of matrices B1 and D1.  M1 >= 0.
-C
-C     P1      (input) INTEGER
-C             The number of output variables in the first system, i.e.
-C             the number of rows of matrices C1 and D1.  P1 >= 0.
-C
-C     N2      (input) INTEGER
-C             The number of state variables in the second system, i.e.
-C             the order of the matrix A2, the number of rows of B2 and
-C             the number of columns of C2.  N2 >= 0.
-C
-C     M2      (input) INTEGER
-C             The number of input variables in the second system, i.e.
-C             the number of columns of matrices B2 and D2.  M2 >= 0.
-C
-C     P2      (input) INTEGER
-C             The number of output variables in the second system, i.e.
-C             the number of rows of matrices C2 and D2.  P2 >= 0.
-C
-C     A1      (input) DOUBLE PRECISION array, dimension (LDA1,N1)
-C             The leading N1-by-N1 part of this array must contain the
-C             state transition matrix A1 for the first system.
-C
-C     LDA1    INTEGER
-C             The leading dimension of array A1.  LDA1 >= MAX(1,N1).
-C
-C     B1      (input) DOUBLE PRECISION array, dimension (LDB1,M1)
-C             The leading N1-by-M1 part of this array must contain the
-C             input/state matrix B1 for the first system.
-C
-C     LDB1    INTEGER
-C             The leading dimension of array B1.  LDB1 >= MAX(1,N1).
-C
-C     C1      (input) DOUBLE PRECISION array, dimension (LDC1,N1)
-C             The leading P1-by-N1 part of this array must contain the
-C             state/output matrix C1 for the first system.
-C
-C     LDC1    INTEGER
-C             The leading dimension of array C1.
-C             LDC1 >= MAX(1,P1) if N1 > 0.
-C             LDC1 >= 1 if N1 = 0.
-C
-C     D1      (input) DOUBLE PRECISION array, dimension (LDD1,M1)
-C             The leading P1-by-M1 part of this array must contain the
-C             input/output matrix D1 for the first system.
-C
-C     LDD1    INTEGER
-C             The leading dimension of array D1.  LDD1 >= MAX(1,P1).
-C
-C     A2      (input) DOUBLE PRECISION array, dimension (LDA2,N2)
-C             The leading N2-by-N2 part of this array must contain the
-C             state transition matrix A2 for the second system.
-C
-C     LDA2    INTEGER
-C             The leading dimension of array A2.  LDA2 >= MAX(1,N2).
-C
-C     B2      (input) DOUBLE PRECISION array, dimension (LDB2,M2)
-C             The leading N2-by-M2 part of this array must contain the
-C             input/state matrix B2 for the second system.
-C
-C     LDB2    INTEGER
-C             The leading dimension of array B2.  LDB2 >= MAX(1,N2).
-C
-C     C2      (input) DOUBLE PRECISION array, dimension (LDC2,N2)
-C             The leading P2-by-N2 part of this array must contain the
-C             state/output matrix C2 for the second system.
-C
-C     LDC2    INTEGER
-C             The leading dimension of array C2.
-C             LDC2 >= MAX(1,P2) if N2 > 0.
-C             LDC2 >= 1 if N2 = 0.
-C
-C     D2      (input) DOUBLE PRECISION array, dimension (LDD2,M2)
-C             The leading P2-by-M2 part of this array must contain the
-C             input/output matrix D2 for the second system.
-C
-C     LDD2    INTEGER
-C             The leading dimension of array D2.  LDD2 >= MAX(1,P2).
-C
-C     N       (output) INTEGER
-C             The number of state variables (N1 + N2) in the resulting
-C             system, i.e. the order of the matrix A, the number of rows
-C             of B and the number of columns of C.
-C
-C     M       (output) INTEGER
-C             The number of input variables (M1 + M2) in the resulting
-C             system, i.e. the number of columns of B and D.
-C
-C     P       (output) INTEGER
-C             The number of output variables (P1 + P2) of the resulting
-C             system, i.e. the number of rows of C and D.
-C
-C     A       (output) DOUBLE PRECISION array, dimension (LDA,N1+N2)
-C             The leading N-by-N part of this array contains the state
-C             transition matrix A for the resulting system.
-C             The array A can overlap A1 if OVER = 'O'.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N1+N2).
-C
-C     B       (output) DOUBLE PRECISION array, dimension (LDB,M1+M2)
-C             The leading N-by-M part of this array contains the
-C             input/state matrix B for the resulting system.
-C             The array B can overlap B1 if OVER = 'O'.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N1+N2).
-C
-C     C       (output) DOUBLE PRECISION array, dimension (LDC,N1+N2)
-C             The leading P-by-N part of this array contains the
-C             state/output matrix C for the resulting system.
-C             The array C can overlap C1 if OVER = 'O'.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.
-C             LDC >= MAX(1,P1+P2) if N1+N2 > 0.
-C             LDC >= 1 if N1+N2 = 0.
-C
-C     D       (output) DOUBLE PRECISION array, dimension (LDD,M1+M2)
-C             The leading P-by-M part of this array contains the
-C             input/output matrix D for the resulting system.
-C             The array D can overlap D1 if OVER = 'O'.
-C
-C     LDD     INTEGER
-C             The leading dimension of array D.  LDD >= MAX(1,P1+P2).
-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
-C     METHOD
-C
-C     The matrices of the resulting systems are determined as:
-C
-C           ( A1   0  )         ( B1  0  )
-C       A = (         ) ,   B = (        ) ,
-C           ( 0    A2 )         ( 0   B2 )
-C
-C           ( C1   0  )         ( D1  0  )
-C       C = (         ) ,   D = (        ) .
-C           ( 0    C2 )         ( 0   D2 )
-C
-C     REFERENCES
-C
-C     None
-C
-C     CONTRIBUTORS
-C
-C     A. Varga, German Aerospace Research Establishment,
-C     Oberpfaffenhofen, Germany, and V. Sima, Katholieke Univ. Leuven,
-C     Belgium, Nov. 1996.
-C
-C     REVISIONS
-C
-C     V. Sima, Research Institute for Informatics, Bucharest, Feb. 2004.
-C
-C     KEYWORDS
-C
-C     Multivariable system, state-space model, state-space
-C     representation.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ZERO
-      PARAMETER         ( ZERO=0.0D0 )
-C     .. Scalar Arguments ..
-      CHARACTER         OVER
-      INTEGER           INFO, LDA, LDA1, LDA2, LDB, LDB1, LDB2, LDC,
-     $                  LDC1, LDC2, LDD, LDD1, LDD2, M, M1, M2, N, N1,
-     $                  N2, P, P1, P2
-C     .. Array Arguments ..
-      DOUBLE PRECISION  A(LDA,*), A1(LDA1,*), A2(LDA2,*), B(LDB,*),
-     $                  B1(LDB1,*), B2(LDB2,*), C(LDC,*), C1(LDC1,*),
-     $                  C2(LDC2,*), D(LDD,*), D1(LDD1,*), D2(LDD2,*)
-C     .. Local Scalars ..
-      LOGICAL           LOVER
-      INTEGER           I, J
-C     .. External Functions ..
-      LOGICAL           LSAME
-      EXTERNAL          LSAME
-C     .. External Subroutines ..
-      EXTERNAL          DLACPY, DLASET, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         MAX, MIN
-C     .. Executable Statements ..
-C
-      LOVER = LSAME( OVER, 'O' )
-      N = N1 + N2
-      M = M1 + M2
-      P = P1 + P2
-      INFO = 0
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT.LOVER .AND. .NOT.LSAME( OVER, 'N' ) ) THEN
-         INFO = -1
-      ELSE IF( N1.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( M1.LT.0 ) THEN
-         INFO = -3
-      ELSE IF( P1.LT.0 ) THEN
-         INFO = -4
-      ELSE IF( N2.LT.0 ) THEN
-         INFO = -5
-      ELSE IF( M2.LT.0 ) THEN
-         INFO = -6
-      ELSE IF( P2.LT.0 ) THEN
-         INFO = -7
-      ELSE IF( LDA1.LT.MAX( 1, N1 ) ) THEN
-         INFO = -9
-      ELSE IF( LDB1.LT.MAX( 1, N1 ) ) THEN
-         INFO = -11
-      ELSE IF( ( N1.GT.0 .AND. LDC1.LT.MAX( 1, P1 ) ) .OR.
-     $         ( N1.EQ.0 .AND. LDC1.LT.1 ) ) THEN
-         INFO = -13
-      ELSE IF( LDD1.LT.MAX( 1, P1 ) ) THEN
-         INFO = -15
-      ELSE IF( LDA2.LT.MAX( 1, N2 ) ) THEN
-         INFO = -17
-      ELSE IF( LDB2.LT.MAX( 1, N2 ) ) THEN
-         INFO = -19
-      ELSE IF( ( N2.GT.0 .AND. LDC2.LT.MAX( 1, P2 ) ) .OR.
-     $         ( N2.EQ.0 .AND. LDC2.LT.1 ) ) THEN
-         INFO = -21
-      ELSE IF( LDD2.LT.MAX( 1, P2 ) ) THEN
-         INFO = -23
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -28
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -30
-      ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, P ) ) .OR.
-     $         ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN
-         INFO = -32
-      ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
-         INFO = -34
-      END IF
-C
-      IF ( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB05QD', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF ( MAX( N, MIN( M, P ) ).EQ.0 )
-     $   RETURN
-C                       ( A1   0  )
-C     Construct     A = (         ) .
-C                       ( 0    A2 )
-C
-      IF ( LOVER .AND. LDA1.LE.LDA ) THEN
-         IF ( LDA1.LT.LDA ) THEN
-C
-            DO 20 J = N1, 1, -1
-               DO 10 I = N1, 1, -1
-                  A(I,J) = A1(I,J)
-   10          CONTINUE
-   20       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', N1, N1, A1, LDA1, A, LDA )
-      END IF
-C
-      IF ( N2.GT.0 ) THEN
-         CALL DLASET( 'F', N1, N2, ZERO, ZERO, A(1,N1+1), LDA )
-         CALL DLASET( 'F', N2, N1, ZERO, ZERO, A(N1+1,1), LDA )
-         CALL DLACPY( 'F', N2, N2, A2, LDA2, A(N1+1,N1+1), LDA )
-      END IF
-C
-C                        ( B1  0  )
-C     Construct      B = (        ) .
-C                        ( 0   B2 )
-C
-      IF ( LOVER .AND. LDB1.LE.LDB ) THEN
-         IF ( LDB1.LT.LDB ) THEN
-C
-            DO 40 J = M1, 1, -1
-               DO 30 I = N1, 1, -1
-                  B(I,J) = B1(I,J)
-   30          CONTINUE
-   40       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', N1, M1, B1, LDB1, B, LDB )
-      END IF
-C
-      IF ( M2.GT.0 )
-     $   CALL DLASET( 'F', N1, M2, ZERO, ZERO, B(1,M1+1), LDB )
-      IF ( N2.GT.0 ) THEN
-         CALL DLASET( 'F', N2, M1, ZERO, ZERO, B(N1+1,1), LDB )
-         IF ( M2.GT.0 )
-     $      CALL DLACPY( 'F', N2, M2, B2, LDB2, B(N1+1,M1+1), LDB )
-      END IF
-C
-C                         ( C1   0  )
-C     Construct      C =  (         ) .
-C                         ( 0    C2 )
-C
-      IF ( LOVER .AND. LDC1.LE.LDC ) THEN
-         IF ( LDC1.LT.LDC ) THEN
-C
-            DO 60 J = N1, 1, -1
-               DO 50 I = P1, 1, -1
-                  C(I,J) = C1(I,J)
-   50          CONTINUE
-   60       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', P1, N1, C1, LDC1, C, LDC )
-      END IF
-C
-      IF ( N2.GT.0 )
-     $   CALL DLASET( 'F', P1, N2, ZERO, ZERO, C(1,N1+1), LDC )
-      IF ( P2.GT.0 ) THEN
-         IF ( N1.GT.0 )
-     $      CALL DLASET( 'F', P2, N1, ZERO, ZERO, C(P1+1,1), LDC )
-         IF ( N2.GT.0 )
-     $      CALL DLACPY( 'F', P2, N2, C2, LDC2, C(P1+1,N1+1), LDC )
-      END IF
-C
-C                          ( D1  0  )
-C     Construct       D =  (        ) .
-C                          ( 0   D2 )
-C
-      IF ( LOVER .AND. LDD1.LE.LDD ) THEN
-         IF ( LDD1.LT.LDD ) THEN
-C
-            DO 80 J = M1, 1, -1
-               DO 70 I = P1, 1, -1
-                  D(I,J) = D1(I,J)
-   70          CONTINUE
-   80       CONTINUE
-C
-         END IF
-      ELSE
-         CALL DLACPY( 'F', P1, M1, D1, LDD1, D, LDD )
-      END IF
-C
-      IF ( M2.GT.0 )
-     $   CALL DLASET( 'F', P1, M2, ZERO, ZERO, D(1,M1+1), LDD )
-      IF ( P2.GT.0 ) THEN
-         CALL DLASET( 'F', P2, M1, ZERO, ZERO, D(P1+1,1), LDD )
-         IF ( M2.GT.0 )
-     $      CALL DLACPY( 'F', P2, M2, D2, LDD2, D(P1+1,M1+1), LDD )
-      END IF
-C
-      RETURN
-C *** Last line of AB05QD ***
-      END

--- a/extra/control-devel/src/AB07MD.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,224 +0,0 @@
-      SUBROUTINE AB07MD( JOBD, N, M, P, A, LDA, B, LDB, C, LDC, D, LDD,
-     $                   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 find the dual of a given state-space representation.
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     JOBD    CHARACTER*1
-C             Specifies whether or not a non-zero matrix D appears in
-C             the given state space model:
-C             = 'D':  D is present;
-C             = 'Z':  D is assumed a zero matrix.
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The order of the state-space representation.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of system outputs.  P >= 0.
-C
-C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C             On entry, the leading N-by-N part of this array must
-C             contain the original state dynamics matrix A.
-C             On exit, the leading N-by-N part of this array contains
-C             the dual state dynamics matrix A'.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N).
-C
-C     B       (input/output) DOUBLE PRECISION array, dimension
-C             (LDB,MAX(M,P))
-C             On entry, the leading N-by-M part of this array must
-C             contain the original input/state matrix B.
-C             On exit, the leading N-by-P part of this array contains
-C             the dual input/state matrix C'.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N).
-C
-C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C             On entry, the leading P-by-N part of this array must
-C             contain the original state/output matrix C.
-C             On exit, the leading M-by-N part of this array contains
-C             the dual state/output matrix B'.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.
-C             LDC >= MAX(1,M,P) if N > 0.
-C             LDC >= 1 if N = 0.
-C
-C     D       (input/output) DOUBLE PRECISION array, dimension
-C             (LDD,MAX(M,P))
-C             On entry, if JOBD = 'D', the leading P-by-M part of this
-C             array must contain the original direct transmission
-C             matrix D.
-C             On exit, if JOBD = 'D', the leading M-by-P part of this
-C             array contains the dual direct transmission matrix D'.
-C             The array D is not referenced if JOBD = 'Z'.
-C
-C     LDD     INTEGER
-C             The leading dimension of array D.
-C             LDD >= MAX(1,M,P) if JOBD = 'D'.
-C             LDD >= 1 if JOBD = 'Z'.
-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
-C     METHOD
-C
-C     If the given state-space representation is the M-input/P-output
-C     (A,B,C,D), its dual is simply the P-input/M-output (A',C',B',D').
-C
-C     REFERENCES
-C
-C     None
-C
-C     NUMERICAL ASPECTS
-C
-C     None
-C
-C     CONTRIBUTOR
-C
-C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Dec. 1996.
-C     Supersedes Release 2.0 routine AB07AD by T.W.C.Williams, Kingston
-C     Polytechnic, United Kingdom, March 1982.
-C
-C     REVISIONS
-C
-C     V. Sima, Research Institute for Informatics, Bucharest, Feb. 2004.
-C
-C     KEYWORDS
-C
-C     Dual system, state-space model, state-space representation.
-C
-C     ******************************************************************
-C
-C     .. Scalar Arguments ..
-      CHARACTER         JOBD
-      INTEGER           INFO, LDA, LDB, LDC, LDD, M, N, P
-C     .. Array Arguments ..
-      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*)
-C     .. Local Scalars ..
-      LOGICAL           LJOBD
-      INTEGER           J, MINMP, MPLIM
-C     .. External functions ..
-      LOGICAL           LSAME
-      EXTERNAL          LSAME
-C     .. External subroutines ..
-      EXTERNAL          DCOPY, DSWAP, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         MAX, MIN
-C     .. Executable Statements ..
-C
-      INFO = 0
-      LJOBD = LSAME( JOBD, 'D' )
-      MPLIM = MAX( M, P )
-      MINMP = MIN( M, P )
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT.LJOBD .AND. .NOT.LSAME( JOBD, 'Z' )  ) THEN
-         INFO = -1
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -3
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -4
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -6
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -8
-      ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, MPLIM ) ) .OR.
-     $         ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN
-         INFO = -10
-      ELSE IF( ( LJOBD .AND. LDD.LT.MAX( 1, MPLIM ) ) .OR.
-     $    ( .NOT.LJOBD .AND. LDD.LT.1 ) ) THEN
-         INFO = -12
-      END IF
-C
-      IF ( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB07MD', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF ( MAX( N, MINMP ).EQ.0 )
-     $   RETURN
-C
-      IF ( N.GT.0 ) THEN
-C
-C        Transpose A, if non-scalar.
-C
-         DO 10 J = 1, N - 1
-            CALL DSWAP( N-J, A(J+1,J), 1, A(J,J+1), LDA )
-   10    CONTINUE
-C
-C        Replace B by C' and C by B'.
-C
-         DO 20 J = 1, MPLIM
-            IF ( J.LE.MINMP ) THEN
-               CALL DSWAP( N, B(1,J), 1, C(J,1), LDC )
-            ELSE IF ( J.GT.P ) THEN
-               CALL DCOPY( N, B(1,J), 1, C(J,1), LDC )
-            ELSE
-               CALL DCOPY( N, C(J,1), LDC, B(1,J), 1 )
-            END IF
-   20    CONTINUE
-C
-      END IF
-C
-      IF ( LJOBD .AND. MINMP.GT.0 ) THEN
-C
-C        Transpose D, if non-scalar.
-C
-         DO 30 J = 1, MPLIM
-            IF ( J.LT.MINMP ) THEN
-               CALL DSWAP( MINMP-J, D(J+1,J), 1, D(J,J+1), LDD )
-            ELSE IF ( J.GT.P ) THEN
-               CALL DCOPY( P, D(1,J), 1, D(J,1), LDD )
-            ELSE IF ( J.GT.M ) THEN
-               CALL DCOPY( M, D(J,1), LDD, D(1,J), 1 )
-            END IF
-   30    CONTINUE
-C
-      END IF
-C
-      RETURN
-C *** Last line of AB07MD ***
-      END

--- a/extra/control-devel/src/AB07ND.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,303 +0,0 @@
-      SUBROUTINE AB07ND( N, M, A, LDA, B, LDB, C, LDC, D, LDD, RCOND,
-     $                   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 compute the inverse (Ai,Bi,Ci,Di) of a given system (A,B,C,D).
-C
-C     ARGUMENTS
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The order of the state matrix A.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs and outputs.  M >= 0.
-C
-C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C             On entry, the leading N-by-N part of this array must
-C             contain the state matrix A of the original system.
-C             On exit, the leading N-by-N part of this array contains
-C             the state matrix Ai of the inverse system.
-C
-C     LDA     INTEGER
-C             The leading dimension of the array A.  LDA >= MAX(1,N).
-C
-C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C             On entry, the leading N-by-M part of this array must
-C             contain the input matrix B of the original system.
-C             On exit, the leading N-by-M part of this array contains
-C             the input matrix Bi of the inverse system.
-C
-C     LDB     INTEGER
-C             The leading dimension of the array B.  LDB >= MAX(1,N).
-C
-C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C             On entry, the leading M-by-N part of this array must
-C             contain the output matrix C of the original system.
-C             On exit, the leading M-by-N part of this array contains
-C             the output matrix Ci of the inverse system.
-C
-C     LDC     INTEGER
-C             The leading dimension of the array C.  LDC >= MAX(1,M).
-C
-C     D       (input/output) DOUBLE PRECISION array, dimension (LDD,M)
-C             On entry, the leading M-by-M part of this array must
-C             contain the feedthrough matrix D of the original system.
-C             On exit, the leading M-by-M part of this array contains
-C             the feedthrough matrix Di of the inverse system.
-C
-C     LDD     INTEGER
-C             The leading dimension of the array D.  LDD >= MAX(1,M).
-C
-C     RCOND   (output) DOUBLE PRECISION
-C             The estimated reciprocal condition number of the
-C             feedthrough matrix D of the original system.
-C
-C     Workspace
-C
-C     IWORK   INTEGER array, dimension (2*M)
-C
-C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
-C             On exit, if INFO = 0 or M+1, DWORK(1) returns the optimal
-C             value of LDWORK.
-C
-C     LDWORK  INTEGER
-C             The length of the array DWORK.  LDWORK >= MAX(1,4*M).
-C             For good performance, LDWORK should be larger.
-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             = i:  the matrix D is exactly singular; the (i,i) diagonal
-C                   element is zero, i <= M; RCOND was set to zero;
-C             = M+1:  the matrix D is numerically singular, i.e., RCOND
-C                   is less than the relative machine precision, EPS
-C                   (see LAPACK Library routine DLAMCH). The
-C                   calculations have been completed, but the results
-C                   could be very inaccurate.
-C
-C     METHOD
-C
-C     The matrices of the inverse system are computed with the formulas:
-C                   -1              -1         -1           -1
-C       Ai = A - B*D  *C,  Bi = -B*D  ,  Ci = D  *C,  Di = D  .
-C
-C     NUMERICAL ASPECTS
-C
-C     The accuracy depends mainly on the condition number of the matrix
-C     D to be inverted. The estimated reciprocal condition number is
-C     returned in RCOND.
-C
-C     CONTRIBUTORS
-C
-C     A. Varga, German Aerospace Center, Oberpfaffenhofen, March 2000.
-C     D. Sima, University of Bucharest, April 2000.
-C     V. Sima, Research Institute for Informatics, Bucharest, Apr. 2000.
-C     Based on the routine SYSINV, A. Varga, 1992.
-C
-C     REVISIONS
-C
-C     A. Varga, German Aerospace Center, Oberpfaffenhofen, July 2000.
-C
-C     KEYWORDS
-C
-C     Inverse system, state-space model, state-space representation.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION   ZERO, ONE
-      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
-C     .. Scalar Arguments ..
-      DOUBLE PRECISION   RCOND
-      INTEGER            INFO, LDA, LDB, LDC, LDD, LDWORK, M, N
-C     .. Array Arguments ..
-      DOUBLE PRECISION   A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
-     $                   DWORK(*)
-      INTEGER            IWORK(*)
-C     .. Local Scalars ..
-      DOUBLE PRECISION   DNORM
-      INTEGER            BL, CHUNK, I, IERR, J, MAXWRK
-      LOGICAL            BLAS3, BLOCK
-C     .. External Functions ..
-      DOUBLE PRECISION   DLAMCH, DLANGE
-      INTEGER            ILAENV
-      EXTERNAL           DLAMCH, DLANGE, ILAENV
-C     .. External Subroutines ..
-      EXTERNAL           DCOPY, DGECON, DGEMM, DGEMV, DGETRF, DGETRI,
-     $                   DLACPY, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC          DBLE, MAX, MIN
-C     .. Executable Statements ..
-C
-      INFO = 0
-C
-C     Test the input scalar arguments.
-C
-      IF( N.LT.0 ) THEN
-         INFO = -1
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -4
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -6
-      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
-         INFO = -8
-      ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
-         INFO = -10
-      ELSE IF( LDWORK.LT.MAX( 1, 4*M ) ) THEN
-         INFO = -14
-      END IF
-C
-      IF ( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB07ND', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF ( M.EQ.0 ) THEN
-         RCOND    = ONE
-         DWORK(1) = ONE
-         RETURN
-      END IF
-C
-C     Factorize D.
-C
-      CALL DGETRF( M, M, D, LDD, IWORK, INFO )
-      IF ( INFO.NE.0 ) THEN
-         RCOND = ZERO
-         RETURN
-      END IF
-C
-C     Compute the reciprocal condition number of the matrix D.
-C     Workspace: need   4*M.
-C     (Note: Comments in the code beginning "Workspace:" describe the
-C      minimal amount of workspace needed at that point in the code,
-C      as well as the preferred amount for good performance.
-C      NB refers to the optimal block size for the immediately
-C      following subroutine, as returned by ILAENV.)
-C
-      DNORM = DLANGE( '1-norm', M, M, D, LDD, DWORK )
-      CALL DGECON( '1-norm', M, D, LDD, DNORM, RCOND, DWORK, IWORK(M+1),
-     $             IERR )
-      IF ( RCOND.LT.DLAMCH( 'Epsilon' ) )
-     $   INFO = M + 1
-C                   -1
-C     Compute Di = D  .
-C     Workspace: need   M;
-C                prefer M*NB.
-C
-      MAXWRK = MAX( 4*M, M*ILAENV( 1, 'DGETRI', ' ', M, -1, -1, -1 ) )
-      CALL DGETRI( M, D, LDD, IWORK, DWORK, LDWORK, IERR )
-      IF ( N.GT.0 ) THEN
-         CHUNK = LDWORK / M
-         BLAS3 = CHUNK.GE.N .AND. M.GT.1
-         BLOCK = MIN( CHUNK, M ).GT.1
-C                          -1
-C        Compute  Bi = -B*D  .
-C
-         IF ( BLAS3 ) THEN
-C
-C           Enough workspace for a fast BLAS 3 algorithm.
-C
-            CALL DLACPY( 'Full', N, M, B, LDB, DWORK, N )
-            CALL DGEMM( 'NoTranspose', 'NoTranspose', N, M, M, -ONE,
-     $                  DWORK, N, D, LDD, ZERO, B, LDB )
-C
-         ELSE IF( BLOCK ) THEN
-C
-C           Use as many rows of B as possible.
-C
-            DO 10 I = 1, N, CHUNK
-               BL = MIN( N-I+1, CHUNK )
-               CALL DLACPY( 'Full', BL, M, B(I,1), LDB, DWORK, BL )
-               CALL DGEMM( 'NoTranspose', 'NoTranspose', BL, M, M, -ONE,
-     $                     DWORK, BL, D, LDD, ZERO, B(I,1), LDB )
-   10       CONTINUE
-C
-         ELSE
-C
-C           Use a BLAS 2 algorithm.
-C
-            DO 20 I = 1, N
-               CALL DCOPY( M, B(I,1), LDB, DWORK, 1 )
-               CALL DGEMV( 'Transpose', M, M, -ONE, D, LDD, DWORK, 1,
-     $                     ZERO, B(I,1), LDB )
-   20       CONTINUE
-C
-         END IF
-C
-C        Compute  Ai = A + Bi*C.
-C
-         CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, M, ONE, B, LDB,
-     $               C, LDC, ONE, A, LDA )
-C                        -1
-C        Compute  C <-- D  *C.
-C
-         IF ( BLAS3 ) THEN
-C
-C           Enough workspace for a fast BLAS 3 algorithm.
-C
-            CALL DLACPY( 'Full', M, N, C, LDC, DWORK, M )
-            CALL DGEMM( 'NoTranspose', 'NoTranspose', M, N, M, ONE,
-     $                  D, LDD, DWORK, M, ZERO, C, LDC )
-C
-         ELSE IF( BLOCK ) THEN
-C
-C           Use as many columns of C as possible.
-C
-            DO 30 J = 1, N, CHUNK
-               BL = MIN( N-J+1, CHUNK )
-               CALL DLACPY( 'Full', M, BL, C(1,J), LDC, DWORK, M )
-               CALL DGEMM( 'NoTranspose', 'NoTranspose', M, BL, M, ONE,
-     $                     D, LDD, DWORK, M, ZERO, C(1,J), LDC )
-   30       CONTINUE
-C
-         ELSE
-C
-C           Use a BLAS 2 algorithm.
-C
-            DO 40 J = 1, N
-               CALL DCOPY( M, C(1,J), 1, DWORK, 1 )
-               CALL DGEMV( 'NoTranspose', M, M, ONE, D, LDD, DWORK, 1,
-     $                     ZERO, C(1,J), 1 )
-   40       CONTINUE
-C
-         END IF
-      END IF
-C
-C     Return optimal workspace in DWORK(1).
-C
-      DWORK(1) = DBLE( MAX( MAXWRK, N*M ) )
-      RETURN
-C
-C *** Last line of AB07ND ***
-      END

--- a/extra/control-devel/src/AB08MD.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,299 +0,0 @@
-      SUBROUTINE AB08MD( EQUIL, N, M, P, A, LDA, B, LDB, C, LDC, D, LDD,
-     $                   RANK, TOL, 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 compute the normal rank of the transfer-function matrix of a
-C     state-space model (A,B,C,D).
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     EQUIL   CHARACTER*1
-C             Specifies whether the user wishes to balance the compound
-C             matrix (see METHOD) as follows:
-C             = 'S':  Perform balancing (scaling);
-C             = 'N':  Do not perform balancing.
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The number of state variables, i.e., the order of the
-C             matrix A.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of system outputs.  P >= 0.
-C
-C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
-C             The leading N-by-N part of this array must contain the
-C             state dynamics matrix A of the system.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N).
-C
-C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
-C             The leading N-by-M part of this array must contain the
-C             input/state matrix B of the system.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N).
-C
-C     C       (input) DOUBLE PRECISION array, dimension (LDC,N)
-C             The leading P-by-N part of this array must contain the
-C             state/output matrix C of the system.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.  LDC >= MAX(1,P).
-C
-C     D       (input) DOUBLE PRECISION array, dimension (LDD,M)
-C             The leading P-by-M part of this array must contain the
-C             direct transmission matrix D of the system.
-C
-C     LDD     INTEGER
-C             The leading dimension of array D.  LDD >= MAX(1,P).
-C
-C     RANK    (output) INTEGER
-C             The normal rank of the transfer-function matrix.
-C
-C     Tolerances
-C
-C     TOL     DOUBLE PRECISION
-C             A tolerance used in rank decisions to determine the
-C             effective rank, which is defined as the order of the
-C             largest leading (or trailing) triangular submatrix in the
-C             QR (or RQ) factorization with column (or row) pivoting
-C             whose estimated condition number is less than 1/TOL.
-C             If the user sets TOL to be less than SQRT((N+P)*(N+M))*EPS
-C             then the tolerance is taken as SQRT((N+P)*(N+M))*EPS,
-C             where EPS is the machine precision (see LAPACK Library
-C             Routine DLAMCH).
-C
-C     Workspace
-C
-C     IWORK   INTEGER array, dimension (2*N+MAX(M,P)+1)
-C
-C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
-C             On exit, if INFO = 0, DWORK(1) returns the optimal value
-C             of LDWORK.
-C
-C     LDWORK  INTEGER
-C             The length of the array DWORK.
-C             LDWORK >= (N+P)*(N+M) +
-C                       MAX( MIN(P,M) + MAX(3*M-1,N), 1,
-C                            MIN(P,N) + MAX(3*P-1,N+P,N+M) )
-C             For optimum performance LDWORK should be larger.
-C
-C             If LDWORK = -1, then a workspace query is assumed;
-C             the routine only calculates the optimal size of the
-C             DWORK array, returns this value as the first entry of
-C             the DWORK array, and no error message related to LDWORK
-C             is issued by XERBLA.
-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
-C     METHOD
-C
-C     The routine reduces the (N+P)-by-(M+N) compound matrix (B  A)
-C                                                            (D  C)
-C
-C     to one with the same invariant zeros and with D of full row rank.
-C     The normal rank of the transfer-function matrix is the rank of D.
-C
-C     REFERENCES
-C
-C     [1] Svaricek, F.
-C         Computation of the Structural Invariants of Linear
-C         Multivariable Systems with an Extended Version of
-C         the Program ZEROS.
-C         System & Control Letters, 6, pp. 261-266, 1985.
-C
-C     [2] Emami-Naeini, A. and Van Dooren, P.
-C         Computation of Zeros of Linear Multivariable Systems.
-C         Automatica, 18, pp. 415-430, 1982.
-C
-C     NUMERICAL ASPECTS
-C
-C     The algorithm is backward stable (see [2] and [1]).
-C
-C     CONTRIBUTOR
-C
-C     A. Varga, German Aerospace Center, Oberpfaffenhofen, May 2001.
-C
-C     REVISIONS
-C
-C     V. Sima, Research Institute for Informatics, Bucharest, June 2001,
-C     Dec. 2003, Jan. 2009, Mar. 2009, Apr. 2009.
-C
-C     KEYWORDS
-C
-C     Multivariable system, orthogonal transformation,
-C     structural invariant.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ZERO, ONE
-      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
-C     .. Scalar Arguments ..
-      CHARACTER         EQUIL
-      INTEGER           INFO, LDA, LDB, LDC, LDD, LDWORK, M, N, P, RANK
-      DOUBLE PRECISION  TOL
-C     .. Array Arguments ..
-      INTEGER           IWORK(*)
-      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), DWORK(*)
-C     .. Local Scalars ..
-      LOGICAL           LEQUIL, LQUERY
-      INTEGER           I, KW, MU, NB, NINFZ, NKROL, NM, NP, NU, RO,
-     $                  SIGMA, WRKOPT
-      DOUBLE PRECISION  MAXRED, SVLMAX, THRESH, TOLER
-C     .. External Functions ..
-      LOGICAL           LSAME
-      DOUBLE PRECISION  DLAMCH, DLANGE
-      EXTERNAL          DLAMCH, DLANGE, LSAME
-C     .. External Subroutines ..
-      EXTERNAL          AB08NX, DLACPY, TB01ID, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         DBLE, INT, MAX, MIN, SQRT
-C     .. Executable Statements ..
-C
-      NP = N + P
-      NM = N + M
-      INFO = 0
-      LEQUIL = LSAME( EQUIL, 'S' )
-      LQUERY = ( LDWORK.EQ.-1 )
-      WRKOPT = NP*NM
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT.LEQUIL .AND. .NOT.LSAME( EQUIL, 'N' ) ) THEN
-         INFO = -1
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -3
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -4
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -6
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -8
-      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
-         INFO = -10
-      ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
-         INFO = -12
-      ELSE
-         KW = WRKOPT + MAX( MIN( P, M ) + MAX( 3*M-1, N ), 1,
-     $                      MIN( P, N ) + MAX( 3*P-1, NP, NM ) )
-         IF( LQUERY ) THEN
-            SVLMAX = ZERO
-            NINFZ  = 0
-            CALL AB08NX( N, M, P, P, 0, SVLMAX, DWORK, MAX( 1, NP ),
-     $                   NINFZ, IWORK, IWORK, MU, NU, NKROL, TOL, IWORK,
-     $                   DWORK, -1, INFO )
-            WRKOPT = MAX( KW, WRKOPT + INT( DWORK(1) ) )
-         ELSE IF( LDWORK.LT.KW ) THEN
-            INFO = -17
-         END IF
-      END IF
-C
-      IF ( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB08MD', -INFO )
-         RETURN
-      ELSE IF( LQUERY ) THEN
-         DWORK(1) = WRKOPT
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF ( MIN( M, P ).EQ.0 ) THEN
-         RANK = 0
-         DWORK(1) = ONE
-         RETURN
-      END IF
-C
-      DO 10 I = 1, 2*N+1
-         IWORK(I) = 0
-   10 CONTINUE
-C
-C     (Note: Comments in the code beginning "Workspace:" describe the
-C     minimal amount of real workspace needed at that point in the
-C     code, as well as the preferred amount for good performance.)
-C
-C     Construct the compound matrix  ( B  A ), dimension (N+P)-by-(M+N).
-C                                    ( D  C )
-C     Workspace: need   (N+P)*(N+M).
-C
-      CALL DLACPY( 'Full', N, M, B, LDB, DWORK, NP )
-      CALL DLACPY( 'Full', P, M, D, LDD, DWORK(N+1), NP )
-      CALL DLACPY( 'Full', N, N, A, LDA, DWORK(NP*M+1), NP )
-      CALL DLACPY( 'Full', P, N, C, LDC, DWORK(NP*M+N+1), NP )
-C
-C     If required, balance the compound matrix (default MAXRED).
-C     Workspace: need   N.
-C
-      KW = WRKOPT + 1
-      IF ( LEQUIL ) THEN
-         MAXRED = ZERO
-         CALL TB01ID( 'A', N, M, P, MAXRED, DWORK(NP*M+1), NP, DWORK,
-     $                NP, DWORK(NP*M+N+1), NP, DWORK(KW), INFO )
-         WRKOPT = WRKOPT + N
-      END IF
-C
-C     If required, set tolerance.
-C
-      THRESH = SQRT( DBLE( NP*NM ) )*DLAMCH( 'Precision' )
-      TOLER = TOL
-      IF ( TOLER.LT.THRESH ) TOLER = THRESH
-      SVLMAX = DLANGE( 'Frobenius', NP, NM, DWORK, NP, DWORK(KW) )
-C
-C     Reduce this system to one with the same invariant zeros and with
-C     D full row rank MU (the normal rank of the original system).
-C     Real workspace: need   (N+P)*(N+M) +
-C                            MAX( 1, MIN(P,M) + MAX(3*M-1,N),
-C                                    MIN(P,N) + MAX(3*P-1,N+P,N+M) );
-C                     prefer larger.
-C     Integer workspace: 2*N+MAX(M,P)+1.
-C
-      RO = P
-      SIGMA = 0
-      NINFZ = 0
-      CALL AB08NX( N, M, P, RO, SIGMA, SVLMAX, DWORK, NP, NINFZ, IWORK,
-     $             IWORK(N+1), MU, NU, NKROL, TOLER, IWORK(2*N+2),
-     $             DWORK(KW), LDWORK-KW+1, INFO )
-      RANK = MU
-C
-      DWORK(1) = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
-      RETURN
-C *** Last line of AB08MD ***
-      END

--- a/extra/control-devel/src/AB08NX.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,446 +0,0 @@
-      SUBROUTINE AB08NX( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD,
-     $                   NINFZ, INFZ, KRONL, MU, NU, NKROL, TOL, 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 extract from the (N+P)-by-(M+N) system
-C                  ( B  A )
-C                  ( D  C )
-C     an (NU+MU)-by-(M+NU) "reduced" system
-C                  ( B' A')
-C                  ( D' C')
-C     having the same transmission zeros but with D' of full row rank.
-C
-C     ARGUMENTS
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The number of state variables.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of system outputs.  P >= 0.
-C
-C     RO      (input/output) INTEGER
-C             On entry,
-C             = P     for the original system;
-C             = MAX(P-M, 0) for the pertransposed system.
-C             On exit, RO contains the last computed rank.
-C
-C     SIGMA   (input/output) INTEGER
-C             On entry,
-C             = 0  for the original system;
-C             = M  for the pertransposed system.
-C             On exit, SIGMA contains the last computed value sigma in
-C             the algorithm.
-C
-C     SVLMAX  (input) DOUBLE PRECISION
-C             During each reduction step, the rank-revealing QR
-C             factorization of a matrix stops when the estimated minimum
-C             singular value is smaller than TOL * MAX(SVLMAX,EMSV),
-C             where EMSV is the estimated maximum singular value.
-C             SVLMAX >= 0.
-C
-C     ABCD    (input/output) DOUBLE PRECISION array, dimension
-C             (LDABCD,M+N)
-C             On entry, the leading (N+P)-by-(M+N) part of this array
-C             must contain the compound input matrix of the system.
-C             On exit, the leading (NU+MU)-by-(M+NU) part of this array
-C             contains the reduced compound input matrix of the system.
-C
-C     LDABCD  INTEGER
-C             The leading dimension of array ABCD.
-C             LDABCD >= MAX(1,N+P).
-C
-C     NINFZ   (input/output) INTEGER
-C             On entry, the currently computed number of infinite zeros.
-C             It should be initialized to zero on the first call.
-C             NINFZ >= 0.
-C             On exit, the number of infinite zeros.
-C
-C     INFZ    (input/output) INTEGER array, dimension (N)
-C             On entry, INFZ(i) must contain the current number of
-C             infinite zeros of degree i, where i = 1,2,...,N, found in
-C             the previous call(s) of the routine. It should be
-C             initialized to zero on the first call.
-C             On exit, INFZ(i) contains the number of infinite zeros of
-C             degree i, where i = 1,2,...,N.
-C
-C     KRONL   (input/output) INTEGER array, dimension (N+1)
-C             On entry, this array must contain the currently computed
-C             left Kronecker (row) indices found in the previous call(s)
-C             of the routine. It should be initialized to zero on the
-C             first call.
-C             On exit, the leading NKROL elements of this array contain
-C             the left Kronecker (row) indices.
-C
-C     MU      (output) INTEGER
-C             The normal rank of the transfer function matrix of the
-C             original system.
-C
-C     NU      (output) INTEGER
-C             The dimension of the reduced system matrix and the number
-C             of (finite) invariant zeros if D' is invertible.
-C
-C     NKROL   (output) INTEGER
-C             The number of left Kronecker indices.
-C
-C     Tolerances
-C
-C     TOL     DOUBLE PRECISION
-C             A tolerance used in rank decisions to determine the
-C             effective rank, which is defined as the order of the
-C             largest leading (or trailing) triangular submatrix in the
-C             QR (or RQ) factorization with column (or row) pivoting
-C             whose estimated condition number is less than 1/TOL.
-C             NOTE that when SVLMAX > 0, the estimated ranks could be
-C             less than those defined above (see SVLMAX).
-C
-C     Workspace
-C
-C     IWORK   INTEGER array, dimension (MAX(M,P))
-C
-C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
-C             On exit, if INFO = 0, DWORK(1) returns the optimal value
-C             of LDWORK.
-C
-C     LDWORK  INTEGER
-C             The length of the array DWORK.
-C             LDWORK >= MAX( 1, MIN(P,M) + MAX(3*M-1,N),
-C                               MIN(P,N) + MAX(3*P-1,N+P,N+M) ).
-C             For optimum performance LDWORK should be larger.
-C
-C             If LDWORK = -1, then a workspace query is assumed;
-C             the routine only calculates the optimal size of the
-C             DWORK array, returns this value as the first entry of
-C             the DWORK array, and no error message related to LDWORK
-C             is issued by XERBLA.
-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
-C     REFERENCES
-C
-C     [1] Svaricek, F.
-C         Computation of the Structural Invariants of Linear
-C         Multivariable Systems with an Extended Version of
-C         the Program ZEROS.
-C         System & Control Letters, 6, pp. 261-266, 1985.
-C
-C     [2] Emami-Naeini, A. and Van Dooren, P.
-C         Computation of Zeros of Linear Multivariable Systems.
-C         Automatica, 18, pp. 415-430, 1982.
-C
-C     NUMERICAL ASPECTS
-C
-C     The algorithm is backward stable.
-C
-C     CONTRIBUTOR
-C
-C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996.
-C     Supersedes Release 2.0 routine AB08BZ by F. Svaricek.
-C
-C     REVISIONS
-C
-C     V. Sima, Oct. 1997; Feb. 1998, Jan. 2009, Apr. 2009.
-C     A. Varga, May 1999; May 2001.
-C
-C     KEYWORDS
-C
-C     Generalized eigenvalue problem, Kronecker indices, multivariable
-C     system, orthogonal transformation, structural invariant.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ZERO
-      PARAMETER         ( ZERO = 0.0D0 )
-C     .. Scalar Arguments ..
-      INTEGER           INFO, LDABCD, LDWORK, M, MU, N, NINFZ, NKROL,
-     $                  NU, P, RO, SIGMA
-      DOUBLE PRECISION  SVLMAX, TOL
-C     .. Array Arguments ..
-      INTEGER           INFZ(*), IWORK(*), KRONL(*)
-      DOUBLE PRECISION  ABCD(LDABCD,*), DWORK(*)
-C     .. Local Scalars ..
-      LOGICAL           LQUERY
-      INTEGER           I1, IK, IROW, ITAU, IZ, JWORK, MM1, MNTAU, MNU,
-     $                  MPM, NB, NP, RANK, RO1, TAU, WRKOPT
-      DOUBLE PRECISION  T
-C     .. Local Arrays ..
-      DOUBLE PRECISION  SVAL(3)
-C     .. External Functions ..
-      INTEGER           ILAENV
-      EXTERNAL          ILAENV
-C     .. External Subroutines ..
-      EXTERNAL          DLAPMT, DLARFG, DLASET, DLATZM, DORMQR, DORMRQ,
-     $                  MB03OY, MB03PY, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         INT, MAX, MIN
-C     .. Executable Statements ..
-C
-      NP   = N + P
-      MPM  = MIN( P, M )
-      INFO = 0
-      LQUERY = ( LDWORK.EQ.-1 )
-C
-C     Test the input scalar arguments.
-C
-      IF( N.LT.0 ) THEN
-         INFO = -1
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -2
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -3
-      ELSE IF( RO.NE.P .AND. RO.NE.MAX( P-M, 0 ) ) THEN
-         INFO = -4
-      ELSE IF( SIGMA.NE.0 .AND. SIGMA.NE.M ) THEN
-         INFO = -5
-      ELSE IF( SVLMAX.LT.ZERO ) THEN
-         INFO = -6
-      ELSE IF( LDABCD.LT.MAX( 1, NP ) ) THEN
-         INFO = -8
-      ELSE IF( NINFZ.LT.0 ) THEN
-         INFO = -9
-      ELSE
-         JWORK = MAX( 1,      MPM + MAX( 3*M - 1, N ),
-     $                MIN( P, N ) + MAX( 3*P - 1, NP, N+M ) )
-         IF( LQUERY ) THEN
-            IF( M.GT.0 ) THEN
-               NB = MIN( 64, ILAENV( 1, 'DORMQR', 'LT', P, N, MPM,
-     $                               -1 ) )
-               WRKOPT = MAX( JWORK, MPM + MAX( 1, N )*NB )
-            ELSE
-               WRKOPT = JWORK
-            END IF
-            NB = MIN( 64, ILAENV( 1, 'DORMRQ', 'RT', NP, N, MIN( P, N ),
-     $                            -1 ) )
-            WRKOPT = MAX( WRKOPT, MIN( P, N ) + MAX( 1, NP )*NB )
-            NB = MIN( 64, ILAENV( 1, 'DORMRQ', 'LN', N, M+N,
-     $                            MIN( P, N ), -1 ) )
-            WRKOPT = MAX( WRKOPT, MIN( P, N ) + MAX( 1, M+N )*NB )
-         ELSE IF( LDWORK.LT.JWORK ) THEN
-            INFO = -18
-         END IF
-      END IF
-C
-      IF ( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB08NX', -INFO )
-         RETURN
-      ELSE IF( LQUERY ) THEN
-         DWORK(1) = WRKOPT
-         RETURN
-      END IF
-C
-      MU = P
-      NU = N
-C
-      IZ = 0
-      IK = 1
-      MM1 = M + 1
-      ITAU = 1
-      NKROL = 0
-      WRKOPT = 1
-C
-C     Main reduction loop:
-C
-C            M   NU                  M     NU
-C      NU  [ B   A ]           NU  [ B     A ]
-C      MU  [ D   C ]  -->    SIGMA [ RD   C1 ]   (SIGMA = rank(D) =
-C                             TAU  [ 0    C2 ]    row size of RD)
-C
-C                                    M   NU-RO  RO
-C                            NU-RO [ B1   A11  A12 ]
-C                     -->      RO  [ B2   A21  A22 ]  (RO = rank(C2) =
-C                            SIGMA [ RD   C11  C12 ]   col size of LC)
-C                             TAU  [ 0     0   LC  ]
-C
-C                                      M   NU-RO
-C                            NU-RO [ B1   A11 ]     NU := NU - RO
-C                                  [----------]     MU := RO + SIGMA
-C                     -->      RO  [ B2   A21 ]      D := [B2;RD]
-C                            SIGMA [ RD   C11 ]      C := [A21;C11]
-C
-   20 IF ( MU.EQ.0 )
-     $   GO TO 80
-C
-C     (Note: Comments in the code beginning "Workspace:" describe the
-C     minimal amount of real workspace needed at that point in the
-C     code, as well as the preferred amount for good performance.)
-C
-      RO1 = RO
-      MNU = M + NU
-      IF ( M.GT.0 ) THEN
-         IF ( SIGMA.NE.0 ) THEN
-            IROW = NU + 1
-C
-C           Compress rows of D.  First exploit triangular shape.
-C           Workspace: need   M+N-1.
-C
-            DO 40 I1 = 1, SIGMA
-               CALL DLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1, T )
-               CALL DLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, T,
-     $                      ABCD(IROW,I1+1), ABCD(IROW+1,I1+1), LDABCD,
-     $                      DWORK )
-               IROW = IROW + 1
-   40       CONTINUE
-            CALL DLASET( 'Lower', RO+SIGMA-1, SIGMA, ZERO, ZERO,
-     $                   ABCD(NU+2,1), LDABCD )
-         END IF
-C
-C        Continue with Householder with column pivoting.
-C
-C        The rank of D is the number of (estimated) singular values
-C        that are greater than TOL * MAX(SVLMAX,EMSV). This number
-C        includes the singular values of the first SIGMA columns.
-C        Integer workspace: need   M;
-C        Workspace: need   min(RO1,M) + 3*M - 1.  RO1 <= P.
-C
-         IF ( SIGMA.LT.M ) THEN
-            JWORK = ITAU + MIN( RO1, M )
-            I1    = SIGMA + 1
-            IROW  = NU + I1
-            CALL MB03OY( RO1, M-SIGMA, ABCD(IROW,I1), LDABCD, TOL,
-     $                   SVLMAX, RANK, SVAL, IWORK, DWORK(ITAU),
-     $                   DWORK(JWORK), INFO )
-            WRKOPT = MAX( WRKOPT, JWORK + 3*M - 2 )
-C
-C           Apply the column permutations to matrices B and part of D.
-C
-            CALL DLAPMT( .TRUE., NU+SIGMA, M-SIGMA, ABCD(1,I1), LDABCD,
-     $                   IWORK )
-C
-            IF ( RANK.GT.0 ) THEN
-C
-C              Apply the Householder transformations to the submatrix C.
-C              Workspace: need   min(RO1,M) + NU;
-C                         prefer min(RO1,M) + NU*NB.
-C
-               CALL DORMQR( 'Left', 'Transpose', RO1, NU, RANK,
-     $                      ABCD(IROW,I1), LDABCD, DWORK(ITAU),
-     $                      ABCD(IROW,MM1), LDABCD, DWORK(JWORK),
-     $                      LDWORK-JWORK+1, INFO )
-               WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
-               IF ( RO1.GT.1 )
-     $            CALL DLASET( 'Lower', RO1-1, MIN( RO1-1, RANK ), ZERO,
-     $                         ZERO, ABCD(IROW+1,I1), LDABCD )
-               RO1 = RO1 - RANK
-            END IF
-         END IF
-      END IF
-C
-      TAU = RO1
-      SIGMA = MU - TAU
-C
-C     Determination of the orders of the infinite zeros.
-C
-      IF ( IZ.GT.0 ) THEN
-         INFZ(IZ) = INFZ(IZ) + RO - TAU
-         NINFZ = NINFZ + IZ*( RO - TAU )
-      END IF
-      IF ( RO1.EQ.0 )
-     $   GO TO 80
-      IZ = IZ + 1
-C
-      IF ( NU.LE.0 ) THEN
-         MU = SIGMA
-         NU = 0
-         RO = 0
-      ELSE
-C
-C        Compress the columns of C2 using RQ factorization with row
-C        pivoting, P * C2 = R * Q.
-C
-         I1 = NU + SIGMA + 1
-         MNTAU = MIN( TAU, NU )
-         JWORK = ITAU + MNTAU
-C
-C        The rank of C2 is the number of (estimated) singular values
-C        greater than TOL * MAX(SVLMAX,EMSV).
-C        Integer Workspace: need TAU;
-C        Workspace: need min(TAU,NU) + 3*TAU - 1.
-C
-         CALL MB03PY( TAU, NU, ABCD(I1,MM1), LDABCD, TOL, SVLMAX, RANK,
-     $                SVAL, IWORK, DWORK(ITAU), DWORK(JWORK), INFO )
-         WRKOPT = MAX( WRKOPT, JWORK + 3*TAU - 1 )
-         IF ( RANK.GT.0 ) THEN
-            IROW = I1 + TAU - RANK
-C
-C           Apply Q' to the first NU columns of [A; C1] from the right.
-C           Workspace: need   min(TAU,NU) + NU + SIGMA; SIGMA <= P;
-C                      prefer min(TAU,NU) + (NU  + SIGMA)*NB.
-C
-            CALL DORMRQ( 'Right', 'Transpose', I1-1, NU, RANK,
-     $                   ABCD(IROW,MM1), LDABCD, DWORK(MNTAU-RANK+1),
-     $                   ABCD(1,MM1), LDABCD, DWORK(JWORK),
-     $                   LDWORK-JWORK+1, INFO )
-            WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
-C
-C           Apply Q to the first NU rows and M + NU columns of [ B  A ]
-C           from the left.
-C           Workspace: need   min(TAU,NU) + M + NU;
-C                      prefer min(TAU,NU) + (M + NU)*NB.
-C
-            CALL DORMRQ( 'Left', 'NoTranspose', NU, MNU, RANK,
-     $                   ABCD(IROW,MM1), LDABCD, DWORK(MNTAU-RANK+1),
-     $                   ABCD, LDABCD, DWORK(JWORK), LDWORK-JWORK+1,
-     $                   INFO )
-            WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
-C
-            CALL DLASET( 'Full', RANK, NU-RANK, ZERO, ZERO,
-     $                   ABCD(IROW,MM1), LDABCD )
-            IF ( RANK.GT.1 )
-     $         CALL DLASET( 'Lower', RANK-1, RANK-1, ZERO, ZERO,
-     $                      ABCD(IROW+1,MM1+NU-RANK), LDABCD )
-         END IF
-C
-         RO = RANK
-      END IF
-C
-C     Determine the left Kronecker indices (row indices).
-C
-      KRONL(IK) = KRONL(IK) + TAU - RO
-      NKROL = NKROL + KRONL(IK)
-      IK = IK + 1
-C
-C     C and D are updated to [A21 ; C11] and [B2 ; RD].
-C
-      NU = NU - RO
-      MU = SIGMA + RO
-      IF ( RO.NE.0 )
-     $   GO TO 20
-C
-   80 CONTINUE
-      DWORK(1) = WRKOPT
-      RETURN
-C *** Last line of AB08NX ***
-      END

--- a/extra/control-devel/src/AB09AD.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,363 +0,0 @@
-      SUBROUTINE AB09AD( DICO, JOB, EQUIL, ORDSEL, N, M, P, NR, A, LDA,
-     $                   B, LDB, C, LDC, HSV, TOL, IWORK, DWORK, LDWORK,
-     $                   IWARN, 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 compute a reduced order model (Ar,Br,Cr) for a stable original
-C     state-space representation (A,B,C) by using either the square-root
-C     or the balancing-free square-root Balance & Truncate (B & T)
-C     model reduction method.
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     DICO    CHARACTER*1
-C             Specifies the type of the original system as follows:
-C             = 'C':  continuous-time system;
-C             = 'D':  discrete-time system.
-C
-C     JOB     CHARACTER*1
-C             Specifies the model reduction approach to be used
-C             as follows:
-C             = 'B':  use the square-root Balance & Truncate method;
-C             = 'N':  use the balancing-free square-root
-C                     Balance & Truncate method.
-C
-C     EQUIL   CHARACTER*1
-C             Specifies whether the user wishes to preliminarily
-C             equilibrate the triplet (A,B,C) as follows:
-C             = 'S':  perform equilibration (scaling);
-C             = 'N':  do not perform equilibration.
-C
-C     ORDSEL  CHARACTER*1
-C             Specifies the order selection method as follows:
-C             = 'F':  the resulting order NR is fixed;
-C             = 'A':  the resulting order NR is automatically determined
-C                     on basis of the given tolerance TOL.
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The order of the original state-space representation, i.e.
-C             the order of the matrix A.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of system outputs.  P >= 0.
-C
-C     NR      (input/output) INTEGER
-C             On entry with ORDSEL = 'F', NR is the desired order of the
-C             resulting reduced order system.  0 <= NR <= N.
-C             On exit, if INFO = 0, NR is the order of the resulting
-C             reduced order model. NR is set as follows:
-C             if ORDSEL = 'F', NR is equal to MIN(NR,NMIN), where NR
-C             is the desired order on entry and NMIN is the order of a
-C             minimal realization of the given system; NMIN is
-C             determined as the number of Hankel singular values greater
-C             than N*EPS*HNORM(A,B,C), where EPS is the machine
-C             precision (see LAPACK Library Routine DLAMCH) and
-C             HNORM(A,B,C) is the Hankel norm of the system (computed
-C             in HSV(1));
-C             if ORDSEL = 'A', NR is equal to the number of Hankel
-C             singular values greater than MAX(TOL,N*EPS*HNORM(A,B,C)).
-C
-C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C             On entry, the leading N-by-N part of this array must
-C             contain the state dynamics matrix A.
-C             On exit, if INFO = 0, the leading NR-by-NR part of this
-C             array contains the state dynamics matrix Ar of the reduced
-C             order system.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N).
-C
-C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C             On entry, the leading N-by-M part of this array must
-C             contain the original input/state matrix B.
-C             On exit, if INFO = 0, the leading NR-by-M part of this
-C             array contains the input/state matrix Br of the reduced
-C             order system.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N).
-C
-C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C             On entry, the leading P-by-N part of this array must
-C             contain the original state/output matrix C.
-C             On exit, if INFO = 0, the leading P-by-NR part of this
-C             array contains the state/output matrix Cr of the reduced
-C             order system.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.  LDC >= MAX(1,P).
-C
-C     HSV     (output) DOUBLE PRECISION array, dimension (N)
-C             If INFO = 0, it contains the Hankel singular values of
-C             the original system ordered decreasingly. HSV(1) is the
-C             Hankel norm of the system.
-C
-C     Tolerances
-C
-C     TOL     DOUBLE PRECISION
-C             If ORDSEL = 'A', TOL contains the tolerance for
-C             determining the order of reduced system.
-C             For model reduction, the recommended value is
-C             TOL = c*HNORM(A,B,C), where c is a constant in the
-C             interval [0.00001,0.001], and HNORM(A,B,C) is the
-C             Hankel-norm of the given system (computed in HSV(1)).
-C             For computing a minimal realization, the recommended
-C             value is TOL = N*EPS*HNORM(A,B,C), where EPS is the
-C             machine precision (see LAPACK Library Routine DLAMCH).
-C             This value is used by default if TOL <= 0 on entry.
-C             If ORDSEL = 'F', the value of TOL is ignored.
-C
-C     Workspace
-C
-C     IWORK   INTEGER array, dimension (LIWORK)
-C             LIWORK = 0, if JOB = 'B';
-C             LIWORK = N, if JOB = 'N'.
-C
-C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
-C             On exit, if INFO = 0, DWORK(1) returns the optimal value
-C             of LDWORK.
-C
-C     LDWORK  INTEGER
-C             The length of the array DWORK.
-C             LDWORK >= MAX(1,N*(2*N+MAX(N,M,P)+5)+N*(N+1)/2).
-C             For optimum performance LDWORK should be larger.
-C
-C     Warning Indicator
-C
-C     IWARN   INTEGER
-C             = 0:  no warning;
-C             = 1:  with ORDSEL = 'F', the selected order NR is greater
-C                   than the order of a minimal realization of the
-C                   given system. In this case, the resulting NR is
-C                   set automatically to a value corresponding to the
-C                   order of a minimal realization of the system.
-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             = 1:  the reduction of A to the real Schur form failed;
-C             = 2:  the state matrix A is not stable (if DICO = 'C')
-C                   or not convergent (if DICO = 'D');
-C             = 3:  the computation of Hankel singular values failed.
-C
-C     METHOD
-C
-C     Let be the stable linear system
-C
-C          d[x(t)] = Ax(t) + Bu(t)
-C          y(t)    = Cx(t)                               (1)
-C
-C     where d[x(t)] is dx(t)/dt for a continuous-time system and x(t+1)
-C     for a discrete-time system. The subroutine AB09AD determines for
-C     the given system (1), the matrices of a reduced order system
-C
-C          d[z(t)] = Ar*z(t) + Br*u(t)
-C          yr(t)   = Cr*z(t)                             (2)
-C
-C     such that
-C
-C           HSV(NR) <= INFNORM(G-Gr) <= 2*[HSV(NR+1) + ... + HSV(N)],
-C
-C     where G and Gr are transfer-function matrices of the systems
-C     (A,B,C) and (Ar,Br,Cr), respectively, and INFNORM(G) is the
-C     infinity-norm of G.
-C
-C     If JOB = 'B', the square-root Balance & Truncate method of [1]
-C     is used and, for DICO = 'C', the resulting model is balanced.
-C     By setting TOL <= 0, the routine can be used to compute balanced
-C     minimal state-space realizations of stable systems.
-C
-C     If JOB = 'N', the balancing-free square-root version of the
-C     Balance & Truncate method [2] is used.
-C     By setting TOL <= 0, the routine can be used to compute minimal
-C     state-space realizations of stable systems.
-C
-C     REFERENCES
-C
-C     [1] Tombs M.S. and Postlethwaite I.
-C         Truncated balanced realization of stable, non-minimal
-C         state-space systems.
-C         Int. J. Control, Vol. 46, pp. 1319-1330, 1987.
-C
-C     [2] Varga A.
-C         Efficient minimal realization procedure based on balancing.
-C         Proc. of IMACS/IFAC Symp. MCTS, Lille, France, May 1991,
-C         A. El Moudui, P. Borne, S. G. Tzafestas (Eds.),
-C         Vol. 2, pp. 42-46.
-C
-C     NUMERICAL ASPECTS
-C
-C     The implemented methods rely on accuracy enhancing square-root or
-C     balancing-free square-root techniques.
-C                                         3
-C     The algorithms require less than 30N  floating point operations.
-C
-C     CONTRIBUTOR
-C
-C     C. Oara and A. Varga, German Aerospace Center,
-C     DLR Oberpfaffenhofen, March 1998.
-C     Based on the RASP routines SRBT and SRBFT.
-C
-C     REVISIONS
-C
-C     May 2, 1998.
-C     November 11, 1998, V. Sima, Research Institute for Informatics,
-C     Bucharest.
-C
-C     KEYWORDS
-C
-C     Balancing, minimal state-space representation, model reduction,
-C     multivariable system, state-space model.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ONE, C100
-      PARAMETER         ( ONE = 1.0D0, C100 = 100.0D0 )
-C     .. Scalar Arguments ..
-      CHARACTER         DICO, EQUIL, JOB, ORDSEL
-      INTEGER           INFO, IWARN, LDA, LDB, LDC, LDWORK, M, N, NR, P
-      DOUBLE PRECISION  TOL
-C     .. Array Arguments ..
-      INTEGER           IWORK(*)
-      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), HSV(*)
-C     .. Local Scalars ..
-      LOGICAL           FIXORD
-      INTEGER           IERR, KI, KR, KT, KTI, KW, NN
-      DOUBLE PRECISION  MAXRED, WRKOPT
-C     .. External Functions ..
-      LOGICAL           LSAME
-      EXTERNAL          LSAME
-C     .. External Subroutines ..
-      EXTERNAL          AB09AX, TB01ID, TB01WD, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         DBLE, MAX, MIN
-C     .. Executable Statements ..
-C
-      INFO   = 0
-      IWARN  = 0
-      FIXORD = LSAME( ORDSEL, 'F' )
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT. ( LSAME( DICO, 'C' ) .OR. LSAME( DICO, 'D' ) ) ) THEN
-         INFO = -1
-      ELSE IF( .NOT. ( LSAME( JOB, 'B' ) .OR. LSAME( JOB, 'N' ) ) ) THEN
-         INFO = -2
-      ELSE IF( .NOT. ( LSAME( EQUIL, 'S' ) .OR.
-     $                 LSAME( EQUIL, 'N' ) ) ) THEN
-         INFO = -3
-      ELSE IF( .NOT. ( FIXORD .OR. LSAME( ORDSEL, 'A' ) ) ) THEN
-         INFO = -4
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -5
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -6
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -7
-      ELSE IF( FIXORD .AND. ( NR.LT.0 .OR. NR.GT.N ) ) THEN
-         INFO = -8
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -10
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -12
-      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
-         INFO = -14
-      ELSE IF( LDWORK.LT.MAX( 1, N*( 2*N + MAX( N, M, P ) + 5 ) +
-     $                         ( N*( N + 1 ) )/2 ) ) THEN
-         INFO = -19
-      END IF
-C
-      IF( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB09AD', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF( MIN( N, M, P ).EQ.0 .OR. ( FIXORD .AND. NR.EQ.0 ) ) THEN
-         NR = 0
-         DWORK(1) = ONE
-         RETURN
-      END IF
-C
-C     Allocate working storage.
-C
-      NN = N*N
-      KT = 1
-      KR = KT + NN
-      KI = KR + N
-      KW = KI + N
-C
-      IF( LSAME( EQUIL, 'S' ) ) THEN
-C
-C        Scale simultaneously the matrices A, B and C:
-C        A <- inv(D)*A*D,  B <- inv(D)*B and C <- C*D, where D is a
-C        diagonal matrix.
-C
-         MAXRED = C100
-         CALL TB01ID( 'A', N, M, P, MAXRED, A, LDA, B, LDB, C, LDC,
-     $                DWORK, INFO )
-      END IF
-C
-C     Reduce A to the real Schur form using an orthogonal similarity
-C     transformation A <- T'*A*T and apply the transformation to
-C     B and C: B <- T'*B and C <- C*T.
-C
-      CALL TB01WD( N, M, P, A, LDA, B, LDB, C, LDC, DWORK(KT), N,
-     $             DWORK(KR), DWORK(KI), DWORK(KW), LDWORK-KW+1, IERR )
-      IF( IERR.NE.0 ) THEN
-         INFO = 1
-         RETURN
-      END IF
-C
-      WRKOPT = DWORK(KW) + DBLE( KW-1 )
-      KTI = KT  + NN
-      KW  = KTI + NN
-C
-      CALL AB09AX( DICO, JOB, ORDSEL, N, M, P, NR, A, LDA, B, LDB, C,
-     $             LDC, HSV, DWORK(KT), N, DWORK(KTI), N, TOL, IWORK,
-     $             DWORK(KW), LDWORK-KW+1, IWARN, IERR )
-C
-      IF( IERR.NE.0 ) THEN
-         INFO = IERR + 1
-         RETURN
-      END IF
-C
-      DWORK(1) = MAX( WRKOPT, DWORK(KW) + DBLE( KW-1 ) )
-C
-      RETURN
-C *** Last line of AB09AD ***
-      END

--- a/extra/control-devel/src/AB09AX.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,564 +0,0 @@
-      SUBROUTINE AB09AX( DICO, JOB, ORDSEL, N, M, P, NR, A, LDA, B, LDB,
-     $                   C, LDC, HSV, T, LDT, TI, LDTI, TOL, IWORK,
-     $                   DWORK, LDWORK, IWARN, 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 compute a reduced order model (Ar,Br,Cr) for a stable original
-C     state-space representation (A,B,C) by using either the square-root
-C     or the balancing-free square-root Balance & Truncate model
-C     reduction method. The state dynamics matrix A of the original
-C     system is an upper quasi-triangular matrix in real Schur canonical
-C     form. The matrices of the reduced order system are computed using
-C     the truncation formulas:
-C
-C          Ar = TI * A * T ,  Br = TI * B ,  Cr = C * T .
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     DICO    CHARACTER*1
-C             Specifies the type of the original system as follows:
-C             = 'C':  continuous-time system;
-C             = 'D':  discrete-time system.
-C
-C     JOB     CHARACTER*1
-C             Specifies the model reduction approach to be used
-C             as follows:
-C             = 'B':  use the square-root Balance & Truncate method;
-C             = 'N':  use the balancing-free square-root
-C                     Balance & Truncate method.
-C
-C     ORDSEL  CHARACTER*1
-C             Specifies the order selection method as follows:
-C             = 'F':  the resulting order NR is fixed;
-C             = 'A':  the resulting order NR is automatically determined
-C                     on basis of the given tolerance TOL.
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The order of the original state-space representation, i.e.
-C             the order of the matrix A.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of system outputs.  P >= 0.
-C
-C     NR      (input/output) INTEGER
-C             On entry with ORDSEL = 'F', NR is the desired order of the
-C             resulting reduced order system.  0 <= NR <= N.
-C             On exit, if INFO = 0, NR is the order of the resulting
-C             reduced order model. NR is set as follows:
-C             if ORDSEL = 'F', NR is equal to MIN(NR,NMIN), where NR
-C             is the desired order on entry and NMIN is the order of a
-C             minimal realization of the given system; NMIN is
-C             determined as the number of Hankel singular values greater
-C             than N*EPS*HNORM(A,B,C), where EPS is the machine
-C             precision (see LAPACK Library Routine DLAMCH) and
-C             HNORM(A,B,C) is the Hankel norm of the system (computed
-C             in HSV(1));
-C             if ORDSEL = 'A', NR is equal to the number of Hankel
-C             singular values greater than MAX(TOL,N*EPS*HNORM(A,B,C)).
-C
-C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C             On entry, the leading N-by-N part of this array must
-C             contain the state dynamics matrix A in a real Schur
-C             canonical form.
-C             On exit, if INFO = 0, the leading NR-by-NR part of this
-C             array contains the state dynamics matrix Ar of the
-C             reduced order system.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N).
-C
-C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C             On entry, the leading N-by-M part of this array must
-C             contain the original input/state matrix B.
-C             On exit, if INFO = 0, the leading NR-by-M part of this
-C             array contains the input/state matrix Br of the reduced
-C             order system.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N).
-C
-C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C             On entry, the leading P-by-N part of this array must
-C             contain the original state/output matrix C.
-C             On exit, if INFO = 0, the leading P-by-NR part of this
-C             array contains the state/output matrix Cr of the reduced
-C             order system.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.  LDC >= MAX(1,P).
-C
-C     HSV     (output) DOUBLE PRECISION array, dimension (N)
-C             If INFO = 0, it contains the Hankel singular values of
-C             the original system ordered decreasingly. HSV(1) is the
-C             Hankel norm of the system.
-C
-C     T       (output) DOUBLE PRECISION array, dimension (LDT,N)
-C             If INFO = 0 and NR > 0, the leading N-by-NR part of this
-C             array contains the right truncation matrix T.
-C
-C     LDT     INTEGER
-C             The leading dimension of array T.  LDT >= MAX(1,N).
-C
-C     TI      (output) DOUBLE PRECISION array, dimension (LDTI,N)
-C             If INFO = 0 and NR > 0, the leading NR-by-N part of this
-C             array contains the left truncation matrix TI.
-C
-C     LDTI    INTEGER
-C             The leading dimension of array TI.  LDTI >= MAX(1,N).
-C
-C     Tolerances
-C
-C     TOL     DOUBLE PRECISION
-C             If ORDSEL = 'A', TOL contains the tolerance for
-C             determining the order of reduced system.
-C             For model reduction, the recommended value is
-C             TOL = c*HNORM(A,B,C), where c is a constant in the
-C             interval [0.00001,0.001], and HNORM(A,B,C) is the
-C             Hankel-norm of the given system (computed in HSV(1)).
-C             For computing a minimal realization, the recommended
-C             value is TOL = N*EPS*HNORM(A,B,C), where EPS is the
-C             machine precision (see LAPACK Library Routine DLAMCH).
-C             This value is used by default if TOL <= 0 on entry.
-C             If ORDSEL = 'F', the value of TOL is ignored.
-C
-C     Workspace
-C
-C     IWORK   INTEGER array, dimension (LIWORK)
-C             LIWORK = 0, if JOB = 'B', or
-C             LIWORK = N, if JOB = 'N'.
-C
-C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
-C             On exit, if INFO = 0, DWORK(1) returns the optimal value
-C             of LDWORK.
-C
-C     LDWORK  INTEGER
-C             The length of the array DWORK.
-C             LDWORK >= MAX(1,N*(MAX(N,M,P)+5) + N*(N+1)/2).
-C             For optimum performance LDWORK should be larger.
-C
-C     Warning Indicator
-C
-C     IWARN   INTEGER
-C             = 0:  no warning;
-C             = 1:  with ORDSEL = 'F', the selected order NR is greater
-C                   than the order of a minimal realization of the
-C                   given system. In this case, the resulting NR is
-C                   set automatically to a value corresponding to the
-C                   order of a minimal realization of the system.
-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             = 1:  the state matrix A is not stable (if DICO = 'C')
-C                   or not convergent (if DICO = 'D');
-C             = 2:  the computation of Hankel singular values failed.
-C
-C     METHOD
-C
-C     Let be the stable linear system
-C
-C          d[x(t)] = Ax(t) + Bu(t)
-C          y(t)    = Cx(t)                               (1)
-C
-C     where d[x(t)] is dx(t)/dt for a continuous-time system and x(t+1)
-C     for a discrete-time system. The subroutine AB09AX determines for
-C     the given system (1), the matrices of a reduced NR order system
-C
-C          d[z(t)] = Ar*z(t) + Br*u(t)
-C          yr(t)   = Cr*z(t)                             (2)
-C
-C     such that
-C
-C           HSV(NR) <= INFNORM(G-Gr) <= 2*[HSV(NR+1) + ... + HSV(N)],
-C
-C     where G and Gr are transfer-function matrices of the systems
-C     (A,B,C) and (Ar,Br,Cr), respectively, and INFNORM(G) is the
-C     infinity-norm of G.
-C
-C     If JOB = 'B', the square-root Balance & Truncate method of [1]
-C     is used and, for DICO = 'C', the resulting model is balanced.
-C     By setting TOL <= 0, the routine can be used to compute balanced
-C     minimal state-space realizations of stable systems.
-C
-C     If JOB = 'N', the balancing-free square-root version of the
-C     Balance & Truncate method [2] is used.
-C     By setting TOL <= 0, the routine can be used to compute minimal
-C     state-space realizations of stable systems.
-C
-C     REFERENCES
-C
-C     [1] Tombs M.S. and Postlethwaite I.
-C         Truncated balanced realization of stable, non-minimal
-C         state-space systems.
-C         Int. J. Control, Vol. 46, pp. 1319-1330, 1987.
-C
-C     [2] Varga A.
-C         Efficient minimal realization procedure based on balancing.
-C         Proc. of IMACS/IFAC Symp. MCTS, Lille, France, May 1991,
-C         A. El Moudui, P. Borne, S. G. Tzafestas (Eds.),
-C         Vol. 2, pp. 42-46.
-C
-C     NUMERICAL ASPECTS
-C
-C     The implemented methods rely on accuracy enhancing square-root or
-C     balancing-free square-root techniques.
-C                                         3
-C     The algorithms require less than 30N  floating point operations.
-C
-C     CONTRIBUTOR
-C
-C     A. Varga, German Aerospace Center,
-C     DLR Oberpfaffenhofen, March 1998.
-C     Based on the RASP routines SRBT1 and SRBFT1.
-C
-C     REVISIONS
-C
-C     May 2, 1998.
-C     November 11, 1998, V. Sima, Research Institute for Informatics,
-C     Bucharest.
-C     December 1998, V. Sima, Katholieke Univ. Leuven, Leuven.
-C     February 14, 1999, A. Varga, German Aerospace Center.
-C     February 22, 1999, V. Sima, Research Institute for Informatics.
-C     February 27, 2000, V. Sima, Research Institute for Informatics.
-C
-C     KEYWORDS
-C
-C     Balancing, minimal state-space representation, model reduction,
-C     multivariable system, state-space model.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ONE, ZERO
-      PARAMETER         ( ONE = 1.0D0, ZERO = 0.0D0 )
-C     .. Scalar Arguments ..
-      CHARACTER         DICO, JOB, ORDSEL
-      INTEGER           INFO, IWARN, LDA, LDB, LDC, LDT, LDTI, LDWORK,
-     $                  M, N, NR, P
-      DOUBLE PRECISION  TOL
-C     .. Array Arguments ..
-      INTEGER           IWORK(*)
-      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), HSV(*),
-     $                  T(LDT,*), TI(LDTI,*)
-C     .. Local Scalars ..
-      LOGICAL           BAL, DISCR, FIXORD, PACKED
-      INTEGER           IERR, IJ, J, K, KTAU, KU, KV, KW, LDW, WRKOPT
-      DOUBLE PRECISION  ATOL, RTOL, SCALEC, SCALEO, TEMP
-C     .. External Functions ..
-      LOGICAL           LSAME
-      DOUBLE PRECISION  DLAMCH
-      EXTERNAL          DLAMCH, LSAME
-C     .. External Subroutines ..
-      EXTERNAL          DGEMM, DGEMV, DGEQRF, DGETRF, DGETRS, DLACPY,
-     $                  DORGQR, DSCAL, DTPMV, DTRMM, DTRMV, MA02AD,
-     $                  MA02DD, MB03UD, SB03OU, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         DBLE, INT, MAX, MIN, SQRT
-C     .. Executable Statements ..
-C
-      INFO   = 0
-      IWARN  = 0
-      DISCR  = LSAME( DICO,   'D' )
-      BAL    = LSAME( JOB,    'B' )
-      FIXORD = LSAME( ORDSEL, 'F' )
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT. ( LSAME( DICO, 'C' ) .OR. DISCR ) ) THEN
-         INFO = -1
-      ELSE IF( .NOT. ( BAL .OR. LSAME( JOB, 'N') ) ) THEN
-         INFO = -2
-      ELSE IF( .NOT. ( FIXORD .OR. LSAME( ORDSEL, 'A' ) ) ) THEN
-         INFO = -3
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -4
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -5
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -6
-      ELSE IF( FIXORD .AND. ( NR.LT.0 .OR. NR.GT.N ) ) THEN
-         INFO = -7
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -9
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -11
-      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
-         INFO = -13
-      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
-         INFO = -16
-      ELSE IF( LDTI.LT.MAX( 1, N ) ) THEN
-         INFO = -18
-      ELSE IF( LDWORK.LT.MAX( 1, N*( MAX( N, M, P ) + 5 ) +
-     $                         ( N*( N + 1 ) )/2 ) ) THEN
-         INFO = -22
-      END IF
-C
-      IF( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB09AX', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF( MIN( N, M, P ).EQ.0 .OR. ( FIXORD .AND. NR.EQ.0 ) ) THEN
-         NR = 0
-         DWORK(1) = ONE
-         RETURN
-      END IF
-C
-      RTOL = DBLE( N )*DLAMCH( 'Epsilon' )
-C
-C     Allocate N*MAX(N,M,P) and N working storage for the matrices U
-C     and TAU, respectively.
-C
-      KU   = 1
-      KTAU = KU + N*MAX( N, M, P )
-      KW   = KTAU + N
-      LDW  = LDWORK - KW + 1
-C
-C     Copy B in U.
-C
-      CALL DLACPY( 'Full', N, M, B, LDB, DWORK(KU), N )
-C
-C     If DISCR = .FALSE., solve for Su the Lyapunov equation
-C                                      2
-C     A*(Su*Su') + (Su*Su')*A' + scalec *B*B' = 0 .
-C
-C     If DISCR = .TRUE., solve for Su the Lyapunov equation
-C                           2
-C     A*(Su*Su')*A' + scalec *B*B' = Su*Su' .
-C
-C     Workspace:  need   N*(MAX(N,M,P) + 5);
-C                 prefer larger.
-C
-      CALL SB03OU( DISCR, .TRUE., N, M, A, LDA, DWORK(KU), N,
-     $             DWORK(KTAU), TI, LDTI, SCALEC, DWORK(KW), LDW, IERR )
-      IF( IERR.NE.0 ) THEN
-         INFO = 1
-         RETURN
-      ENDIF
-      WRKOPT = INT( DWORK(KW) ) + KW - 1
-C
-C     Copy C in U.
-C
-      CALL DLACPY( 'Full', P, N, C, LDC, DWORK(KU), P )
-C
-C     If DISCR = .FALSE., solve for Ru the Lyapunov equation
-C                                      2
-C     A'*(Ru'*Ru) + (Ru'*Ru)*A + scaleo  * C'*C = 0 .
-C
-C     If DISCR = .TRUE., solve for Ru the Lyapunov equation
-C                           2
-C     A'*(Ru'*Ru)*A + scaleo  * C'*C = Ru'*Ru .
-C
-C     Workspace:  need   N*(MAX(N,M,P) + 5);
-C                 prefer larger.
-C
-      CALL SB03OU( DISCR, .FALSE., N, P, A, LDA, DWORK(KU), P,
-     $             DWORK(KTAU), T, LDT, SCALEO, DWORK(KW), LDW, IERR )
-      WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
-C
-C     Allocate N*(N+1)/2 (or, if possible, N*N) working storage for the
-C     matrix V, a packed (or unpacked) copy of Su, and save Su in V.
-C     (The locations for TAU are reused here.)
-C
-      KV = KTAU
-      IF ( LDWORK-KV+1.LT.N*( N + 5 ) ) THEN
-         PACKED = .TRUE.
-         CALL MA02DD( 'Pack', 'Upper', N, TI, LDTI, DWORK(KV) )
-         KW = KV + ( N*( N + 1 ) )/2
-      ELSE
-         PACKED = .FALSE.
-         CALL DLACPY( 'Upper', N, N, TI, LDTI, DWORK(KV), N )
-         KW = KV + N*N
-      END IF
-C                               | x x |
-C     Compute Ru*Su in the form | 0 x | in TI.
-C
-      DO 10 J = 1, N
-         CALL DTRMV( 'Upper', 'NoTranspose', 'NonUnit', J, T, LDT,
-     $               TI(1,J), 1 )
-   10 CONTINUE
-C
-C     Compute the singular value decomposition Ru*Su = V*S*UT
-C     of the upper triangular matrix Ru*Su, with UT in TI and V in U.
-C
-C     Workspace:  need   N*MAX(N,M,P) + N*(N+1)/2 + 5*N;
-C                 prefer larger.
-C
-      CALL MB03UD( 'Vectors', 'Vectors', N, TI, LDTI, DWORK(KU), N, HSV,
-     $             DWORK(KW), LDWORK-KW+1, IERR )
-      IF( IERR.NE.0 ) THEN
-         INFO = 2
-         RETURN
-      ENDIF
-      WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
-C
-C     Scale singular values.
-C
-      CALL DSCAL( N, ONE / SCALEC / SCALEO, HSV, 1 )
-C
-C     Partition S, U and V conformally as:
-C
-C     S = diag(S1,S2),  U = [U1,U2] (U' in TI) and V = [V1,V2] (in U).
-C
-C     Compute the order of reduced system, as the order of S1.
-C
-      ATOL = RTOL*HSV(1)
-      IF( FIXORD ) THEN
-         IF( NR.GT.0 ) THEN
-            IF( HSV(NR).LE.ATOL ) THEN
-               NR = 0
-               IWARN = 1
-               FIXORD = .FALSE.
-            ENDIF
-         ENDIF
-      ELSE
-         ATOL = MAX( TOL, ATOL )
-         NR = 0
-      ENDIF
-      IF( .NOT.FIXORD ) THEN
-         DO 20 J = 1, N
-            IF( HSV(J).LE.ATOL ) GO TO 30
-            NR = NR + 1
-   20    CONTINUE
-   30    CONTINUE
-      ENDIF
-C
-      IF( NR.EQ.0 ) THEN
-         DWORK(1) = WRKOPT
-         RETURN
-      END IF
-C
-C     Compute the truncation matrices.
-C
-C     Compute TI' =  Ru'*V1 in U.
-C
-      CALL DTRMM( 'Left', 'Upper', 'Transpose', 'NonUnit', N, NR, ONE,
-     $             T, LDT, DWORK(KU), N )
-C
-C     Compute T = Su*U1 (with Su packed, if not enough workspace).
-C
-      CALL MA02AD( 'Full', NR, N, TI, LDTI, T, LDT )
-      IF ( PACKED ) THEN
-         DO 40 J = 1, NR
-            CALL DTPMV( 'Upper', 'NoTranspose', 'NonUnit', N, DWORK(KV),
-     $                  T(1,J), 1 )
-   40    CONTINUE
-      ELSE
-         CALL DTRMM( 'Left', 'Upper', 'NoTranspose', 'NonUnit', N, NR,
-     $               ONE, DWORK(KV), N, T, LDT )
-      END IF
-C
-      IF( BAL ) THEN
-         IJ = KU
-C
-C        Square-Root B & T method.
-C
-C        Compute the truncation matrices for balancing
-C                    -1/2           -1/2
-C                T*S1     and TI'*S1
-C
-         DO 50 J = 1, NR
-            TEMP = ONE/SQRT( HSV(J) )
-            CALL DSCAL( N, TEMP, T(1,J), 1 )
-            CALL DSCAL( N, TEMP, DWORK(IJ), 1 )
-            IJ = IJ + N
-   50    CONTINUE
-      ELSE
-C
-C        Balancing-Free B & T method.
-C
-C        Compute orthogonal bases for the images of matrices T and TI'.
-C
-C        Workspace:  need   N*MAX(N,M,P) + 2*NR;
-C                    prefer N*MAX(N,M,P) + NR*(NB+1)
-C                           (NB determined by ILAENV for DGEQRF).
-C
-         KW  = KTAU + NR
-         LDW = LDWORK - KW + 1
-         CALL DGEQRF( N, NR, T, LDT, DWORK(KTAU), DWORK(KW), LDW, IERR )
-         CALL DORGQR( N, NR, NR, T, LDT, DWORK(KTAU), DWORK(KW), LDW,
-     $                IERR )
-         CALL DGEQRF( N, NR, DWORK(KU), N, DWORK(KTAU), DWORK(KW), LDW,
-     $                IERR )
-         WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
-         CALL DORGQR( N, NR, NR, DWORK(KU), N, DWORK(KTAU), DWORK(KW),
-     $                LDW, IERR )
-         WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
-      END IF
-C
-C     Transpose TI' to obtain TI.
-C
-      CALL MA02AD( 'Full', N, NR, DWORK(KU), N, TI, LDTI )
-C
-      IF( .NOT.BAL ) THEN
-C                      -1
-C        Compute (TI*T)  *TI in TI.
-C
-         CALL DGEMM( 'NoTranspose', 'NoTranspose', NR, NR, N, ONE, TI,
-     $               LDTI, T, LDT, ZERO, DWORK(KU), N )
-         CALL DGETRF( NR, NR, DWORK(KU), N, IWORK, IERR )
-         CALL DGETRS( 'NoTranspose', NR, N, DWORK(KU), N, IWORK, TI,
-     $                LDTI, IERR )
-      END IF
-C
-C     Compute TI*A*T (A is in RSF).
-C
-      IJ = KU
-      DO 60 J = 1, N
-         K = MIN( J+1, N )
-         CALL DGEMV( 'NoTranspose', NR, K, ONE, TI, LDTI, A(1,J), 1,
-     $               ZERO, DWORK(IJ), 1 )
-         IJ = IJ + N
-   60 CONTINUE
-      CALL DGEMM( 'NoTranspose', 'NoTranspose', NR, NR, N, ONE,
-     $            DWORK(KU), N, T, LDT, ZERO, A, LDA )
-C
-C     Compute TI*B and C*T.
-C
-      CALL DLACPY( 'Full', N, M, B, LDB, DWORK(KU), N )
-      CALL DGEMM( 'NoTranspose', 'NoTranspose', NR, M, N, ONE, TI, LDTI,
-     $            DWORK(KU), N, ZERO, B, LDB )
-C
-      CALL DLACPY( 'Full', P, N, C, LDC, DWORK(KU), P )
-      CALL DGEMM( 'NoTranspose', 'NoTranspose', P, NR, N, ONE,
-     $            DWORK(KU), P, T, LDT, ZERO, C, LDC )
-C
-      DWORK(1) = WRKOPT
-C
-      RETURN
-C *** Last line of AB09AX ***
-      END

--- a/extra/control-devel/src/AB09BD.f	Wed Feb 22 15:13:45 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,385 +0,0 @@
-      SUBROUTINE AB09BD( DICO, JOB, EQUIL, ORDSEL, N, M, P, NR, A, LDA,
-     $                   B, LDB, C, LDC, D, LDD, HSV, TOL1, TOL2, IWORK,
-     $                   DWORK, LDWORK, IWARN, 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 compute a reduced order model (Ar,Br,Cr,Dr) for a stable
-C     original state-space representation (A,B,C,D) by using either the
-C     square-root or the balancing-free square-root Singular
-C     Perturbation Approximation (SPA) model reduction method.
-C
-C     ARGUMENTS
-C
-C     Mode Parameters
-C
-C     DICO    CHARACTER*1
-C             Specifies the type of the original system as follows:
-C             = 'C':  continuous-time system;
-C             = 'D':  discrete-time system.
-C
-C     JOB     CHARACTER*1
-C             Specifies the model reduction approach to be used
-C             as follows:
-C             = 'B':  use the square-root SPA method;
-C             = 'N':  use the balancing-free square-root SPA method.
-C
-C     EQUIL   CHARACTER*1
-C             Specifies whether the user wishes to preliminarily
-C             equilibrate the triplet (A,B,C) as follows:
-C             = 'S':  perform equilibration (scaling);
-C             = 'N':  do not perform equilibration.
-C
-C     ORDSEL  CHARACTER*1
-C             Specifies the order selection method as follows:
-C             = 'F':  the resulting order NR is fixed;
-C             = 'A':  the resulting order NR is automatically determined
-C                     on basis of the given tolerance TOL1.
-C
-C     Input/Output Parameters
-C
-C     N       (input) INTEGER
-C             The order of the original state-space representation, i.e.
-C             the order of the matrix A.  N >= 0.
-C
-C     M       (input) INTEGER
-C             The number of system inputs.  M >= 0.
-C
-C     P       (input) INTEGER
-C             The number of system outputs.  P >= 0.
-C
-C     NR      (input/output) INTEGER
-C             On entry with ORDSEL = 'F', NR is the desired order of
-C             the resulting reduced order system.  0 <= NR <= N.
-C             On exit, if INFO = 0, NR is the order of the resulting
-C             reduced order model. NR is set as follows:
-C             if ORDSEL = 'F', NR is equal to MIN(NR,NMIN), where NR
-C             is the desired order on entry and NMIN is the order of a
-C             minimal realization of the given system; NMIN is
-C             determined as the number of Hankel singular values greater
-C             than N*EPS*HNORM(A,B,C), where EPS is the machine
-C             precision (see LAPACK Library Routine DLAMCH) and
-C             HNORM(A,B,C) is the Hankel norm of the system (computed
-C             in HSV(1));
-C             if ORDSEL = 'A', NR is equal to the number of Hankel
-C             singular values greater than MAX(TOL1,N*EPS*HNORM(A,B,C)).
-C
-C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C             On entry, the leading N-by-N part of this array must
-C             contain the state dynamics matrix A.
-C             On exit, if INFO = 0, the leading NR-by-NR part of this
-C             array contains the state dynamics matrix Ar of the
-C             reduced order system.
-C
-C     LDA     INTEGER
-C             The leading dimension of array A.  LDA >= MAX(1,N).
-C
-C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C             On entry, the leading N-by-M part of this array must
-C             contain the original input/state matrix B.
-C             On exit, if INFO = 0, the leading NR-by-M part of this
-C             array contains the input/state matrix Br of the reduced
-C             order system.
-C
-C     LDB     INTEGER
-C             The leading dimension of array B.  LDB >= MAX(1,N).
-C
-C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C             On entry, the leading P-by-N part of this array must
-C             contain the original state/output matrix C.
-C             On exit, if INFO = 0, the leading P-by-NR part of this
-C             array contains the state/output matrix Cr of the reduced
-C             order system.
-C
-C     LDC     INTEGER
-C             The leading dimension of array C.  LDC >= MAX(1,P).
-C
-C     D       (input/output) DOUBLE PRECISION array, dimension (LDD,M)
-C             On entry, the leading P-by-M part of this array must
-C             contain the original input/output matrix D.
-C             On exit, if INFO = 0, the leading P-by-M part of this
-C             array contains the input/output matrix Dr of the reduced
-C             order system.
-C
-C     LDD     INTEGER
-C             The leading dimension of array D.  LDD >= MAX(1,P).
-C
-C     HSV     (output) DOUBLE PRECISION array, dimension (N)
-C             If INFO = 0, it contains the Hankel singular values of
-C             the original system ordered decreasingly. HSV(1) is the
-C             Hankel norm of the system.
-C
-C     Tolerances
-C
-C     TOL1    DOUBLE PRECISION
-C             If ORDSEL = 'A', TOL1 contains the tolerance for
-C             determining the order of reduced system.
-C             For model reduction, the recommended value is
-C             TOL1 = c*HNORM(A,B,C), where c is a constant in the
-C             interval [0.00001,0.001], and HNORM(A,B,C) is the
-C             Hankel-norm of the given system (computed in HSV(1)).
-C             For computing a minimal realization, the recommended
-C             value is TOL1 = N*EPS*HNORM(A,B,C), where EPS is the
-C             machine precision (see LAPACK Library Routine DLAMCH).
-C             This value is used by default if TOL1 <= 0 on entry.
-C             If ORDSEL = 'F', the value of TOL1 is ignored.
-C
-C     TOL2    DOUBLE PRECISION
-C             The tolerance for determining the order of a minimal
-C             realization of the given system. The recommended value is
-C             TOL2 = N*EPS*HNORM(A,B,C). This value is used by default
-C             if TOL2 <= 0 on entry.
-C             If TOL2 > 0, then TOL2 <= TOL1.
-C
-C     Workspace
-C
-C     IWORK   INTEGER array, dimension MAX(1,2*N)
-C             On exit with INFO = 0, IWORK(1) contains the order of the
-C             minimal realization of the system.
-C
-C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
-C             On exit, if INFO = 0, DWORK(1) returns the optimal value
-C             of LDWORK.
-C
-C     LDWORK  INTEGER
-C             The length of the array DWORK.
-C             LDWORK >= MAX(1,N*(2*N+MAX(N,M,P)+5)+N*(N+1)/2).
-C             For optimum performance LDWORK should be larger.
-C
-C     Warning Indicator
-C
-C     IWARN   INTEGER
-C             = 0:  no warning;
-C             = 1:  with ORDSEL = 'F', the selected order NR is greater
-C                   than the order of a minimal realization of the
-C                   given system. In this case, the resulting NR is
-C                   set automatically to a value corresponding to the
-C                   order of a minimal realization of the system.
-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             = 1:  the reduction of A to the real Schur form failed;
-C             = 2:  the state matrix A is not stable (if DICO = 'C')
-C                   or not convergent (if DICO = 'D');
-C             = 3:  the computation of Hankel singular values failed.
-C
-C     METHOD
-C
-C     Let be the stable linear system
-C
-C          d[x(t)] = Ax(t) + Bu(t)
-C          y(t)    = Cx(t) + Du(t)                           (1)
-C
-C     where d[x(t)] is dx(t)/dt for a continuous-time system and x(t+1)
-C     for a discrete-time system. The subroutine AB09BD determines for
-C     the given system (1), the matrices of a reduced order system
-C
-C          d[z(t)] = Ar*z(t) + Br*u(t)
-C          yr(t)   = Cr*z(t) + Dr*u(t)                       (2)
-C
-C     such that
-C
-C           HSV(NR) <= INFNORM(G-Gr) <= 2*[HSV(NR+1) + ... + HSV(N)],
-C
-C     where G and Gr are transfer-function matrices of the systems
-C     (A,B,C,D) and (Ar,Br,Cr,Dr), respectively, and INFNORM(G) is the
-C     infinity-norm of G.
-C
-C     If JOB = 'B', the balancing-based square-root SPA method of [1]
-C     is used and the resulting model is balanced.
-C
-C     If JOB = 'N', the balancing-free square-root SPA method of [2]
-C     is used.
-C     By setting TOL1 = TOL2, the routine can be used to compute
-C     Balance & Truncate approximations.
-C
-C     REFERENCES
-C
-C     [1] Liu Y. and Anderson B.D.O.
-C         Singular Perturbation Approximation of Balanced Systems,
-C         Int. J. Control, Vol. 50, pp. 1379-1405, 1989.
-C
-C     [2] Varga A.
-C         Balancing-free square-root algorithm for computing singular
-C         perturbation approximations.
-C         Proc. 30-th IEEE CDC,  Brighton, Dec. 11-13, 1991,
-C         Vol. 2, pp. 1062-1065.
-C
-C     NUMERICAL ASPECTS
-C
-C     The implemented methods rely on accuracy enhancing square-root or
-C     balancing-free square-root techniques.
-C                                         3
-C     The algorithms require less than 30N  floating point operations.
-C
-C     CONTRIBUTOR
-C
-C     C. Oara and A. Varga, German Aerospace Center,
-C     DLR Oberpfaffenhofen, March 1998.
-C     Based on the RASP routine SRBFSP.
-C
-C     REVISIONS
-C
-C     May 2, 1998.
-C     November 11, 1998, V. Sima, Research Institute for Informatics,
-C     Bucharest.
-C
-C     KEYWORDS
-C
-C     Balancing, minimal state-space representation, model reduction,
-C     multivariable system, singular perturbation approximation,
-C     state-space model.
-C
-C     ******************************************************************
-C
-C     .. Parameters ..
-      DOUBLE PRECISION  ZERO, ONE, C100
-      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0, C100 = 100.0D0 )
-C     .. Scalar Arguments ..
-      CHARACTER         DICO, EQUIL, JOB, ORDSEL
-      INTEGER           INFO, IWARN, LDA, LDB, LDC, LDD, LDWORK,
-     $                  M, N, NR, P
-      DOUBLE PRECISION  TOL1, TOL2
-C     .. Array Arguments ..
-      INTEGER           IWORK(*)
-      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
-     $                  DWORK(*), HSV(*)
-C     .. Local Scalars ..
-      LOGICAL           FIXORD
-      INTEGER           IERR, KI, KR, KT, KTI, KW, NN
-      DOUBLE PRECISION  MAXRED, WRKOPT
-C     .. External Functions ..
-      LOGICAL           LSAME
-      EXTERNAL          LSAME
-C     .. External Subroutines ..
-      EXTERNAL          AB09BX, TB01ID, TB01WD, XERBLA
-C     .. Intrinsic Functions ..
-      INTRINSIC         DBLE, MAX, MIN
-C     .. Executable Statements ..
-C
-      INFO   = 0
-      IWARN  = 0
-      FIXORD = LSAME( ORDSEL, 'F' )
-C
-C     Test the input scalar arguments.
-C
-      IF( .NOT. ( LSAME( DICO, 'C' ) .OR. LSAME( DICO, 'D' ) ) ) THEN
-         INFO = -1
-      ELSE IF( .NOT. ( LSAME( JOB, 'B' ) .OR. LSAME( JOB, 'N' ) ) ) THEN
-         INFO = -2
-      ELSE IF( .NOT. ( LSAME( EQUIL, 'S' ) .OR.
-     $                 LSAME( EQUIL, 'N' ) ) ) THEN
-         INFO = -3
-      ELSE IF( .NOT. ( FIXORD .OR. LSAME( ORDSEL, 'A' ) ) ) THEN
-         INFO = -4
-      ELSE IF( N.LT.0 ) THEN
-         INFO = -5
-      ELSE IF( M.LT.0 ) THEN
-         INFO = -6
-      ELSE IF( P.LT.0 ) THEN
-         INFO = -7
-      ELSE IF( FIXORD .AND. ( NR.LT.0 .OR. NR.GT.N ) ) THEN
-         INFO = -8
-      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
-         INFO = -10
-      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
-         INFO = -12
-      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
-         INFO = -14
-      ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
-         INFO = -16
-      ELSE IF( TOL2.GT.ZERO .AND. TOL2.GT.TOL1 ) THEN
-         INFO = -19
-      ELSE IF( LDWORK.LT.MAX( 1, N*( 2*N + MAX( N, M, P ) + 5 ) +
-     $                         ( N*( N + 1 ) )/2 ) ) THEN
-         INFO = -22
-      END IF
-C
-      IF( INFO.NE.0 ) THEN
-C
-C        Error return.
-C
-         CALL XERBLA( 'AB09BD', -INFO )
-         RETURN
-      END IF
-C
-C     Quick return if possible.
-C
-      IF( MIN( N, M, P ).EQ.0 ) THEN
-         NR = 0
-         IWORK(1) = 0
-         DWORK(1) = ONE
-         RETURN
-      END IF
-C
-C     Allocate working storage.
-C
-      NN = N*N
-      KT = 1
-      KR = KT + NN
-      KI = KR + N
-      KW = KI + N
-C
-      IF( LSAME( EQUIL, 'S' ) ) THEN
-C
-C        Scale simultaneously the matrices A, B and C:
-C        A <- inv(D)*A*D,  B <- inv(D)*B and C <- C*D, where D is a
-C        diagonal matrix.
-C
-         MAXRED = C100
-         CALL TB01ID( 'All', N, M, P, MAXRED, A, LDA, B, LDB, C, LDC,
-     $                DWORK, INFO )
-      END IF
-C
-C     Reduce A to the real Schur form using an orthogonal similarity
-C     transformation A <- T'*A*T and apply the transformation to
-C     B and C: B <- T'*B and C <- C*T.
-C
-      CALL TB01WD( N, M, P, A, LDA, B, LDB, C, LDC, DWORK(KT), N,
-     $             DWORK(KR), DWORK(KI), DWORK(KW), LDWORK-KW+1, IERR )
-      IF( IERR.NE.0 ) THEN
-         INFO = 1
-         RETURN
-      END IF
-C
-      WRKOPT = DWORK(KW) + DBLE( KW-1 )
-C
-      KTI = KT  + NN
-      KW  = KTI + NN
-      CALL AB09BX( DICO, JOB, ORDSEL, N, M, P, NR, A, LDA, B, LDB,
-     $             C, LDC, D, LDD, HSV, DWORK(KT), N, DWORK(KTI), N,
-     $             TOL1, TOL2, IWORK, DWORK(KW), LDWORK-KW+1, IWARN,
-     $             IERR )
-C
-      IF( IERR.NE.0 ) THEN
-         INFO = IERR + 1
-         RETURN
-      END IF
-C
-      DWORK(1) = MAX( WRKOPT, DWORK(KW) + DBLE( KW-1 ) )
-C
-      RETURN
-C *** Last line of AB09BD ***
-      END