**Purpose**

To find the singular value decomposition (SVD) giving the system order, using the triangular factor of the concatenated block Hankel matrices. Related preliminary calculations needed for computing the system matrices are also performed.

SUBROUTINE IB01ND( METH, JOBD, NOBR, M, L, R, LDR, SV, TOL, IWORK, $ DWORK, LDWORK, IWARN, INFO ) C .. Scalar Arguments .. DOUBLE PRECISION TOL INTEGER INFO, IWARN, L, LDR, LDWORK, M, NOBR CHARACTER JOBD, METH C .. Array Arguments .. DOUBLE PRECISION DWORK(*), R(LDR, *), SV(*) INTEGER IWORK(*)

**Mode Parameters**

METH CHARACTER*1 Specifies the subspace identification method to be used, as follows: = 'M': MOESP algorithm with past inputs and outputs; = 'N': N4SID algorithm. JOBD CHARACTER*1 Specifies whether or not the matrices B and D should later be computed using the MOESP approach, as follows: = 'M': the matrices B and D should later be computed using the MOESP approach; = 'N': the matrices B and D should not be computed using the MOESP approach. This parameter is not relevant for METH = 'N'.

NOBR (input) INTEGER The number of block rows, s, in the input and output block Hankel matrices. NOBR > 0. M (input) INTEGER The number of system inputs. M >= 0. L (input) INTEGER The number of system outputs. L > 0. R (input/output) DOUBLE PRECISION array, dimension ( LDR,2*(M+L)*NOBR ) On entry, the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of this array must contain the upper triangular factor R from the QR factorization of the concatenated block Hankel matrices. Denote R_ij, i,j = 1:4, the ij submatrix of R, partitioned by M*NOBR, M*NOBR, L*NOBR, and L*NOBR rows and columns. On exit, if INFO = 0, the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of this array contains the matrix S, the processed upper triangular factor R, as required by other subroutines. Specifically, let S_ij, i,j = 1:4, be the ij submatrix of S, partitioned by M*NOBR, L*NOBR, M*NOBR, and L*NOBR rows and columns. The submatrix S_22 contains the matrix of left singular vectors needed subsequently. Useful information is stored in S_11 and in the block-column S_14 : S_44. For METH = 'M' and JOBD = 'M', the upper triangular part of S_31 contains the upper triangular factor in the QR factorization of the matrix R_1c = [ R_12' R_22' R_11' ]', and S_12 contains the corresponding leading part of the transformed matrix R_2c = [ R_13' R_23' R_14' ]'. For METH = 'N', the subarray S_41 : S_43 contains the transpose of the matrix contained in S_14 : S_34. LDR INTEGER The leading dimension of the array R. LDR >= MAX( 2*(M+L)*NOBR, 3*M*NOBR ), for METH = 'M' and JOBD = 'M'; LDR >= 2*(M+L)*NOBR, for METH = 'M' and JOBD = 'N' or for METH = 'N'. SV (output) DOUBLE PRECISION array, dimension ( L*NOBR ) The singular values of the relevant part of the triangular factor from the QR factorization of the concatenated block Hankel matrices.

TOL DOUBLE PRECISION The tolerance to be used for estimating the rank of matrices. If the user sets TOL > 0, then the given value of TOL is used as a lower bound for the reciprocal condition number; an m-by-n matrix whose estimated condition number is less than 1/TOL is considered to be of full rank. If the user sets TOL <= 0, then an implicitly computed, default tolerance, defined by TOLDEF = m*n*EPS, is used instead, where EPS is the relative machine precision (see LAPACK Library routine DLAMCH). This parameter is not used for METH = 'M'.

IWORK INTEGER array, dimension ((M+L)*NOBR) This parameter is not referenced for METH = 'M'. DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK, and, for METH = 'N', DWORK(2) and DWORK(3) contain the reciprocal condition numbers of the triangular factors of the matrices U_f and r_1 [6]. On exit, if INFO = -12, DWORK(1) returns the minimum value of LDWORK. LDWORK INTEGER The length of the array DWORK. LDWORK >= max( (2*M-1)*NOBR, (M+L)*NOBR, 5*L*NOBR ), if METH = 'M' and JOBD = 'M'; LDWORK >= 5*L*NOBR, if METH = 'M' and JOBD = 'N'; LDWORK >= 5*(M+L)*NOBR+1, if METH = 'N'. For good performance, LDWORK should be larger.

IWARN INTEGER = 0: no warning; = 4: the least squares problems with coefficient matrix U_f, used for computing the weighted oblique projection (for METH = 'N'), have a rank-deficient coefficient matrix; = 5: the least squares problem with coefficient matrix r_1 [6], used for computing the weighted oblique projection (for METH = 'N'), has a rank-deficient coefficient matrix.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 2: the singular value decomposition (SVD) algorithm did not converge.

A singular value decomposition (SVD) of a certain matrix is computed, which reveals the order n of the system as the number of "non-zero" singular values. For the MOESP approach, this matrix is [ R_24' R_34' ]' := R(ms+1:(2m+l)s,(2m+l)s+1:2(m+l)s), where R is the upper triangular factor R constructed by SLICOT Library routine IB01MD. For the N4SID approach, a weighted oblique projection is computed from the upper triangular factor R and its SVD is then found.

[1] Verhaegen M., and Dewilde, P. Subspace Model Identification. Part 1: The output-error state-space model identification class of algorithms. Int. J. Control, 56, pp. 1187-1210, 1992. [2] Verhaegen M. Subspace Model Identification. Part 3: Analysis of the ordinary output-error state-space model identification algorithm. Int. J. Control, 58, pp. 555-586, 1993. [3] Verhaegen M. Identification of the deterministic part of MIMO state space models given in innovations form from input-output data. Automatica, Vol.30, No.1, pp.61-74, 1994. [4] Van Overschee, P., and De Moor, B. N4SID: Subspace Algorithms for the Identification of Combined Deterministic-Stochastic Systems. Automatica, Vol.30, No.1, pp. 75-93, 1994. [5] Van Overschee, P., and De Moor, B. Subspace Identification for Linear Systems: Theory - Implementation - Applications. Kluwer Academic Publishers, Boston/London/Dordrecht, 1996. [6] Sima, V. Subspace-based Algorithms for Multivariable System Identification. Studies in Informatics and Control, 5, pp. 335-344, 1996.

The implemented method is numerically stable. 3 The algorithm requires 0(((m+l)s) ) floating point operations.

None

**Program Text**

None

None

None