**Purpose**

To compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using the frequency weighted optimal Hankel-norm approximation method. The Hankel norm of the weighted error V*(G-Gr)*W or conj(V)*(G-Gr)*conj(W) is minimized, where G and Gr are the transfer-function matrices of the original and reduced systems, respectively, and V and W are the transfer-function matrices of the left and right frequency weights, specified by their state space realizations (AV,BV,CV,DV) and (AW,BW,CW,DW), respectively. When minimizing the weighted error V*(G-Gr)*W, V and W must be antistable transfer-function matrices. When minimizing conj(V)*(G-Gr)*conj(W), V and W must be stable transfer-function matrices. Additionally, V and W must be invertible transfer-function matrices, with the feedthrough matrices DV and DW invertible. If the original system is unstable, then the frequency weighted Hankel-norm approximation is computed only for the ALPHA-stable part of the system. For a transfer-function matrix G, conj(G) denotes the conjugate of G given by G'(-s) for a continuous-time system or G'(1/z) for a discrete-time system.

SUBROUTINE AB09KD( JOB, DICO, WEIGHT, EQUIL, ORDSEL, N, NV, NW, M, $ P, NR, ALPHA, A, LDA, B, LDB, C, LDC, D, LDD, $ AV, LDAV, BV, LDBV, CV, LDCV, DV, LDDV, $ AW, LDAW, BW, LDBW, CW, LDCW, DW, LDDW, $ NS, HSV, TOL1, TOL2, IWORK, DWORK, LDWORK, $ IWARN, INFO ) C .. Scalar Arguments .. CHARACTER DICO, EQUIL, JOB, ORDSEL, WEIGHT INTEGER INFO, IWARN, LDA, LDAV, LDAW, LDB, LDBV, LDBW, $ LDC, LDCV, LDCW, LDD, LDDV, LDDW, LDWORK, M, N, $ NR, NS, NV, NW, P DOUBLE PRECISION ALPHA, TOL1, TOL2 C .. Array Arguments .. INTEGER IWORK(*) DOUBLE PRECISION A(LDA,*), AV(LDAV,*), AW(LDAW,*), $ B(LDB,*), BV(LDBV,*), BW(LDBW,*), $ C(LDC,*), CV(LDCV,*), CW(LDCW,*), $ D(LDD,*), DV(LDDV,*), DW(LDDW,*), DWORK(*), $ HSV(*)

**Mode Parameters**

JOB CHARACTER*1 Specifies the frequency-weighting problem as follows: = 'N': solve min||V*(G-Gr)*W||_H; = 'C': solve min||conj(V)*(G-Gr)*conj(W)||_H. DICO CHARACTER*1 Specifies the type of the original system as follows: = 'C': continuous-time system; = 'D': discrete-time system. WEIGHT CHARACTER*1 Specifies the type of frequency weighting, as follows: = 'N': no weightings are used (V = I, W = I); = 'L': only left weighting V is used (W = I); = 'R': only right weighting W is used (V = I); = 'B': both left and right weightings V and W are used. EQUIL CHARACTER*1 Specifies whether the user wishes to preliminarily equilibrate the triplet (A,B,C) as follows: = 'S': perform equilibration (scaling); = 'N': do not perform equilibration. ORDSEL CHARACTER*1 Specifies the order selection method as follows: = 'F': the resulting order NR is fixed; = 'A': the resulting order NR is automatically determined on basis of the given tolerance TOL1.

