**Purpose**

To calculate a QR factorization of the first block column and apply the orthogonal transformations (from the left) also to the second block column of a structured matrix, as follows _ _ [ R B ] [ R B ] Q' * [ ] = [ _ ] [ A C ] [ 0 C ] _ where R and R are upper triangular. The matrix A can be full or upper trapezoidal/triangular. The problem structure is exploited.

SUBROUTINE MB04OD( UPLO, N, M, P, R, LDR, A, LDA, B, LDB, C, LDC, $ TAU, DWORK ) C .. Scalar Arguments .. CHARACTER UPLO INTEGER LDA, LDB, LDC, LDR, M, N, P C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), $ R(LDR,*), TAU(*)

**Mode Parameters**

UPLO CHARACTER*1 Indicates if the matrix A is or not triangular as follows: = 'U': Matrix A is upper trapezoidal/triangular; = 'F': Matrix A is full.

N (input) INTEGER _ The order of the matrices R and R. N >= 0. M (input) INTEGER The number of columns of the matrices B and C. M >= 0. P (input) INTEGER The number of rows of the matrices A and C. P >= 0. R (input/output) DOUBLE PRECISION array, dimension (LDR,N) On entry, the leading N-by-N upper triangular part of this array must contain the upper triangular matrix R. On exit, the leading N-by-N upper triangular part of this _ array contains the upper triangular matrix R. The strict lower triangular part of this array is not referenced. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,N). A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, if UPLO = 'F', the leading P-by-N part of this array must contain the matrix A. If UPLO = 'U', the leading MIN(P,N)-by-N part of this array must contain the upper trapezoidal (upper triangular if P >= N) matrix A, and the elements below the diagonal are not referenced. On exit, the leading P-by-N part (upper trapezoidal or triangular, if UPLO = 'U') of this array contains the trailing components (the vectors v, see Method) of the elementary reflectors used in the factorization. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,P). B (input/output) DOUBLE PRECISION array, dimension (LDB,M) On entry, the leading N-by-M part of this array must contain the matrix B. On exit, the leading N-by-M part of this array contains _ the computed matrix B. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N). C (input/output) DOUBLE PRECISION array, dimension (LDC,M) On entry, the leading P-by-M part of this array must contain the matrix C. On exit, the leading P-by-M part of this array contains _ the computed matrix C. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,P). TAU (output) DOUBLE PRECISION array, dimension (N) The scalar factors of the elementary reflectors used.

DWORK DOUBLE PRECISION array, dimension (MAX(N-1,M))

The routine uses N Householder transformations exploiting the zero pattern of the block matrix. A Householder matrix has the form ( 1 ) H = I - tau *u *u', u = ( v ), i i i i i ( i) where v is a P-vector, if UPLO = 'F', or a min(i,P)-vector, if i UPLO = 'U'. The components of v are stored in the i-th column i of A, and tau is stored in TAU(i). i In-line code for applying Householder transformations is used whenever possible (see MB04OY routine).

The algorithm is backward stable.

None

**Program Text**

* MB04OD EXAMPLE PROGRAM TEXT. * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO = 0.0D0 ) INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER MMAX, NMAX, PMAX PARAMETER ( MMAX = 20, NMAX = 20, PMAX = 20 ) INTEGER LDA, LDB, LDC, LDR PARAMETER ( LDA = PMAX, LDB = NMAX, LDC = PMAX, $ LDR = NMAX ) INTEGER LDWORK PARAMETER ( LDWORK = MAX( NMAX-1,MMAX ) ) * .. Local Scalars .. CHARACTER*1 UPLO INTEGER I, J, M, N, P * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,MMAX), $ DWORK(LDWORK), R(LDR,NMAX), TAU(NMAX) * .. External Subroutines .. EXTERNAL MB04OD * .. Intrinsic Functions .. INTRINSIC MAX * .. Executable Statements .. * WRITE ( NOUT, FMT = 99999 ) * Skip the heading in the data file and read the data. READ ( NIN, FMT = '()' ) READ ( NIN, FMT = * ) N, M, P, UPLO IF ( N.LT.0 .OR. N.GT.NMAX ) THEN WRITE ( NOUT, FMT = 99993 ) N ELSE IF ( M.LT.0 .OR. M.GT.MMAX ) THEN WRITE ( NOUT, FMT = 99992 ) M ELSE IF ( P.LT.0 .OR. P.GT.PMAX ) THEN WRITE ( NOUT, FMT = 99991 ) P ELSE READ ( NIN, FMT = * ) ( ( R(I,J), J = 1,N ), I = 1,N ) READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,P ) READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N ) READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,M ), I = 1,P ) * Compute and apply QR factorization. CALL MB04OD( UPLO, N, M, P, R, LDR, A, LDA, B, LDB, C, $ LDC, TAU, DWORK ) * WRITE ( NOUT, FMT = 99997 ) DO 40 I = 1, N DO 20 J = 1, I-1 R(I,J) = ZERO 20 CONTINUE WRITE ( NOUT, FMT = 99996 ) ( R(I,J), J = 1,N ) 40 CONTINUE IF ( M.GT.0 ) THEN WRITE ( NOUT, FMT = 99995 ) DO 60 I = 1, N WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,M ) 60 CONTINUE IF ( P.GT.0 ) THEN WRITE ( NOUT, FMT = 99994 ) DO 80 I = 1, P WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,M ) 80 CONTINUE END IF END IF END IF END IF END IF STOP * 99999 FORMAT (' MB04OD EXAMPLE PROGRAM RESULTS',/1X) 99997 FORMAT (' The updated matrix R is ') 99996 FORMAT (20(1X,F10.4)) 99995 FORMAT (' The updated matrix B is ') 99994 FORMAT (' The updated matrix C is ') 99993 FORMAT (/' N is out of range.',/' N = ',I5) 99992 FORMAT (/' M is out of range.',/' M = ',I5) 99991 FORMAT (/' P is out of range.',/' P = ',I5) END

MB04OD EXAMPLE PROGRAM DATA 3 2 2 F 3. 2. 1. 0. 2. 1. 0. 0. 1. 2. 3. 1. 4. 6. 5. 3. 2. 1. 3. 3. 2. 1. 3. 3. 2.

MB04OD EXAMPLE PROGRAM RESULTS The updated matrix R is -5.3852 -6.6850 -4.6424 0.0000 -2.8828 -2.0694 0.0000 0.0000 -1.7793 The updated matrix B is -4.2710 -3.7139 -0.1555 -2.1411 -1.6021 0.9398 The updated matrix C is 0.5850 1.0141 -2.7974 -3.1162