NF01BS

QR factorization with column pivoting for Wiener system identification

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