**Purpose**

To construct a state-space representation (A,BS,CS,DS) of the projection of G*W or G*conj(W) containing the poles of G, from the state-space representations (A,B,C,D) and (AW-lambda*EW,BW,CW,DW), of the transfer-function matrices G and W, respectively. G is assumed to be a stable transfer-function matrix and the state matrix A must be in a real Schur form. When computing the stable projection of G*W, it is assumed that G and W have completely distinct poles. When computing the stable projection of G*conj(W), it is assumed that G and conj(W) have completely distinct poles. Note: 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 AB09JW( JOB, DICO, JOBEW, STBCHK, N, M, P, NW, MW, $ A, LDA, B, LDB, C, LDC, D, LDD, AW, LDAW, $ EW, LDEW, BW, LDBW, CW, LDCW, DW, LDDW, IWORK, $ DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER DICO, JOB, JOBEW, STBCHK INTEGER INFO, LDA, LDAW, LDB, LDBW, LDC, LDCW, $ LDD, LDDW, LDEW, LDWORK, M, MW, N, NW, P C .. Array Arguments .. INTEGER IWORK(*) DOUBLE PRECISION A(LDA,*), AW(LDAW,*), B(LDB,*), BW(LDBW,*), $ C(LDC,*), CW(LDCW,*), D(LDD,*), DW(LDDW,*), $ DWORK(*), EW(LDEW,*)

**Mode Parameters**

JOB CHARACTER*1 Specifies the projection to be computed as follows: = 'W': compute the projection of G*W containing the poles of G; = 'C': compute the projection of G*conj(W) containing the poles of G. DICO CHARACTER*1 Specifies the type of the systems as follows: = 'C': G and W are continuous-time systems; = 'D': G and W are discrete-time systems. JOBEW CHARACTER*1 Specifies whether EW is a general square or an identity matrix as follows: = 'G': EW is a general square matrix; = 'I': EW is the identity matrix. STBCHK CHARACTER*1 Specifies whether stability/antistability of W is to be checked as follows: = 'C': check stability if JOB = 'C' or antistability if JOB = 'W'; = 'N': do not check stability or antistability.

N (input) INTEGER The dimension of the state vector of the system with the transfer-function matrix G. N >= 0. M (input) INTEGER The dimension of the input vector of the system with the transfer-function matrix G, and also the dimension of the output vector if JOB = 'W', or of the input vector if JOB = 'C', of the system with the transfer-function matrix W. M >= 0. P (input) INTEGER The dimension of the output vector of the system with the transfer-function matrix G. P >= 0. NW (input) INTEGER The dimension of the state vector of the system with the transfer-function matrix W. NW >= 0. MW (input) INTEGER The dimension of the input vector, if JOB = 'W', or of the output vector, if JOB = 'C', of the system with the transfer-function matrix W. MW >= 0. A (input) DOUBLE PRECISION array, dimension (LDA,N) The leading N-by-N part of this array must contain the state matrix A of the system with the transfer-function matrix G in a real Schur form. LDA INTEGER The leading dimension of the array A. LDA >= MAX(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB,MAX(M,MW)) On entry, the leading N-by-M part of this array must contain the input matrix B of the system with the transfer-function matrix G. On exit, if INFO = 0, the leading N-by-MW part of this array contains the input matrix BS of the projection of G*W, if JOB = 'W', or of G*conj(W), if JOB = 'C'. LDB INTEGER The leading dimension of the array B. LDB >= MAX(1,N). C (input) DOUBLE PRECISION array, dimension (LDC,N) The leading P-by-N part of this array must contain the output/state matrix C of the system with the transfer-function matrix G. The matrix CS is equal to C. LDC INTEGER The leading dimension of the array C. LDC >= MAX(1,P). D (input/output) DOUBLE PRECISION array, dimension (LDB,MAX(M,MW)) On entry, the leading P-by-M part of this array must contain the feedthrough matrix D of the system with the transfer-function matrix G. On exit, if INFO = 0, the leading P-by-MW part of this array contains the feedthrough matrix DS of the projection of G*W, if JOB = 'W', or of G*conj(W), if JOB = 'C'. LDD INTEGER The leading dimension of the array D. LDD >= MAX(1,P). AW (input/output) DOUBLE PRECISION array, dimension (LDAW,NW) On entry, the leading NW-by-NW part of this array must contain the state matrix AW of the system with the transfer-function matrix W. On exit, if INFO = 0, the leading NW-by-NW part of this array contains a condensed matrix as follows: if JOBEW = 'I', it contains the real Schur form of AW; if JOBEW = 'G' and JOB = 'W', it contains a quasi-upper triangular matrix representing the real Schur matrix in the real generalized Schur form of the pair (AW,EW); if JOBEW = 'G', JOB = 'C' and DICO = 'C', it contains a quasi-upper triangular matrix corresponding to the generalized real Schur form of the pair (AW',EW'); if JOBEW = 'G', JOB = 'C' and DICO = 'D', it contains an upper triangular matrix corresponding to the generalized real Schur form of the pair (EW',AW'). LDAW INTEGER The leading dimension of the array AW. LDAW >= MAX(1,NW). EW (input/output) DOUBLE PRECISION array, dimension (LDEW,NW) On entry, if JOBEW = 'G', the leading NW-by-NW part of this array must contain the descriptor matrix EW of the system with the transfer-function matrix W. If JOBEW = 'I', EW is assumed to be an identity matrix and is not referenced. On exit, if INFO = 0 and JOBEW = 'G', the leading NW-by-NW part of this array contains a condensed matrix as follows: if JOB = 'W', it contains an upper triangular matrix corresponding to the real generalized Schur form of the pair (AW,EW); if JOB = 'C' and DICO = 'C', it contains an upper triangular matrix corresponding to the generalized real Schur form of the pair (AW',EW'); if JOB = 'C' and DICO = 'D', it contains a quasi-upper triangular matrix corresponding to the generalized real Schur form of the pair (EW',AW'). LDEW INTEGER The leading dimension of the array EW. LDEW >= MAX(1,NW), if JOBEW = 'G'; LDEW >= 1, if JOBEW = 'I'. BW (input/output) DOUBLE PRECISION array, dimension (LDBW,MBW), where MBW = MW, if JOB = 'W', and MBW = M, if JOB = 'C'. On entry, the leading NW-by-MBW part of this array must contain the input matrix BW of the system with the transfer-function matrix W. On exit, if INFO = 0, the leading NW-by-MBW part of this array contains Q'*BW, where Q is the orthogonal matrix that reduces AW to the real Schur form or the left orthogonal matrix used to reduce the pair (AW,EW), (AW',EW') or (EW',AW') to the generalized real Schur form. LDBW INTEGER The leading dimension of the array BW. LDBW >= MAX(1,NW). CW (input/output) DOUBLE PRECISION array, dimension (LDCW,NW) On entry, the leading PCW-by-NW part of this array must contain the output matrix CW of the system with the transfer-function matrix W, where PCW = M if JOB = 'W' or PCW = MW if JOB = 'C'. On exit, if INFO = 0, the leading PCW-by-NW part of this array contains CW*Q, where Q is the orthogonal matrix that reduces AW to the real Schur form, or CW*Z, where Z is the right orthogonal matrix used to reduce the pair (AW,EW), (AW',EW') or (EW',AW') to the generalized real Schur form. LDCW INTEGER The leading dimension of the array CW. LDCW >= MAX(1,PCW), where PCW = M if JOB = 'W', or PCW = MW if JOB = 'C'. DW (input) DOUBLE PRECISION array, dimension (LDDW,MBW), where MBW = MW if JOB = 'W', and MBW = M if JOB = 'C'. The leading PCW-by-MBW part of this array must contain the feedthrough matrix DW of the system with the transfer-function matrix W, where PCW = M if JOB = 'W', or PCW = MW if JOB = 'C'. LDDW INTEGER LDDW >= MAX(1,PCW), where PCW = M if JOB = 'W', or PCW = MW if JOB = 'C'.

IWORK INTEGER array, dimension (LIWORK) LIWORK = 0, if JOBEW = 'I'; LIWORK = NW+N+6, if JOBEW = 'G'. 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 >= LW1, if JOBEW = 'I', LDWORK >= LW2, if JOBEW = 'G', where LW1 = MAX( 1, NW*(NW+5), NW*N + MAX( a, N*MW, P*MW ) ) a = 0, if DICO = 'C' or JOB = 'W', a = 2*NW, if DICO = 'D' and JOB = 'C'; LW2 = MAX( 2*NW*NW + MAX( 11*NW+16, NW*M, MW*NW ), NW*N + MAX( NW*N+N*N, MW*N, P*MW ) ). For good performance, LDWORK should be larger.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 1: the reduction of the pair (AW,EW) to the real generalized Schur form failed (JOBEW = 'G'), or the reduction of the matrix AW to the real Schur form failed (JOBEW = 'I); = 2: the solution of the Sylvester equation failed because the matrix A and the pencil AW-lambda*EW have common eigenvalues (if JOB = 'W'), or the pencil -AW-lambda*EW and A have common eigenvalues (if JOB = 'C' and DICO = 'C'), or the pencil AW-lambda*EW has an eigenvalue which is the reciprocal of one of eigenvalues of A (if JOB = 'C' and DICO = 'D'); = 3: the solution of the Sylvester equation failed because the matrices A and AW have common eigenvalues (if JOB = 'W'), or the matrices A and -AW have common eigenvalues (if JOB = 'C' and DICO = 'C'), or the matrix A has an eigenvalue which is the reciprocal of one of eigenvalues of AW (if JOB = 'C' and DICO = 'D'); = 4: JOB = 'W' and the pair (AW,EW) has not completely unstable generalized eigenvalues, or JOB = 'C' and the pair (AW,EW) has not completely stable generalized eigenvalues.

If JOB = 'W', the matrices of the stable projection of G*W are computed as BS = B*DW + Y*BW, CS = C, DS = D*DW, where Y satisfies the generalized Sylvester equation -A*Y*EW + Y*AW + B*CW = 0. If JOB = 'C', the matrices of the stable projection of G*conj(W) are computed using the following formulas: - for a continuous-time system, the matrices BS, CS and DS of the stable projection are computed as BS = B*DW' + Y*CW', CS = C, DS = D*DW', where Y satisfies the generalized Sylvester equation A*Y*EW' + Y*AW' + B*BW' = 0. - for a discrete-time system, the matrices BS, CS and DS of the stable projection are computed as BS = B*DW' + A*Y*CW', CS = C, DS = D*DW' + C*Y*CW', where Y satisfies the generalized Sylvester equation Y*EW' - A*Y*AW' = B*BW'.

[1] Varga, A. Efficient and numerically reliable implementation of the frequency-weighted Hankel-norm approximation model reduction approach. Proc. 2001 ECC, Porto, Portugal, 2001. [2] Zhou, K. Frequency-weighted H-infinity norm and optimal Hankel norm model reduction. IEEE Trans. Autom. Control, vol. 40, pp. 1687-1699, 1995.

The implemented methods rely on numerically stable algorithms.

None

**Program Text**

None

None

None