**Purpose**

To calculate the Jacobian dy/dX of the Wiener system x(t+1) = A*x(t) + B*u(t) z(t) = C*x(t) + D*u(t), y(t,i) = sum( ws(k, i)*f(w(k, i)*z(t) + b(k,i)) ) + b(k+1,i), where t = 1, 2, ..., NSMP, i = 1, 2, ..., L, k = 1, 2, ..., NN. NN is arbitrary eligible and has to be provided in IPAR(2), and X = ( wb(1), ..., wb(L), theta ) is described below. Denoting y(j) = y(1:NSMP,j), the Jacobian J has the block form dy(1)/dwb(1) 0 ..... 0 dy(1)/dtheta 0 dy(2)/dwb(2) ..... 0 dy(2)/dtheta ..... ..... ..... ..... ..... 0 ..... 0 dy(L)/dwb(L) dy(L)/dtheta but it will be returned without the zero blocks, in the form dy(1)/dwb(1) dy(1)/dtheta ... dy(L)/dwb(L) dy(L)/dtheta. dy(i)/dwb(i) depends on f and is calculated by the routine NF01BY; dy(i)/dtheta is computed by a forward-difference approximation.

SUBROUTINE NF01BD( CJTE, NSMP, M, L, IPAR, LIPAR, X, LX, U, LDU, $ E, J, LDJ, JTE, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER CJTE INTEGER INFO, L, LDJ, LDU, LDWORK, LX, LIPAR, M, NSMP C .. Array Arguments .. INTEGER IPAR(*) DOUBLE PRECISION DWORK(*), E(*), J(LDJ, *), JTE(*), U(LDU,*), $ X(*)

**Mode Parameters**

CJTE CHARACTER*1 Specifies whether the matrix-vector product J'*e should be computed or not, as follows: = 'C' : compute J'*e; = 'N' : do not compute J'*e.

NSMP (input) INTEGER The number of training samples. NSMP >= 0. M (input) INTEGER The length of each input sample. M >= 0. L (input) INTEGER The length of each output sample. L >= 0. IPAR (input/output) INTEGER array, dimension (LIPAR) On entry, the first entries of this array must contain the integer parameters needed; specifically, IPAR(1) must contain the order of the linear part, N; actually, N = abs(IPAR(1)), since setting IPAR(1) < 0 has a special meaning (see below); IPAR(2) must contain the number of neurons for the nonlinear part, NN, NN >= 0. On exit, if IPAR(1) < 0 on entry, then no computations are performed, except the needed tests on input parameters, but the following values are returned: IPAR(1) contains the length of the array J, LJ; LDJ contains the leading dimension of array J. Otherwise, IPAR(1) and LDJ are unchanged on exit. LIPAR (input) INTEGER The length of the array IPAR. LIPAR >= 2. X (input) DOUBLE PRECISION array, dimension (LX) The leading LPAR entries of this array must contain the set of system parameters, where LPAR = (NN*(L + 2) + 1)*L + N*(M + L + 1) + L*M. X has the form (wb(1), ..., wb(L), theta), where the vectors wb(i) have the structure (w(1,1), ..., w(1,L), ..., w(NN,1), ..., w(NN,L), ws(1), ..., ws(NN), b(1), ..., b(NN+1) ), and the vector theta represents the matrices A, B, C, D and x(1), and it can be retrieved from these matrices by SLICOT Library routine TB01VD and retranslated by TB01VY. LX (input) INTEGER The length of X. LX >= (NN*(L + 2) + 1)*L + N*(M + L + 1) + L*M. U (input) DOUBLE PRECISION array, dimension (LDU, M) The leading NSMP-by-M part of this array must contain the set of input samples, U = ( U(1,1),...,U(1,M); ...; U(NSMP,1),...,U(NSMP,M) ). LDU INTEGER The leading dimension of array U. LDU >= MAX(1,NSMP). E (input) DOUBLE PRECISION array, dimension (NSMP*L) If CJTE = 'C', this array must contain a vector e, which will be premultiplied with J', e = vec( Y - y ), where Y is set of output samples, and vec denotes the concatenation of the columns of a matrix. If CJTE = 'N', this array is not referenced. J (output) DOUBLE PRECISION array, dimension (LDJ, *) The leading NSMP*L-by-NCOLJ part of this array contains the Jacobian of the error function stored in a compressed form, as described above, where NCOLJ = NN*(L + 2) + 1 + N*(M + L + 1) + L*M. LDJ INTEGER The leading dimension of array J. LDJ >= MAX(1,NSMP*L). Note that LDJ is an input parameter, except for IPAR(1) < 0 on entry, when it is an output parameter. JTE (output) DOUBLE PRECISION array, dimension (LPAR) If CJTE = 'C', this array contains the matrix-vector product J'*e. If CJTE = 'N', this array is not referenced.

DWORK DOUBLE PRECISION array, dimension (LDWORK) LDWORK INTEGER The length of the array DWORK. LDWORK >= 2*NSMP*L + MAX( 2*NN, (N + L)*(N + M) + 2*N + MAX( N*(N + L), N + M + L ) ) if M > 0; LDWORK >= 2*NSMP*L + MAX( 2*NN, (N + L)*N + 2*N + MAX( N*(N + L), L ) ), if M = 0. A larger value of LDWORK could improve the efficiency.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.

BLAS routines are used for the matrix-vector multiplications, and the SLICOT Library routine TB01VY is called for the conversion of the output normal form parameters to an LTI-system; the routine NF01AD is then used for the simulation of the system with given parameters, and the routine NF01BY is called for the (analytically performed) calculation of the parts referring to the parameters of the static nonlinearity.

None

**Program Text**

None

None

None