Purpose
To perform the D-step in the D-K iteration. It handles continuous-time case.Specification
SUBROUTINE SB10MD( NC, MP, LENDAT, F, ORD, MNB, NBLOCK, ITYPE,
$ QUTOL, A, LDA, B, LDB, C, LDC, D, LDD, OMEGA,
$ TOTORD, AD, LDAD, BD, LDBD, CD, LDCD, DD, LDDD,
$ MJU, IWORK, LIWORK, DWORK, LDWORK, ZWORK,
$ LZWORK, INFO )
C .. Scalar Arguments ..
INTEGER F, INFO, LDA, LDAD, LDB, LDBD, LDC, LDCD, LDD,
$ LDDD, LDWORK, LENDAT, LIWORK, LZWORK, MNB, MP,
$ NC, ORD, TOTORD
DOUBLE PRECISION QUTOL
C .. Array Arguments ..
INTEGER ITYPE(*), IWORK(*), NBLOCK(*)
DOUBLE PRECISION A(LDA, *), AD(LDAD, *), B(LDB, *), BD(LDBD, *),
$ C(LDC, *), CD(LDCD, *), D(LDD, *), DD(LDDD, *),
$ DWORK(*), MJU(*), OMEGA(*)
COMPLEX*16 ZWORK(*)
Arguments
Input/Output Parameters
NC (input) INTEGER
The order of the matrix A. NC >= 0.
MP (input) INTEGER
The order of the matrix D. MP >= 0.
LENDAT (input) INTEGER
The length of the vector OMEGA. LENDAT >= 2.
F (input) INTEGER
The number of the measurements and controls, i.e.,
the size of the block I_f in the D-scaling system.
F >= 0.
ORD (input/output) INTEGER
The MAX order of EACH block in the fitting procedure.
ORD <= LENDAT-1.
On exit, if ORD < 1 then ORD = 1.
MNB (input) INTEGER
The number of diagonal blocks in the block structure of
the uncertainty, and the length of the vectors NBLOCK
and ITYPE. 1 <= MNB <= MP.
NBLOCK (input) INTEGER array, dimension (MNB)
The vector of length MNB containing the block structure
of the uncertainty. NBLOCK(I), I = 1:MNB, is the size of
each block.
ITYPE (input) INTEGER array, dimension (MNB)
The vector of length MNB indicating the type of each
block.
For I = 1 : MNB,
ITYPE(I) = 1 indicates that the corresponding block is a
real block. IN THIS CASE ONLY MJU(JW) WILL BE ESTIMATED
CORRECTLY, BUT NOT D(S)!
ITYPE(I) = 2 indicates that the corresponding block is a
complex block. THIS IS THE ONLY ALLOWED VALUE NOW!
NBLOCK(I) must be equal to 1 if ITYPE(I) is equal to 1.
QUTOL (input) DOUBLE PRECISION
The acceptable mean relative error between the D(jw) and
the frequency responce of the estimated block
[ADi,BDi;CDi,DDi]. When it is reached, the result is
taken as good enough.
A good value is QUTOL = 2.0.
If QUTOL < 0 then only mju(jw) is being estimated,
not D(s).
A (input/output) DOUBLE PRECISION array, dimension (LDA,NC)
On entry, the leading NC-by-NC part of this array must
contain the A matrix of the closed-loop system.
On exit, if MP > 0, the leading NC-by-NC part of this
array contains an upper Hessenberg matrix similar to A.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,NC).
B (input/output) DOUBLE PRECISION array, dimension (LDB,MP)
On entry, the leading NC-by-MP part of this array must
contain the B matrix of the closed-loop system.
On exit, the leading NC-by-MP part of this array contains
the transformed B matrix corresponding to the Hessenberg
form of A.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1,NC).
C (input/output) DOUBLE PRECISION array, dimension (LDC,NC)
On entry, the leading MP-by-NC part of this array must
contain the C matrix of the closed-loop system.
On exit, the leading MP-by-NC part of this array contains
the transformed C matrix corresponding to the Hessenberg
form of A.
LDC INTEGER
The leading dimension of the array C. LDC >= MAX(1,MP).
D (input) DOUBLE PRECISION array, dimension (LDD,MP)
The leading MP-by-MP part of this array must contain the
D matrix of the closed-loop system.
LDD INTEGER
The leading dimension of the array D. LDD >= MAX(1,MP).
OMEGA (input) DOUBLE PRECISION array, dimension (LENDAT)
The vector with the frequencies.
TOTORD (output) INTEGER
The TOTAL order of the D-scaling system.
TOTORD is set to zero, if QUTOL < 0.
AD (output) DOUBLE PRECISION array, dimension (LDAD,MP*ORD)
The leading TOTORD-by-TOTORD part of this array contains
the A matrix of the D-scaling system.
Not referenced if QUTOL < 0.
LDAD INTEGER
The leading dimension of the array AD.
LDAD >= MAX(1,MP*ORD), if QUTOL >= 0;
LDAD >= 1, if QUTOL < 0.
BD (output) DOUBLE PRECISION array, dimension (LDBD,MP+F)
The leading TOTORD-by-(MP+F) part of this array contains
the B matrix of the D-scaling system.
Not referenced if QUTOL < 0.
LDBD INTEGER
The leading dimension of the array BD.
LDBD >= MAX(1,MP*ORD), if QUTOL >= 0;
LDBD >= 1, if QUTOL < 0.
CD (output) DOUBLE PRECISION array, dimension (LDCD,MP*ORD)
The leading (MP+F)-by-TOTORD part of this array contains
the C matrix of the D-scaling system.
Not referenced if QUTOL < 0.
LDCD INTEGER
The leading dimension of the array CD.
LDCD >= MAX(1,MP+F), if QUTOL >= 0;
LDCD >= 1, if QUTOL < 0.
DD (output) DOUBLE PRECISION array, dimension (LDDD,MP+F)
The leading (MP+F)-by-(MP+F) part of this array contains
the D matrix of the D-scaling system.
Not referenced if QUTOL < 0.
LDDD INTEGER
The leading dimension of the array DD.
LDDD >= MAX(1,MP+F), if QUTOL >= 0;
LDDD >= 1, if QUTOL < 0.
MJU (output) DOUBLE PRECISION array, dimension (LENDAT)
The vector with the upper bound of the structured
singular value (mju) for each frequency in OMEGA.
Workspace
IWORK INTEGER array, dimension (LIWORK)
LIWORK INTEGER
The length of the array IWORK.
LIWORK >= MAX( NC, 4*MNB-2, MP, 2*ORD+1 ), if QUTOL >= 0;
LIWORK >= MAX( NC, 4*MNB-2, MP ), if QUTOL < 0.
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK, DWORK(2) returns the optimal value of LZWORK,
and DWORK(3) returns an estimate of the minimum reciprocal
of the condition numbers (with respect to inversion) of
the generated Hessenberg matrices.
LDWORK INTEGER
The length of the array DWORK.
LDWORK >= MAX( 3, LWM, LWD ), where
LWM = LWA + MAX( NC + MAX( NC, MP-1 ),
2*MP*MP*MNB - MP*MP + 9*MNB*MNB +
MP*MNB + 11*MP + 33*MNB - 11 );
LWD = LWB + MAX( 2, LW1, LW2, LW3, LW4, 2*ORD ),
if QUTOL >= 0;
LWD = 0, if QUTOL < 0;
LWA = MP*LENDAT + 2*MNB + MP - 1;
LWB = LENDAT*(MP + 2) + ORD*(ORD + 2) + 1;
LW1 = 2*LENDAT + 4*HNPTS; HNPTS = 2048;
LW2 = LENDAT + 6*HNPTS; MN = MIN( 2*LENDAT, 2*ORD+1 );
LW3 = 2*LENDAT*(2*ORD + 1) + MAX( 2*LENDAT, 2*ORD + 1 ) +
MAX( MN + 6*ORD + 4, 2*MN + 1 );
LW4 = MAX( ORD*ORD + 5*ORD, 6*ORD + 1 + MIN( 1, ORD ) ).
ZWORK COMPLEX*16 array, dimension (LZWORK)
LZWORK INTEGER
The length of the array ZWORK.
LZWORK >= MAX( LZM, LZD ), where
LZM = MAX( MP*MP + NC*MP + NC*NC + 2*NC,
6*MP*MP*MNB + 13*MP*MP + 6*MNB + 6*MP - 3 );
LZD = MAX( LENDAT*(2*ORD + 3), ORD*ORD + 3*ORD + 1 ),
if QUTOL >= 0;
LZD = 0, if QUTOL < 0.
Error Indicator
INFO (output) INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: if one or more values w in OMEGA are (close to
some) poles of the closed-loop system, i.e., the
matrix jw*I - A is (numerically) singular;
= 2: the block sizes must be positive integers;
= 3: the sum of block sizes must be equal to MP;
= 4: the size of a real block must be equal to 1;
= 5: the block type must be either 1 or 2;
= 6: errors in solving linear equations or in matrix
inversion;
= 7: errors in computing eigenvalues or singular values.
= 1i: INFO on exit from SB10YD is i. (1i means 10 + i.)
Method
I. First, W(jw) for the given closed-loop system is being
estimated.
II. Now, AB13MD SLICOT subroutine can obtain the D(jw) scaling
system with respect to NBLOCK and ITYPE, and colaterally,
mju(jw).
If QUTOL < 0 then the estimations stop and the routine exits.
III. Now that we have D(jw), SB10YD subroutine can do block-by-
block fit. For each block it tries with an increasing order
of the fit, starting with 1 until the
(mean quadratic error + max quadratic error)/2
between the Dii(jw) and the estimated frequency responce
of the block becomes less than or equal to the routine
argument QUTOL, or the order becomes equal to ORD.
IV. Arrange the obtained blocks in the AD, BD, CD and DD
matrices and estimate the total order of D(s), TOTORD.
V. Add the system I_f to the system obtained in IV.
References
[1] Balas, G., Doyle, J., Glover, K., Packard, A. and Smith, R.
Mu-analysis and Synthesis toolbox - User's Guide,
The Mathworks Inc., Natick, MA, USA, 1998.
Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None