N (input) INTEGER The order of the original state-space representation, i.e., the order of the matrix A. N >= 0. NV (input) INTEGER The order of the realization of the left frequency weighting V, i.e., the order of the matrix AV. NV >= 0. NW (input) INTEGER The order of the realization of the right frequency weighting W, i.e., the order of the matrix AW. NW >= 0. M (input) INTEGER The number of system inputs. M >= 0. P (input) INTEGER The number of system outputs. P >= 0. NR (input/output) INTEGER On entry with ORDSEL = 'F', NR is the desired order of the resulting reduced order system. 0 <= NR <= N. On exit, if INFO = 0, NR is the order of the resulting reduced order model. For a system with NU ALPHA-unstable eigenvalues and NS ALPHA-stable eigenvalues (NU+NS = N), NR is set as follows: if ORDSEL = 'F', NR is equal to NU+MIN(MAX(0,NR-NU-KR+1),NMIN), where KR is the multiplicity of the Hankel singular value HSV(NR-NU+1), NR is the desired order on entry, and NMIN is the order of a minimal realization of the ALPHA-stable part of the given system; NMIN is determined as the number of Hankel singular values greater than NS*EPS*HNORM(As,Bs,Cs), where EPS is the machine precision (see LAPACK Library Routine DLAMCH) and HNORM(As,Bs,Cs) is the Hankel norm of the ALPHA-stable part of the weighted system (computed in HSV(1)); if ORDSEL = 'A', NR is the sum of NU and the number of Hankel singular values greater than MAX(TOL1,NS*EPS*HNORM(As,Bs,Cs)). ALPHA (input) DOUBLE PRECISION Specifies the ALPHA-stability boundary for the eigenvalues of the state dynamics matrix A. For a continuous-time system (DICO = 'C'), ALPHA <= 0 is the boundary value for the real parts of eigenvalues, while for a discrete-time system (DICO = 'D'), 0 <= ALPHA <= 1 represents the boundary value for the moduli of eigenvalues. The ALPHA-stability domain does not include the boundary. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading N-by-N part of this array must contain the state dynamics matrix A. On exit, if INFO = 0, the leading NR-by-NR part of this array contains the state dynamics matrix Ar of the reduced order system in a real Schur form. The resulting A has a block-diagonal form with two blocks. For a system with NU ALPHA-unstable eigenvalues and NS ALPHA-stable eigenvalues (NU+NS = N), the leading NU-by-NU block contains the unreduced part of A corresponding to ALPHA-unstable eigenvalues. The trailing (NR+NS-N)-by-(NR+NS-N) block contains the reduced part of A corresponding to ALPHA-stable eigenvalues. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB,M) On entry, the leading N-by-M part of this array must contain the original input/state matrix B. On exit, if INFO = 0, the leading NR-by-M part of this array contains the input/state matrix Br of the reduced order system. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N). C (input/output) DOUBLE PRECISION array, dimension (LDC,N) On entry, the leading P-by-N part of this array must contain the original state/output matrix C. On exit, if INFO = 0, the leading P-by-NR part of this array contains the state/output matrix Cr of the reduced order system. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,P). D (input/output) DOUBLE PRECISION array, dimension (LDD,M) On entry, the leading P-by-M part of this array must contain the original input/output matrix D. On exit, if INFO = 0, the leading P-by-M part of this array contains the input/output matrix Dr of the reduced order system. LDD INTEGER The leading dimension of array D. LDD >= MAX(1,P). AV (input/output) DOUBLE PRECISION array, dimension (LDAV,NV) On entry, if WEIGHT = 'L' or 'B', the leading NV-by-NV part of this array must contain the state matrix AV of a state space realization of the left frequency weighting V. On exit, if WEIGHT = 'L' or 'B', and INFO = 0, the leading NV-by-NV part of this array contains a real Schur form of the state matrix of a state space realization of the inverse of V. AV is not referenced if WEIGHT = 'R' or 'N'. LDAV INTEGER The leading dimension of the array AV. LDAV >= MAX(1,NV), if WEIGHT = 'L' or 'B'; LDAV >= 1, if WEIGHT = 'R' or 'N'. BV (input/output) DOUBLE PRECISION array, dimension (LDBV,P) On entry, if WEIGHT = 'L' or 'B', the leading NV-by-P part of this array must contain the input matrix BV of a state space realization of the left frequency weighting V. On exit, if WEIGHT = 'L' or 'B', and INFO = 0, the leading NV-by-P part of this array contains the input matrix of a state space realization of the inverse of V. BV is not referenced if WEIGHT = 'R' or 'N'. LDBV INTEGER The leading dimension of the array BV. LDBV >= MAX(1,NV), if WEIGHT = 'L' or 'B'; LDBV >= 1, if WEIGHT = 'R' or 'N'. CV (input/output) DOUBLE PRECISION array, dimension (LDCV,NV) On entry, if WEIGHT = 'L' or 'B', the leading P-by-NV part of this array must contain the output matrix CV of a state space realization of the left frequency weighting V. On exit, if WEIGHT = 'L' or 'B', and INFO = 0, the leading P-by-NV part of this array contains the output matrix of a state space realization of the inverse of V. CV is not referenced if WEIGHT = 'R' or 'N'. LDCV INTEGER The leading dimension of the array CV. LDCV >= MAX(1,P), if WEIGHT = 'L' or 'B'; LDCV >= 1, if WEIGHT = 'R' or 'N'. DV (input/output) DOUBLE PRECISION array, dimension (LDDV,P) On entry, if WEIGHT = 'L' or 'B', the leading P-by-P part of this array must contain the feedthrough matrix DV of a state space realization of the left frequency weighting V. On exit, if WEIGHT = 'L' or 'B', and INFO = 0, the leading P-by-P part of this array contains the feedthrough matrix of a state space realization of the inverse of V. DV is not referenced if WEIGHT = 'R' or 'N'. LDDV INTEGER The leading dimension of the array DV. LDDV >= MAX(1,P), if WEIGHT = 'L' or 'B'; LDDV >= 1, if WEIGHT = 'R' or 'N'. AW (input/output) DOUBLE PRECISION array, dimension (LDAW,NW) On entry, if WEIGHT = 'R' or 'B', the leading NW-by-NW part of this array must contain the state matrix AW of a state space realization of the right frequency weighting W. On exit, if WEIGHT = 'R' or 'B', and INFO = 0, the leading NW-by-NW part of this array contains a real Schur form of the state matrix of a state space realization of the inverse of W. AW is not referenced if WEIGHT = 'L' or 'N'. LDAW INTEGER The leading dimension of the array AW. LDAW >= MAX(1,NW), if WEIGHT = 'R' or 'B'; LDAW >= 1, if WEIGHT = 'L' or 'N'. BW (input/output) DOUBLE PRECISION array, dimension (LDBW,M) On entry, if WEIGHT = 'R' or 'B', the leading NW-by-M part of this array must contain the input matrix BW of a state space realization of the right frequency weighting W. On exit, if WEIGHT = 'R' or 'B', and INFO = 0, the leading NW-by-M part of this array contains the input matrix of a state space realization of the inverse of W. BW is not referenced if WEIGHT = 'L' or 'N'. LDBW INTEGER The leading dimension of the array BW. LDBW >= MAX(1,NW), if WEIGHT = 'R' or 'B'; LDBW >= 1, if WEIGHT = 'L' or 'N'. CW (input/output) DOUBLE PRECISION array, dimension (LDCW,NW) On entry, if WEIGHT = 'R' or 'B', the leading M-by-NW part of this array must contain the output matrix CW of a state space realization of the right frequency weighting W. On exit, if WEIGHT = 'R' or 'B', and INFO = 0, the leading M-by-NW part of this array contains the output matrix of a state space realization of the inverse of W. CW is not referenced if WEIGHT = 'L' or 'N'. LDCW INTEGER The leading dimension of the array CW. LDCW >= MAX(1,M), if WEIGHT = 'R' or 'B'; LDCW >= 1, if WEIGHT = 'L' or 'N'. DW (input/output) DOUBLE PRECISION array, dimension (LDDW,M) On entry, if WEIGHT = 'R' or 'B', the leading M-by-M part of this array must contain the feedthrough matrix DW of a state space realization of the right frequency weighting W. On exit, if WEIGHT = 'R' or 'B', and INFO = 0, the leading M-by-M part of this array contains the feedthrough matrix of a state space realization of the inverse of W. DW is not referenced if WEIGHT = 'L' or 'N'. LDDW INTEGER The leading dimension of the array DW. LDDW >= MAX(1,M), if WEIGHT = 'R' or 'B'; LDDW >= 1, if WEIGHT = 'L' or 'N'. NS (output) INTEGER The dimension of the ALPHA-stable subsystem. HSV (output) DOUBLE PRECISION array, dimension (N) If INFO = 0, the leading NS elements of this array contain the Hankel singular values, ordered decreasingly, of the ALPHA-stable part of the weighted original system. HSV(1) is the Hankel norm of the ALPHA-stable weighted subsystem.

