**Purpose**

To solve one of the systems of linear equations R*x = b , or R'*x = b , in the least squares sense, where R is an n-by-n block upper triangular matrix, with the structure / R_1 0 .. 0 | L_1 \ | 0 R_2 .. 0 | L_2 | | : : .. : | : | , | 0 0 .. R_l | L_l | \ 0 0 .. 0 | R_l+1 / with the upper triangular submatrices R_k, k = 1:l+1, square, and the first l of the same order, BSN. The diagonal elements of each block R_k have nonincreasing magnitude. The matrix R is stored in the compressed form, as returned by SLICOT Library routine NF01BS, / R_1 | L_1 \ | R_2 | L_2 | Rc = | : | : | , | R_l | L_l | \ X | R_l+1 / where the submatrix X is irrelevant. If the matrix R does not have full rank, then a least squares solution is obtained. If l <= 1, then R is an upper triangular matrix and its full upper triangle is stored. Optionally, the transpose of the matrix R can be stored in the strict lower triangles of the submatrices R_k, k = 1:l+1, and in the arrays SDIAG and S, as described at the parameter UPLO below.

SUBROUTINE NF01BR( COND, UPLO, TRANS, N, IPAR, LIPAR, R, LDR, $ SDIAG, S, LDS, B, RANKS, TOL, DWORK, LDWORK, $ INFO ) C .. Scalar Arguments .. CHARACTER COND, TRANS, UPLO INTEGER INFO, LDR, LDS, LDWORK, LIPAR, N DOUBLE PRECISION TOL C .. Array Arguments .. INTEGER IPAR(*), RANKS(*) DOUBLE PRECISION B(*), DWORK(*), R(LDR,*), S(LDS,*), SDIAG(*)

**Mode Parameters**

COND CHARACTER*1 Specifies whether the condition of submatrices R_k should be estimated, as follows: = 'E' : use incremental condition estimation and store the numerical rank of R_k in the array entry RANKS(k), for k = 1:l+1; = 'N' : do not use condition estimation, but check the diagonal entries of R_k for zero values; = 'U' : use the ranks already stored in RANKS(1:l+1). UPLO CHARACTER*1 Specifies the storage scheme for the matrix R, as follows: = 'U' : the upper triangular part is stored as in Rc; = 'L' : the lower triangular part is stored, namely, - the transpose of the strict upper triangle of R_k is stored in the strict lower triangle of R_k, for k = 1:l+1; - the diagonal elements of R_k, k = 1:l+1, are stored in the array SDIAG; - the transpose of the last block column in R (without R_l+1) is stored in the array S. TRANS CHARACTER*1 Specifies the form of the system of equations, as follows: = 'N': R*x = b (No transpose); = 'T': R'*x = b (Transpose); = 'C': R'*x = b (Transpose).

