**Purpose**

To compute the QR factorization of the Jacobian matrix J, as received in compressed form from SLICOT Library routine NF01BD, / dy(1)/dwb(1) | dy(1)/ dtheta \ Jc = | : | : | , \ dy(L)/dwb(L) | dy(L)/ dtheta / and to apply the transformation Q on the error vector e (in-situ). The factorization is J*P = Q*R, where Q is a matrix with orthogonal columns, P a permutation matrix, and R an upper trapezoidal matrix with diagonal elements of nonincreasing magnitude for each block column (see below). The 1-norm of the scaled gradient is also returned. Actually, 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 the zero blocks are omitted. The diagonal blocks have the same size and correspond to the nonlinear part. The last block column corresponds to the linear part. It is assumed that the Jacobian matrix has at least as many rows as columns. The linear or nonlinear parts can be empty. If L <= 1, the Jacobian is represented as a full matrix.

SUBROUTINE NF01BS( N, IPAR, LIPAR, FNORM, J, LDJ, E, JNORMS, $ GNORM, IPVT, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. INTEGER INFO, LDJ, LDWORK, LIPAR, N DOUBLE PRECISION FNORM, GNORM C .. Array Arguments .. INTEGER IPAR(*), IPVT(*) DOUBLE PRECISION DWORK(*), E(*), J(*), JNORMS(*)

**Input/Output Parameters**

N (input) INTEGER The number of columns of the Jacobian matrix J. N = BN*BSN + ST >= 0. (See parameter description below.) IPAR (input) INTEGER array, dimension (LIPAR) The integer parameters describing the structure of the matrix J, as follows: IPAR(1) must contain ST, the number of parameters corresponding to the linear part. ST >= 0. IPAR(2) must contain BN, the number of blocks, BN = L, for the parameters corresponding to the nonlinear part. BN >= 0. IPAR(3) must contain BSM, the number of rows of the blocks J_k = dy(k)/dwb(k), k = 1:BN, if BN > 0, or the number of rows of the matrix J, if BN <= 1. BN*BSM >= N, if BN > 0; BSM >= N, if BN = 0. IPAR(4) must contain BSN, the number of columns of the blocks J_k, k = 1:BN. BSN >= 0. LIPAR (input) INTEGER The length of the array IPAR. LIPAR >= 4. FNORM (input) DOUBLE PRECISION The Euclidean norm of the vector e. FNORM >= 0. J (input/output) DOUBLE PRECISION array, dimension (LDJ, NC) where NC = N if BN <= 1, and NC = BSN+ST, if BN > 1. On entry, the leading NR-by-NC part of this array must contain the (compressed) representation (Jc) of the Jacobian matrix J, where NR = BSM if BN <= 1, and NR = BN*BSM, if BN > 1. On exit, the leading N-by-NC part of this array contains a (compressed) representation of the upper triangular factor R of the Jacobian matrix. The matrix R has the same structure as the Jacobian matrix J, but with an additional diagonal block. Note that for efficiency of the later calculations, the matrix R is delivered with the leading dimension MAX(1,N), possibly much smaller than the value of LDJ on entry. LDJ (input/output) INTEGER The leading dimension of array J. On entry, LDJ >= MAX(1,NR). On exit, LDJ >= MAX(1,N). E (input/output) DOUBLE PRECISION array, dimension (NR) On entry, this array contains the vector e, e = vec( Y - y ), where Y is set of output samples, and vec denotes the concatenation of the columns of a matrix. On exit, this array contains the updated vector Z*Q'*e, where Z is the block row permutation matrix used in the QR factorization of J (see METHOD). JNORMS (output) DOUBLE PRECISION array, dimension (N) This array contains the Euclidean norms of the columns of the Jacobian matrix, considered in the initial order. GNORM (output) DOUBLE PRECISION If FNORM > 0, the 1-norm of the scaled vector J'*e/FNORM, with each element i further divided by JNORMS(i) (if JNORMS(i) is nonzero). If FNORM = 0, the returned value of GNORM is 0. IPVT (output) INTEGER array, dimension (N) This array defines the permutation matrix P such that J*P = Q*R. Column j of P is column IPVT(j) of the identity matrix.

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 >= 1, if N = 0 or BN <= 1 and BSM = N = 1; otherwise, LDWORK >= 4*N+1, if BN <= 1 or BSN = 0; LDWORK >= JWORK, if BN > 1 and BSN > 0, where JWORK is given by the following procedure: JWORK = BSN + MAX(3*BSN+1,ST); JWORK = MAX(JWORK,4*ST+1), if BSM > BSN; JWORK = MAX(JWORK,(BSM-BSN)*(BN-1)), if BSN < BSM < 2*BSN. For optimum performance LDWORK should be larger.

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

A QR factorization with column pivoting of the matrix J is computed, J*P = Q*R. If l = L > 1, the R factor of the QR factorization has the same structure as the Jacobian, but with an additional diagonal block. Denote / J_1 0 .. 0 | L_1 \ | 0 J_2 .. 0 | L_2 | J = | : : .. : | : | . | : : .. : | : | \ 0 0 .. J_l | L_l / The algorithm consists in two phases. In the first phase, the algorithm uses QR factorizations with column pivoting for each block J_k, k = 1:l, and applies the orthogonal matrix Q'_k to the corresponding part of the last block column and of e. After all block rows have been processed, the block rows are interchanged so that the zeroed submatrices in the first l block columns are moved to the bottom part. The same block row permutation Z is also applied to the vector e. At the end of the first phase, the structure of the processed matrix J is / R_1 0 .. 0 | L^1_1 \ | 0 R_2 .. 0 | L^1_2 | | : : .. : | : | . | : : .. : | : | | 0 0 .. R_l | L^1_l | | 0 0 .. 0 | L^2_1 | | : : .. : | : | \ 0 0 .. 0 | L^2_l / In the second phase, the submatrix L^2_1:l is triangularized using an additional QR factorization with pivoting. (The columns of L^1_1:l are also permuted accordingly.) Therefore, the column pivoting is restricted to each such local block column. If l <= 1, the matrix J is triangularized in one phase, by one QR factorization with pivoting. In this case, the column pivoting is global.

None

**Program Text**

None

None

None