TOL1 DOUBLE PRECISION If ORDSEL = 'A', TOL1 contains the tolerance for determining the order of reduced system. For model reduction, the recommended value is TOL1 = c*HNORM(As,Bs,Cs), where c is a constant in the interval [0.00001,0.001], and HNORM(As,Bs,Cs) is the Hankel-norm of the ALPHA-stable part of the weighted original system (computed in HSV(1)). If TOL1 <= 0 on entry, the used default value is TOL1 = NS*EPS*HNORM(As,Bs,Cs), where NS is the number of ALPHA-stable eigenvalues of A and EPS is the machine precision (see LAPACK Library Routine DLAMCH). If ORDSEL = 'F', the value of TOL1 is ignored. TOL2 DOUBLE PRECISION The tolerance for determining the order of a minimal realization of the ALPHA-stable part of the given system. The recommended value is TOL2 = NS*EPS*HNORM(As,Bs,Cs). This value is used by default if TOL2 <= 0 on entry. If TOL2 > 0 and ORDSEL = 'A', then TOL2 <= TOL1.

IWORK INTEGER array, dimension (LIWORK) LIWORK = MAX(1,M,c), if DICO = 'C', LIWORK = MAX(1,N,M,c), if DICO = 'D', where c = 0, if WEIGHT = 'N', c = 2*P, if WEIGHT = 'L', c = 2*M, if WEIGHT = 'R', c = MAX(2*M,2*P), if WEIGHT = 'B'. On exit, if INFO = 0, IWORK(1) contains NMIN, the order of the computed minimal realization. DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK. LDWORK INTEGER The length of the array DWORK. LDWORK >= MAX( LDW1, LDW2, LDW3, LDW4 ), where LDW1 = 0 if WEIGHT = 'R' or 'N' and LDW1 = MAX( NV*(NV+5), NV*N + MAX( a, P*N, P*M ) ) if WEIGHT = 'L' or WEIGHT = 'B', LDW2 = 0 if WEIGHT = 'L' or 'N' and LDW2 = MAX( NW*(NW+5), NW*N + MAX( b, M*N, P*M ) ) if WEIGHT = 'R' or WEIGHT = 'B', with a = 0, b = 0, if DICO = 'C' or JOB = 'N', a = 2*NV, b = 2*NW, if DICO = 'D' and JOB = 'C'; LDW3 = N*(2*N + MAX(N,M,P) + 5) + N*(N+1)/2, LDW4 = N*(M+P+2) + 2*M*P + MIN(N,M) + MAX( 3*M+1, MIN(N,M)+P ). For optimum performance LDWORK should be larger.

