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 square-root or balancing-free square-root
Balance & Truncate (B&T) or Singular Perturbation Approximation
(SPA) model reduction methods. The algorithm tries to minimize
the norm of the frequency-weighted error
||V*(G-Gr)*W||
where G and Gr are the transfer-function matrices of the original
and reduced order models, respectively, and V and W are
frequency-weighting transfer-function matrices. V and W must not
have poles on the imaginary axis for a continuous-time
system or on the unit circle for a discrete-time system.
If G is unstable, only the ALPHA-stable part of G is reduced.
In case of possible pole-zero cancellations in V*G and/or G*W,
the absolute values of parameters ALPHAO and/or ALPHAC must be
different from 1.
Specification
SUBROUTINE AB09ID( DICO, JOBC, JOBO, JOB, WEIGHT, EQUIL, ORDSEL,
$ N, M, P, NV, PV, NW, MW, NR, ALPHA, ALPHAC,
$ ALPHAO, 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, JOBC, JOBO, ORDSEL, WEIGHT
INTEGER INFO, IWARN, LDA, LDAV, LDAW, LDB, LDBV, LDBW,
$ LDC, LDCV, LDCW, LDD, LDDV, LDDW, LDWORK, M, MW,
$ N, NR, NS, NV, NW, P, PV
DOUBLE PRECISION ALPHA, ALPHAC, ALPHAO, 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(*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of the original system as follows:
= 'C': continuous-time system;
= 'D': discrete-time system.
JOBC CHARACTER*1
Specifies the choice of frequency-weighted controllability
Grammian as follows:
= 'S': choice corresponding to a combination method [4]
of the approaches of Enns [1] and Lin-Chiu [2,3];
= 'E': choice corresponding to the stability enhanced
modified combination method of [4].
JOBO CHARACTER*1
Specifies the choice of frequency-weighted observability
Grammian as follows:
= 'S': choice corresponding to a combination method [4]
of the approaches of Enns [1] and Lin-Chiu [2,3];
= 'E': choice corresponding to the stability enhanced
modified combination method of [4].
JOB CHARACTER*1
Specifies the model reduction approach to be used
as follows:
= 'B': use the square-root Balance & Truncate method;
= 'F': use the balancing-free square-root
Balance & Truncate method;
= 'S': use the square-root Singular Perturbation
Approximation method;
= 'P': use the balancing-free square-root
Singular Perturbation Approximation method.
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.
Input/Output Parameters
N (input) INTEGER
The order of the original state-space representation,
i.e., the order of the matrix A. N >= 0.
M (input) INTEGER
The number of system inputs. M >= 0.
P (input) INTEGER
The number of system outputs. P >= 0.
NV (input) INTEGER
The order of the matrix AV. Also the number of rows of
the matrix BV and the number of columns of the matrix CV.
NV represents the dimension of the state vector of the
system with the transfer-function matrix V. NV >= 0.
PV (input) INTEGER
The number of rows of the matrices CV and DV. PV >= 0.
PV represents the dimension of the output vector of the
system with the transfer-function matrix V.
NW (input) INTEGER
The order of the matrix AW. Also the number of rows of
the matrix BW and the number of columns of the matrix CW.
NW represents the dimension of the state vector of the
system with the transfer-function matrix W. NW >= 0.
MW (input) INTEGER
The number of columns of the matrices BW and DW. MW >= 0.
MW represents the dimension of the input vector of the
system with the transfer-function matrix W.
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),NMIN), where NR is the desired order
on entry, NMIN is the number of frequency-weighted Hankel
singular values greater than NS*EPS*S1, EPS is the
machine precision (see LAPACK Library Routine DLAMCH)
and S1 is the largest Hankel singular value (computed
in HSV(1)); NR can be further reduced to ensure
HSV(NR-NU) > HSV(NR+1-NU);
if ORDSEL = 'A', NR is the sum of NU and the number of
Hankel singular values greater than MAX(TOL1,NS*EPS*S1).
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.
ALPHAC (input) DOUBLE PRECISION
Combination method parameter for defining the
frequency-weighted controllability Grammian (see METHOD);
ABS(ALPHAC) <= 1.
ALPHAO (input) DOUBLE PRECISION
Combination method parameter for defining the
frequency-weighted observability Grammian (see METHOD);
ABS(ALPHAO) <= 1.
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.
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
the system with the transfer-function matrix V.
On exit, if WEIGHT = 'L' or 'B', MIN(N,M,P) > 0 and
INFO = 0, the leading NVR-by-NVR part of this array
contains the state matrix of a minimal realization of V
in a real Schur form. NVR is returned in IWORK(2).
AV is not referenced if WEIGHT = 'R' or 'N',
or MIN(N,M,P) = 0.
LDAV INTEGER
The leading dimension of 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 the
system with the transfer-function matrix V.
On exit, if WEIGHT = 'L' or 'B', MIN(N,M,P) > 0 and
INFO = 0, the leading NVR-by-P part of this array contains
the input matrix of a minimal realization of V.
BV is not referenced if WEIGHT = 'R' or 'N',
or MIN(N,M,P) = 0.
LDBV INTEGER
The leading dimension of 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 PV-by-NV
part of this array must contain the output matrix CV of
the system with the transfer-function matrix V.
On exit, if WEIGHT = 'L' or 'B', MIN(N,M,P) > 0 and
INFO = 0, the leading PV-by-NVR part of this array
contains the output matrix of a minimal realization of V.
CV is not referenced if WEIGHT = 'R' or 'N',
or MIN(N,M,P) = 0.
LDCV INTEGER
The leading dimension of array CV.
LDCV >= MAX(1,PV), if WEIGHT = 'L' or 'B';
LDCV >= 1, if WEIGHT = 'R' or 'N'.
DV (input) DOUBLE PRECISION array, dimension (LDDV,P)
If WEIGHT = 'L' or 'B', the leading PV-by-P part of this
array must contain the feedthrough matrix DV of the system
with the transfer-function matrix V.
DV is not referenced if WEIGHT = 'R' or 'N',
or MIN(N,M,P) = 0.
LDDV INTEGER
The leading dimension of array DV.
LDDV >= MAX(1,PV), 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
the system with the transfer-function matrix W.
On exit, if WEIGHT = 'R' or 'B', MIN(N,M,P) > 0 and
INFO = 0, the leading NWR-by-NWR part of this array
contains the state matrix of a minimal realization of W
in a real Schur form. NWR is returned in IWORK(3).
AW is not referenced if WEIGHT = 'L' or 'N',
or MIN(N,M,P) = 0.
LDAW INTEGER
The leading dimension of 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,MW)
On entry, if WEIGHT = 'R' or 'B', the leading NW-by-MW
part of this array must contain the input matrix BW of the
system with the transfer-function matrix W.
On exit, if WEIGHT = 'R' or 'B', MIN(N,M,P) > 0 and
INFO = 0, the leading NWR-by-MW part of this array
contains the input matrix of a minimal realization of W.
BW is not referenced if WEIGHT = 'L' or 'N',
or MIN(N,M,P) = 0.
LDBW INTEGER
The leading dimension of 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 the
system with the transfer-function matrix W.
On exit, if WEIGHT = 'R' or 'B', MIN(N,M,P) > 0 and
INFO = 0, the leading M-by-NWR part of this array contains
the output matrix of a minimal realization of W.
CW is not referenced if WEIGHT = 'L' or 'N',
or MIN(N,M,P) = 0.
LDCW INTEGER
The leading dimension of array CW.
LDCW >= MAX(1,M), if WEIGHT = 'R' or 'B';
LDCW >= 1, if WEIGHT = 'L' or 'N'.
DW (input) DOUBLE PRECISION array, dimension (LDDW,MW)
If WEIGHT = 'R' or 'B', the leading M-by-MW part of this
array must contain the feedthrough matrix DW of the system
with the transfer-function matrix W.
DW is not referenced if WEIGHT = 'L' or 'N',
or MIN(N,M,P) = 0.
LDDW INTEGER
The leading dimension of 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 frequency-weighted Hankel singular values, ordered
decreasingly, of the ALPHA-stable part of the original
system.
Tolerances
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*S1, where c is a constant in the
interval [0.00001,0.001], and S1 is the largest
frequency-weighted Hankel singular value of the
ALPHA-stable part of the original system (computed
in HSV(1)).
If TOL1 <= 0 on entry, the used default value is
TOL1 = NS*EPS*S1, 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*S1.
This value is used by default if TOL2 <= 0 on entry.
If TOL2 > 0 and ORDSEL = 'A', then TOL2 <= TOL1.
Workspace
IWORK INTEGER array, dimension
( MAX( 3, LIWRK1, LIWRK2, LIWRK3 ) ), where
LIWRK1 = 0, if JOB = 'B';
LIWRK1 = N, if JOB = 'F';
LIWRK1 = 2*N, if JOB = 'S' or 'P';
LIWRK2 = 0, if WEIGHT = 'R' or 'N' or NV = 0;
LIWRK2 = NV+MAX(P,PV), if WEIGHT = 'L' or 'B' and NV > 0;
LIWRK3 = 0, if WEIGHT = 'L' or 'N' or NW = 0;
LIWRK3 = NW+MAX(M,MW), if WEIGHT = 'R' or 'B' and NW > 0.
On exit, if INFO = 0, IWORK(1) contains the order of a
minimal realization of the stable part of the system,
IWORK(2) and IWORK(3) contain the actual orders
of the state space realizations of V and W, respectively.
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( LMINL, LMINR, LRCF,
2*N*N + MAX( 1, LLEFT, LRIGHT, 2*N*N+5*N,
N*MAX(M,P) ) ),
where
LMINL = 0, if WEIGHT = 'R' or 'N' or NV = 0; otherwise,
LMINL = MAX(LLCF,NV+MAX(NV,3*P)) if P = PV;
LMINL = MAX(P,PV)*(2*NV+MAX(P,PV))+
MAX(LLCF,NV+MAX(NV,3*P,3*PV)) if P <> PV;
LRCF = 0, and
LMINR = 0, if WEIGHT = 'L' or 'N' or NW = 0; otherwise,
LMINR = NW+MAX(NW,3*M) if M = MW;
LMINR = 2*NW*MAX(M,MW)+NW+MAX(NW,3*M,3*MW) if M <> MW;
LLCF = PV*(NV+PV)+PV*NV+MAX(NV*(NV+5), PV*(PV+2),
4*PV, 4*P);
LRCF = MW*(NW+MW)+MAX(NW*(NW+5),MW*(MW+2),4*MW,4*M)
LLEFT = (N+NV)*(N+NV+MAX(N+NV,PV)+5)
if WEIGHT = 'L' or 'B' and PV > 0;
LLEFT = N*(P+5) if WEIGHT = 'R' or 'N' or PV = 0;
LRIGHT = (N+NW)*(N+NW+MAX(N+NW,MW)+5)
if WEIGHT = 'R' or 'B' and MW > 0;
LRIGHT = N*(M+5) if WEIGHT = 'L' or 'N' or MW = 0.
For optimum performance LDWORK should be larger.
Warning Indicator
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 corresponds
to repeated singular values for the ALPHA-stable
part, which are neither all included nor all
excluded from the reduced model; in this case, the
resulting NR is automatically decreased to exclude
all repeated singular values;
= 3: 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.
= 10+K: K violations of the numerical stability condition
occured during the assignment of eigenvalues in the
SLICOT Library routines SB08CD and/or SB08DD.
Error Indicator
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 to a real Schur form of the state
matrix of a minimal realization of V failed;
= 4: a failure was detected during the ordering of the
real Schur form of the state matrix of a minimal
realization of V or in the iterative process to
compute a left coprime factorization with inner
denominator;
= 5: if DICO = 'C' and the matrix AV has an observable
eigenvalue on the imaginary axis, or DICO = 'D' and
AV has an observable eigenvalue on the unit circle;
= 6: the reduction to a real Schur form of the state
matrix of a minimal realization of W failed;
= 7: a failure was detected during the ordering of the
real Schur form of the state matrix of a minimal
realization of W or in the iterative process to
compute a right coprime factorization with inner
denominator;
= 8: if DICO = 'C' and the matrix AW has a controllable
eigenvalue on the imaginary axis, or DICO = 'D' and
AW has a controllable eigenvalue on the unit circle;
= 9: the computation of eigenvalues failed;
= 10: the computation of Hankel singular values failed.
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 AB09ID 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 norm of the frequency-weighted error
V*(G-Gr)*W, (3)
where V and W are transfer-function matrices without poles on the
imaginary axis in continuous-time case or on the unit circle in
discrete-time case.
The following procedure is used to reduce G:
1) Decompose additively G, of order N, as
G = G1 + G2,
such that G1 = (A1,B1,C1,D) has only ALPHA-stable poles and
G2 = (A2,B2,C2,0), of order NU, has only ALPHA-unstable poles.
2) Compute for G1 a B&T or SPA frequency-weighted approximation
G1r of order NR-NU using the combination method or the
modified combination method of [4].
3) Assemble the reduced model Gr as
Gr = G1r + G2.
For the frequency-weighted reduction of the ALPHA-stable part,
several methods described in [4] can be employed in conjunction
with the combination method and modified combination method
proposed in [4].
If JOB = 'B', the square-root B&T method is used.
If JOB = 'F', the balancing-free square-root version of the
B&T method is used.
If JOB = 'S', the square-root version of the SPA method is used.
If JOB = 'P', the balancing-free square-root version of the
SPA method is used.
For each of these methods, left and right truncation matrices
are determined using the Cholesky factors of an input
frequency-weighted controllability Grammian P and an output
frequency-weighted observability Grammian Q.
P and Q are computed from the controllability Grammian Pi of G*W
and the observability Grammian Qo of V*G. Using special
realizations of G*W and V*G, Pi and Qo are computed in the
partitioned forms
Pi = ( P11 P12 ) and Qo = ( Q11 Q12 ) ,
( P12' P22 ) ( Q12' Q22 )
where P11 and Q11 are the leading N-by-N parts of Pi and Qo,
respectively. Let P0 and Q0 be non-negative definite matrices
defined below
-1
P0 = P11 - ALPHAC**2*P12*P22 *P21 ,
-1
Q0 = Q11 - ALPHAO**2*Q12*Q22 *Q21.
The frequency-weighted controllability and observability
Grammians, P and Q, respectively, are defined as follows:
P = P0 if JOBC = 'S' (standard combination method [4]);
P = P1 >= P0 if JOBC = 'E', where P1 is the controllability
Grammian defined to enforce stability for a modified combination
method of [4];
Q = Q0 if JOBO = 'S' (standard combination method [4]);
Q = Q1 >= Q0 if JOBO = 'E', where Q1 is the observability
Grammian defined to enforce stability for a modified combination
method of [4].
If JOBC = JOBO = 'S' and ALPHAC = ALPHAO = 0, the choice of
Grammians corresponds to the method of Enns [1], while if
ALPHAC = ALPHAO = 1, the choice of Grammians corresponds
to the method of Lin and Chiu [2,3].
If JOBC = 'S' and ALPHAC = 1, no pole-zero cancellations must
occur in G*W. If JOBO = 'S' and ALPHAO = 1, no pole-zero
cancellations must occur in V*G. The presence of pole-zero
cancellations leads to meaningless results and must be avoided.
The frequency-weighted Hankel singular values HSV(1), ....,
HSV(N) are computed as the square roots of the eigenvalues
of the product P*Q.
References
[1] Enns, D.
Model reduction with balanced realizations: An error bound
and a frequency weighted generalization.
Proc. 23-th CDC, Las Vegas, pp. 127-132, 1984.
[2] Lin, C.-A. and Chiu, T.-Y.
Model reduction via frequency-weighted balanced realization.
Control Theory and Advanced Technology, vol. 8,
pp. 341-351, 1992.
[3] Sreeram, V., Anderson, B.D.O and Madievski, A.G.
New results on frequency weighted balanced reduction
technique.
Proc. ACC, Seattle, Washington, pp. 4004-4009, 1995.
[4] Varga, A. and Anderson, B.D.O.
Square-root balancing-free methods for the frequency-weighted
balancing related model reduction.
(report in preparation)
Numerical Aspects
The implemented methods rely on accuracy enhancing square-root techniques.Further Comments
NoneExample
Program Text
* AB09ID EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER MMAX, MWMAX, NMAX, NVMAX, NWMAX, PMAX, PVMAX
PARAMETER ( MMAX = 20, MWMAX = 20,
$ NMAX = 20, NVMAX = 20, NWMAX = 20,
$ PMAX = 20, PVMAX = 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 = PVMAX, LDCW = MMAX,
$ LDD = PMAX, LDDV = PVMAX, LDDW = MMAX )
INTEGER LIWORK
PARAMETER ( LIWORK = MAX( 2*NMAX,
$ NVMAX + MAX( PMAX, PVMAX ),
$ NWMAX + MAX( MMAX, MWMAX ) ) )
INTEGER LDW1, LDW2, LDW3, LDW4, LDW5, LDW6, LDW7, LDW8,
$ LDWORK
PARAMETER ( LDW1 = NMAX + NVMAX, LDW2 = NMAX + NWMAX,
$ LDW3 = MAX( LDW1*( LDW1 + MAX( LDW1, PVMAX ) +
$ 5 ), NMAX*( PMAX + 5 ) ),
$ LDW4 = MAX( LDW2*( LDW2 + MAX( LDW2, MWMAX ) +
$ 5 ), NMAX*( MMAX + 5 ) ),
$ LDW5 = PVMAX*( NVMAX + PVMAX ) + PVMAX*NVMAX +
$ MAX( NVMAX*( NVMAX + 5 ), 4*PVMAX,
$ PVMAX*( PVMAX + 2 ), 4*PMAX ),
$ LDW6 = MAX( PMAX, PVMAX )*( 2*NVMAX +
$ MAX( PMAX, PVMAX ) ) +
$ MAX( LDW5, NVMAX +
$ MAX( NVMAX, 3*PMAX, 3*PVMAX )
$ ),
$ LDW7 = MAX( NWMAX + MAX( NWMAX, 3*MMAX ),
$ 2*NWMAX*MAX( MMAX, MWMAX ) +
$ NWMAX + MAX( NWMAX, 3*MMAX,
$ 3*MWMAX ) ),
$ LDW8 = MWMAX*( NWMAX + MWMAX ) +
$ MAX( NWMAX*( NWMAX + 5 ), 4*MWMAX,
$ MWMAX*( MWMAX + 2 ), 4*MMAX ) )
PARAMETER ( LDWORK = MAX( LDW6, LDW7, LDW8,
$ 2*NMAX*NMAX +
$ MAX( 1, LDW3, LDW4,
$ 2*NMAX*NMAX + 5*NMAX,
$ NMAX*MAX( MMAX, PMAX ) ) )
$ )
* .. Local Scalars ..
LOGICAL LEFTW, RIGHTW
DOUBLE PRECISION ALPHA, ALPHAC, ALPHAO, TOL1, TOL2
INTEGER I, INFO, IWARN, J, M, MW, N, NR, NS, NV, NW, P,
$ PV
CHARACTER*1 DICO, EQUIL, JOB, JOBC, JOBO, ORDSEL, WEIGHT
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), AV(LDAV,NVMAX), AW(LDAW,NWMAX),
$ B(LDB,MMAX), BV(LDBV,PMAX), BW(LDBW,MWMAX),
$ C(LDC,NMAX), CV(LDCV,NVMAX), CW(LDCW,NWMAX),
$ D(LDD,MMAX), DV(LDDV,PMAX), DW(LDDW,MWMAX),
$ DWORK(LDWORK), HSV(NMAX)
INTEGER IWORK(LIWORK)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL AB09ID
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. 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, PV, NW, MW, NR,
$ ALPHA, ALPHAC, ALPHAO, TOL1, TOL2,
$ DICO, JOBC, JOBO, JOB, 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.LE.0 .OR. 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 ) 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 )
IF( PV.LE.0 .OR. PV.GT.PVMAX ) THEN
WRITE ( NOUT, FMT = 99985 ) PV
ELSE
READ ( NIN, FMT = * )
$ ( ( CV(I,J), J = 1,NV ), I = 1,PV )
END IF
END IF
IF( PV.LE.0 .OR. PV.GT.PVMAX ) THEN
WRITE ( NOUT, FMT = 99985 ) PV
ELSE
READ ( NIN, FMT = * )
$ ( ( DV(I,J), J = 1,P ), I = 1,PV )
END IF
END IF
END IF
IF( RIGHTW ) THEN
IF( NW.LT.0 .OR. NW.GT.NWMAX ) THEN
WRITE ( NOUT, FMT = 99984 ) NW
ELSE
IF( NW.GT.0 ) THEN
READ ( NIN, FMT = * )
$ ( ( AW(I,J), J = 1,NW ), I = 1,NW )
IF( MW.LE.0 .OR. MW.GT.MWMAX ) THEN
WRITE ( NOUT, FMT = 99983 ) MW
ELSE
READ ( NIN, FMT = * )
$ ( ( BW(I,J), J = 1,MW ), I = 1,NW )
END IF
READ ( NIN, FMT = * )
$ ( ( CW(I,J), J = 1,NW ), I = 1,M )
END IF
IF( MW.LE.0 .OR. MW.GT.MWMAX ) THEN
WRITE ( NOUT, FMT = 99983 ) MW
ELSE
READ ( NIN, FMT = * )
$ ( ( DW(I,J), J = 1,MW ), I = 1,M )
END IF
END IF
END IF
* Find a reduced ssr for (A,B,C,D).
CALL AB09ID( DICO, JOBC, JOBO, JOB, WEIGHT, EQUIL,
$ ORDSEL, N, M, P, NV, PV, NW, MW, NR, ALPHA,
$ ALPHAC, ALPHAO, 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 = 99982 ) 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 (' AB09ID EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB09ID = ',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 (/' PV is out of range.',/' PV = ',I5)
99984 FORMAT (/' NW is out of range.',/' NW = ',I5)
99983 FORMAT (/' MW is out of range.',/' MW = ',I5)
99982 FORMAT (' IWARN on exit from AB09ID = ',I2)
END
Program Data
AB09ID EXAMPLE PROGRAM DATA (Continuous system)
3 1 1 6 1 0 0 2 0.0 0.0 0.0 0.1E0 0.0 C S S F L S F
-26.4000 6.4023 4.3868
32.0000 0 0
0 8.0000 0
16
0
0
9.2994 1.1624 0.1090
0
-1.0000 0 4.0000 -9.2994 -1.1624 -0.1090
0 2.0000 0 -9.2994 -1.1624 -0.1090
0 0 -3.0000 -9.2994 -1.1624 -0.1090
16.0000 16.0000 16.0000 -26.4000 6.4023 4.3868
0 0 0 32.0000 0 0
0 0 0 0 8.0000 0
1
1
1
0
0
0
1 1 1 0 0 0
0
Program Results
AB09ID EXAMPLE PROGRAM RESULTS The order of reduced model = 2 The Hankel singular values of weighted ALPHA-stable part are 3.8253 0.2005 The reduced state dynamics matrix Ar is 9.1900 0.0000 0.0000 -34.5297 The reduced input/state matrix Br is 11.9593 16.9329 The reduced state/output matrix Cr is 2.8955 6.9152 The reduced input/output matrix Dr is 0.0000