**Purpose**

To perform the D-step in the D-K iteration. It handles continuous-time case.

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(*)

**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.

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.

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.)

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.

[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.

None

**Program Text**

None

None

None