IWARN INTEGER = 0: no warning; = 1: with ORDSEL = 'F', the selected order NR is greater than NSMIN, the sum of the order of the ALPHA-unstable part and the order of a minimal realization of the ALPHA-stable part of the given system; in this case, the resulting NR is set equal to NSMIN; = 2: with ORDSEL = 'F', the selected order NR is less than the order of the ALPHA-unstable part of the given system; in this case NR is set equal to the order of the ALPHA-unstable part.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 1: the computation of the ordered real Schur form of A failed; = 2: the separation of the ALPHA-stable/unstable diagonal blocks failed because of very close eigenvalues; = 3: the reduction of AV or AV-BV*inv(DV)*CV to a real Schur form failed; = 4: the reduction of AW or AW-BW*inv(DW)*CW to a real Schur form failed; = 5: JOB = 'N' and AV is not antistable, or JOB = 'C' and AV is not stable; = 6: JOB = 'N' and AW is not antistable, or JOB = 'C' and AW is not stable; = 7: the computation of Hankel singular values failed; = 8: the computation of stable projection in the Hankel-norm approximation algorithm failed; = 9: the order of computed stable projection in the Hankel-norm approximation algorithm differs from the order of Hankel-norm approximation; = 10: DV is singular; = 11: DW is singular; = 12: the solution of the Sylvester equation failed because the zeros of V (if JOB = 'N') or of conj(V) (if JOB = 'C') are not distinct from the poles of G1sr (see METHOD); = 13: the solution of the Sylvester equation failed because the zeros of W (if JOB = 'N') or of conj(W) (if JOB = 'C') are not distinct from the poles of G1sr (see METHOD).

