MB04OD

QR factorization of a special structured block matrix (variant)

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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.

```
Specification
```      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(*)

```
Arguments

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.

```
Input/Output Parameters
```  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.

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

```
Method
```  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).

```
Numerical Aspects
```  The algorithm is backward stable.

```
```  None
```
Example

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
```
Program Data
``` 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.
```
Program Results
``` 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
```