N (input) INTEGER The order of the matrix R. N = BN*BSN + ST >= 0. (See parameter description below.) IPAR (input) INTEGER array, dimension (LIPAR) The integer parameters describing the structure of the matrix R, as follows: IPAR(1) must contain ST, the number of columns of the submatrices L_k and the order of R_l+1. ST >= 0. IPAR(2) must contain BN, the number of blocks, l, in the block diagonal part of R. BN >= 0. IPAR(3) must contain BSM, the number of rows of the blocks R_k, k = 1:l. BSM >= 0. IPAR(4) must contain BSN, the number of columns of the blocks R_k, k = 1:l. BSN >= 0. BSM is not used by this routine, but assumed equal to BSN. LIPAR (input) INTEGER The length of the array IPAR. LIPAR >= 4. R (input) DOUBLE PRECISION array, dimension (LDR, NC) where NC = N if BN <= 1, and NC = BSN+ST, if BN > 1. If UPLO = 'U', the leading N-by-NC part of this array must contain the (compressed) representation (Rc) of the upper triangular matrix R. The submatrix X in Rc and the strict lower triangular parts of the diagonal blocks R_k, k = 1:l+1, are not referenced. If BN <= 1 or BSN = 0, then the full upper triangle of R must be stored. If UPLO = 'L', BN > 1 and BSN > 0, the leading (N-ST)-by-BSN part of this array must contain the transposes of the strict upper triangles of R_k, k = 1:l, stored in the strict lower triangles of R_k, and the strict lower triangle of R_l+1 must contain the transpose of the strict upper triangle of R_l+1. The submatrix X in Rc is not referenced. The diagonal elements of R_k, and, if COND = 'E', the upper triangular parts of R_k, k = 1:l+1, are modified internally, but are restored on exit. If UPLO = 'L' and BN <= 1 or BSN = 0, the leading N-by-N strict lower triangular part of this array must contain the transpose of the strict upper triangular part of R. The diagonal elements and, if COND = 'E', the upper triangular elements are modified internally, but are restored on exit. LDR INTEGER The leading dimension of the array R. LDR >= MAX(1,N). SDIAG (input) DOUBLE PRECISION array, dimension (N) If UPLO = 'L', this array must contain the diagonal entries of R_k, k = 1:l+1. This array is modified internally, but is restored on exit. This parameter is not referenced if UPLO = 'U'. S (input) DOUBLE PRECISION array, dimension (LDS,N-ST) If UPLO = 'L', BN > 1, and BSN > 0, the leading ST-by-(N-ST) part of this array must contain the transpose of the rectangular part of the last block column in R, that is [ L_1' L_2' ... L_l' ] . If COND = 'E', S is modified internally, but is restored on exit. This parameter is not referenced if UPLO = 'U', or BN <= 1, or BSN = 0. LDS INTEGER The leading dimension of the array S. LDS >= 1, if UPLO = 'U', or BN <= 1, or BSN = 0; LDS >= MAX(1,ST), if UPLO = 'L', BN > 1, and BSN > 0. B (input/output) DOUBLE PRECISION array, dimension (N) On entry, this array must contain the right hand side vector b. On exit, this array contains the (least squares) solution of the system R*x = b or R'*x = b. RANKS (input or output) INTEGER array, dimension (r), where r = BN + 1, if ST > 0, BSN > 0, and BN > 1; r = BN, if ST = 0 and BSN > 0; r = 1, if ST > 0 and ( BSN = 0 or BN <= 1 ); r = 0, if ST = 0 and BSN = 0. On entry, if COND = 'U' and N > 0, this array must contain the numerical ranks of the submatrices R_k, k = 1:l(+1). On exit, if COND = 'E' or 'N' and N > 0, this array contains the numerical ranks of the submatrices R_k, k = 1:l(+1), estimated according to the value of COND.

TOL DOUBLE PRECISION If COND = 'E', the tolerance to be used for finding the ranks of the submatrices R_k. If the user sets TOL > 0, then the given value of TOL is used as a lower bound for the reciprocal condition number; a (sub)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 = N*EPS, is used instead, where EPS is the machine precision (see LAPACK Library routine DLAMCH). This parameter is not relevant if COND = 'U' or 'N'.

DWORK DOUBLE PRECISION array, dimension (LDWORK) LDWORK INTEGER The length of the array DWORK. Denote Full = ( BN <= 1 or BSN = 0 ); Comp = ( BN > 1 and BSN > 0 ). LDWORK >= 2*N, if Full and COND = 'E'; LDWORK >= 2*MAX(BSN,ST), if Comp and COND = 'E'; LDWORK >= 0, in the remaining cases.

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

Block back or forward substitution is used (depending on TRANS and UPLO), exploiting the special structure and storage scheme of the matrix R. If a submatrix R_k, k = 1:l+1, is singular, a local basic least squares solution is computed. Therefore, the returned result is not the basic least squares solution for the whole problem, but a concatenation of (least squares) solutions of the individual subproblems involving R_k, k = 1:l+1 (with adapted right hand sides).

2 2 The algorithm requires 0(BN*BSN + ST + N*ST) operations and is backward stable, if R is nonsingular.

None

**Program Text**

None

None

None