Let G be the transfer-function matrix of the original linear system d[x(t)] = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t), (1) where d[x(t)] is dx(t)/dt for a continuous-time system and x(t+1) for a discrete-time system. The subroutine AB09KD determines the matrices of a reduced order system d[z(t)] = Ar*z(t) + Br*u(t) yr(t) = Cr*z(t) + Dr*u(t), (2) such that the corresponding transfer-function matrix Gr minimizes the Hankel-norm of the frequency-weighted error V*(G-Gr)*W, (3) or conj(V)*(G-Gr)*conj(W). (4) For minimizing (3), V and W are assumed to be antistable, while for minimizing (4), V and W are assumed to be stable transfer- function matrices. Note: conj(G) = G'(-s) for a continuous-time system and conj(G) = G'(1/z) for a discrete-time system. The following procedure is used to reduce G (see [1]): 1) Decompose additively G as G = G1 + G2, such that G1 = (A1,B1,C1,D) has only ALPHA-stable poles and G2 = (A2,B2,C2,0) has only ALPHA-unstable poles. 2) Compute G1s, the stable projection of V*G1*W or conj(V)*G1*conj(W), using explicit formulas [4]. 3) Determine G1sr, the optimal Hankel-norm approximation of G1s of order r. 4) Compute G1r, the stable projection of either inv(V)*G1sr*inv(W) or conj(inv(V))*G1sr*conj(inv(W)), using explicit formulas [4]. 5) Assemble the reduced model Gr as Gr = G1r + G2. To reduce the weighted ALPHA-stable part G1s at step 3, the optimal Hankel-norm approximation method of [2], based on the square-root balancing projection formulas of [3], is employed. The optimal weighted approximation error satisfies HNORM[V*(G-Gr)*W] = S(r+1), or HNORM[conj(V)*(G-Gr)*conj(W)] = S(r+1), where S(r+1) is the (r+1)-th Hankel singular value of G1s, the transfer-function matrix computed at step 2 of the above procedure, and HNORM(.) denotes the Hankel-norm.

[1] Latham, G.A. and Anderson, B.D.O. Frequency-weighted optimal Hankel-norm approximation of stable transfer functions. Systems & Control Letters, Vol. 5, pp. 229-236, 1985. [2] Glover, K. All optimal Hankel norm approximation of linear multivariable systems and their L-infinity error bounds. Int. J. Control, Vol. 36, pp. 1145-1193, 1984. [3] Tombs M.S. and Postlethwaite I. Truncated balanced realization of stable, non-minimal state-space systems. Int. J. Control, Vol. 46, pp. 1319-1330, 1987. [4] Varga A. Explicit formulas for an efficient implementation of the frequency-weighting model reduction approach. Proc. 1993 European Control Conference, Groningen, NL, pp. 693-696, 1993.

The implemented methods rely on an accuracy enhancing square-root technique. 3 The algorithms require less than 30N floating point operations.

None

**Program Text**

