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.Further Comments
NoneExample
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