* AB09KD EXAMPLE PROGRAM TEXT * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER MMAX, NMAX, NVMAX, NWMAX, PMAX PARAMETER ( MMAX = 20, NMAX = 20, NVMAX = 10, NWMAX = 10, $ PMAX = 20 ) INTEGER LDA, LDAV, LDAW, LDB, LDBV, LDBW, $ LDC, LDCV, LDCW, LDD, LDDV, LDDW PARAMETER ( LDA = NMAX, LDAV = NVMAX, LDAW = NWMAX, $ LDB = NMAX, LDBV = NVMAX, LDBW = NWMAX, $ LDC = PMAX, LDCV = PMAX, LDCW = MMAX, $ LDD = PMAX, LDDV = PMAX, LDDW = MMAX ) INTEGER LIWORK PARAMETER ( LIWORK = 2*MAX( MMAX, PMAX ) ) INTEGER LDW1, LDW2, LDW3, LDW4, LDWORK PARAMETER ( LDW1 = MAX( NVMAX*( NVMAX + 5 ), NVMAX*NMAX + $ MAX( 2*NVMAX, PMAX*NMAX, PMAX*MMAX ) )) PARAMETER ( LDW2 = MAX( NWMAX*( NWMAX + 5 ), NWMAX*NMAX + $ MAX( 2*NWMAX, MMAX*NMAX, PMAX*MMAX ) )) PARAMETER ( LDW3 = NMAX*( 2*NMAX + MAX( NMAX, MMAX, PMAX ) $ + 5 ) + ( NMAX*( NMAX + 1 ) )/2 ) PARAMETER ( LDW4 = NMAX*( MMAX + PMAX + 2 ) + 2*MMAX*PMAX + $ MIN( NMAX, MMAX ) + $ MAX( 3*MMAX + 1, $ MIN( NMAX, MMAX ) + PMAX ) ) PARAMETER ( LDWORK = MAX( LDW1, LDW2, LDW3, LDW4 ) ) * .. Local Scalars .. LOGICAL LEFTW, RIGHTW DOUBLE PRECISION ALPHA, TOL1, TOL2 INTEGER I, INFO, IWARN, J, M, N, NR, NS, NV, NW, P CHARACTER*1 DICO, EQUIL, JOB, ORDSEL, WEIGHT * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), AV(LDAV,NVMAX), AW(LDAW,NWMAX), $ B(LDB,MMAX), BV(LDBV,PMAX), BW(LDBW,MMAX), $ C(LDC,NMAX), CV(LDCV,NVMAX), CW(LDCW,NWMAX), $ D(LDD,MMAX), DV(LDDV,PMAX), DW(LDDW,MMAX), $ DWORK(LDWORK), HSV(NMAX) INTEGER IWORK(LIWORK) * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL AB09KD * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. Executable Statements .. * WRITE ( NOUT, FMT = 99999 ) * Skip the heading in the data file and read the data. READ ( NIN, FMT = '()' ) READ ( NIN, FMT = * ) N, M, P, NV, NW, NR, ALPHA, TOL1, TOL2, $ JOB, DICO, WEIGHT, EQUIL, ORDSEL LEFTW = LSAME( WEIGHT, 'L' ) .OR. LSAME( WEIGHT, 'B' ) RIGHTW = LSAME( WEIGHT, 'R' ) .OR. LSAME( WEIGHT, 'B' ) IF( N.LE.0 .OR. N.GT.NMAX ) THEN WRITE ( NOUT, FMT = 99990 ) N ELSE READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N ) IF( M.LE.0 .OR. M.GT.MMAX ) THEN WRITE ( NOUT, FMT = 99989 ) M ELSE READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1, N ) IF( P.GT.PMAX ) THEN WRITE ( NOUT, FMT = 99988 ) P ELSE READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P ) READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,P ) IF( LEFTW .OR. NV.GT.0 ) THEN IF( NV.LT.0 .OR. NV.GT.NVMAX ) THEN WRITE ( NOUT, FMT = 99986 ) NV ELSE IF( NV.GT.0 ) THEN READ ( NIN, FMT = * ) $ ( ( AV(I,J), J = 1,NV ), I = 1,NV ) READ ( NIN, FMT = * ) $ ( ( BV(I,J), J = 1,P ), I = 1, NV ) READ ( NIN, FMT = * ) $ ( ( CV(I,J), J = 1,NV ), I = 1,P ) END IF IF( LEFTW ) READ ( NIN, FMT = * ) $ ( ( DV(I,J), J = 1,P ), I = 1,P ) END IF END IF IF( RIGHTW ) THEN IF( NW.LT.0 .OR. NW.GT.NWMAX ) THEN WRITE ( NOUT, FMT = 99985 ) NW ELSE IF( NW.GT.0 ) THEN READ ( NIN, FMT = * ) $ ( ( AW(I,J), J = 1,NW ), I = 1,NW ) READ ( NIN, FMT = * ) $ ( ( BW(I,J), J = 1,M ), I = 1, NW ) READ ( NIN, FMT = * ) $ ( ( CW(I,J), J = 1,NW ), I = 1,M ) END IF READ ( NIN, FMT = * ) $ ( ( DW(I,J), J = 1,M ), I = 1,M ) END IF END IF * Find a reduced ssr for (A,B,C,D). CALL AB09KD( JOB, DICO, WEIGHT, EQUIL, ORDSEL, N, NV, NW, $ M, P, NR, ALPHA, A, LDA, B, LDB, C, LDC, $ D, LDD, AV, LDAV, BV, LDBV, CV, LDCV, $ DV, LDDV, AW, LDAW, BW, LDBW, CW, LDCW, $ DW, LDDW, NS, HSV, TOL1, TOL2, IWORK, $ DWORK, LDWORK, IWARN, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE IF( IWARN.NE.0) WRITE ( NOUT, FMT = 99984 ) IWARN WRITE ( NOUT, FMT = 99997 ) NR WRITE ( NOUT, FMT = 99987 ) WRITE ( NOUT, FMT = 99995 ) ( HSV(J), J = 1, NS ) IF( NR.GT.0 ) WRITE ( NOUT, FMT = 99996 ) DO 20 I = 1, NR WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NR ) 20 CONTINUE IF( NR.GT.0 ) WRITE ( NOUT, FMT = 99993 ) DO 40 I = 1, NR WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M ) 40 CONTINUE IF( NR.GT.0 ) WRITE ( NOUT, FMT = 99992 ) DO 60 I = 1, P WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,NR ) 60 CONTINUE WRITE ( NOUT, FMT = 99991 ) DO 70 I = 1, P WRITE ( NOUT, FMT = 99995 ) ( D(I,J), J = 1,M ) 70 CONTINUE END IF END IF END IF END IF STOP * 99999 FORMAT (' AB09KD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from AB09KD = ',I2) 99997 FORMAT (/' The order of reduced model = ',I2) 99996 FORMAT (/' The reduced state dynamics matrix Ar is ') 99995 FORMAT (20(1X,F8.4)) 99993 FORMAT (/' The reduced input/state matrix Br is ') 99992 FORMAT (/' The reduced state/output matrix Cr is ') 99991 FORMAT (/' The reduced input/output matrix Dr is ') 99990 FORMAT (/' N is out of range.',/' N = ',I5) 99989 FORMAT (/' M is out of range.',/' M = ',I5) 99988 FORMAT (/' P is out of range.',/' P = ',I5) 99987 FORMAT (/' The Hankel singular values of weighted ALPHA-stable', $ ' part are') 99986 FORMAT (/' NV is out of range.',/' NV = ',I5) 99985 FORMAT (/' NW is out of range.',/' NW = ',I5) 99984 FORMAT (' IWARN on exit from AB09KD = ',I2) END

AB09KD EXAMPLE PROGRAM DATA (Continuous system) 6 1 1 2 0 0 0.0 1.E-1 1.E-14 N C L S A -3.8637 -7.4641 -9.1416 -7.4641 -3.8637 -1.0000 1.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 1.0000 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0.2000 -1.0000 1.0000 0 1 0 -1.8000 0 1

AB09KD EXAMPLE PROGRAM RESULTS The order of reduced model = 4 The Hankel singular values of weighted ALPHA-stable part are 2.6790 2.1589 0.8424 0.1929 0.0219 0.0011 The reduced state dynamics matrix Ar is -0.2391 0.3072 1.1630 1.1967 -2.9709 -0.2391 2.6270 3.1027 0.0000 0.0000 -0.5137 -1.2842 0.0000 0.0000 0.1519 -0.5137 The reduced input/state matrix Br is -1.0497 -3.7052 0.8223 0.7435 The reduced state/output matrix Cr is -0.4466 0.0143 -0.4780 -0.2013 The reduced input/output matrix Dr is 0.0219