## MB02JD

### Full QR factorization of a block Toeplitz matrix of full rank

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

Purpose

```  To compute a lower triangular matrix R and a matrix Q with
Q^T Q = I such that
T
T  =  Q R ,

where T is a K*M-by-L*N block Toeplitz matrix with blocks of size
(K,L). The first column of T will be denoted by TC and the first
row by TR. It is assumed that the first MIN(M*K, N*L) columns of T
have full rank.

By subsequent calls of this routine the factors Q and R can be
computed block column by block column.

```
Specification
```      SUBROUTINE MB02JD( JOB, K, L, M, N, P, S, TC, LDTC, TR, LDTR, Q,
\$                   LDQ, R, LDR, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER         JOB
INTEGER           INFO, K, L, LDQ, LDR, LDTC, LDTR, LDWORK,
\$                  M, N, P, S
C     .. Array Arguments ..
DOUBLE PRECISION  DWORK(LDWORK), Q(LDQ,*), R(LDR,*), TC(LDTC,*),
\$                  TR(LDTR,*)

```
Arguments

Mode Parameters

```  JOB     CHARACTER*1
Specifies the output of the routine as follows:
= 'Q':  computes Q and R;
= 'R':  only computes R.

```
Input/Output Parameters
```  K       (input)  INTEGER
The number of rows in one block of T.  K >= 0.

L       (input)  INTEGER
The number of columns in one block of T.  L >= 0.

M       (input)  INTEGER
The number of blocks in one block column of T.  M >= 0.

N       (input)  INTEGER
The number of blocks in one block row of T.  N >= 0.

P       (input)  INTEGER
The number of previously computed block columns of R.
P*L < MIN( M*K,N*L ) + L and P >= 0.

S       (input)  INTEGER
The number of block columns of R to compute.
(P+S)*L < MIN( M*K,N*L ) + L and S >= 0.

TC      (input) DOUBLE PRECISION array, dimension (LDTC, L)
On entry, if P = 0, the leading M*K-by-L part of this
array must contain the first block column of T.

LDTC    INTEGER
The leading dimension of the array TC.
LDTC >= MAX(1,M*K).

TR      (input)  DOUBLE PRECISION array, dimension (LDTR,(N-1)*L)
On entry, if P = 0, the leading K-by-(N-1)*L part of this
array must contain the first block row of T without the

LDTR    INTEGER
The leading dimension of the array TR.
LDTR >= MAX(1,K).

Q       (input/output)  DOUBLE PRECISION array, dimension
(LDQ,MIN( S*L, MIN( M*K,N*L )-P*L ))
On entry, if JOB = 'Q'  and  P > 0, the leading M*K-by-L
part of this array must contain the last block column of Q
from a previous call of this routine.
On exit, if JOB = 'Q'  and  INFO = 0, the leading
M*K-by-MIN( S*L, MIN( M*K,N*L )-P*L ) part of this array
contains the P-th to (P+S)-th block columns of the factor
Q.

LDQ     INTEGER
The leading dimension of the array Q.
LDQ >= MAX(1,M*K), if JOB = 'Q';
LDQ >= 1,          if JOB = 'R'.

R       (input/output)  DOUBLE PRECISION array, dimension
(LDR,MIN( S*L, MIN( M*K,N*L )-P*L ))
On entry, if P > 0, the leading (N-P+1)*L-by-L
part of this array must contain the nozero part of the
last block column of R from a previous call of this
routine.
One exit, if INFO = 0, the leading
MIN( N, N-P+1 )*L-by-MIN( S*L, MIN( M*K,N*L )-P*L )
part of this array contains the nonzero parts of the P-th
to (P+S)-th block columns of the lower triangular
factor R.
Note that elements in the strictly upper triangular part
will not be referenced.

LDR     INTEGER
The leading dimension of the array R.
LDR >= MAX( 1, MIN( N, N-P+1 )*L )

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
On exit, if INFO = -17,  DWORK(1) returns the minimum
value of LDWORK.
If JOB = 'Q', the first 1 + ( (N-1)*L + M*K )*( 2*K + L )
elements of DWORK should be preserved during successive
calls of the routine.
If JOB = 'R', the first 1 + (N-1)*L*( 2*K + L ) elements
of DWORK should be preserved during successive calls of
the routine.

LDWORK  INTEGER
The length of the array DWORK.
JOB = 'Q':
LDWORK >= 1 + ( M*K + ( N - 1 )*L )*( L + 2*K ) + 6*L
+ MAX( M*K,( N - MAX( 1,P )*L ) );
JOB = 'R':
If P = 0,
LDWORK >= MAX( 1 + ( N - 1 )*L*( L + 2*K ) + 6*L
+ (N-1)*L, M*K*( L + 1 ) + L );
If P > 0,
LDWORK >= 1 + (N-1)*L*( L + 2*K ) + 6*L + (N-P)*L.
For optimum performance LDWORK should be larger.

If LDWORK = -1, then a workspace query is assumed;
the routine only calculates the optimal size of the
DWORK array, returns this value as the first entry of
the DWORK array, and no error message related to LDWORK
is issued by XERBLA.

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  the full rank condition for the first MIN(M*K, N*L)
columns of T is (numerically) violated.

```
Method
```  Block Householder transformations and modified hyperbolic
rotations are used in the Schur algorithm , .

```
References
```   Kailath, T. and Sayed, A.
Fast Reliable Algorithms for Matrices with Structure.

 Kressner, D. and Van Dooren, P.
Factorizations and linear system solvers for matrices with
Toeplitz structure.
SLICOT Working Note 2000-2, 2000.

```
Numerical Aspects
```  The implemented method yields a factor R which has comparable
accuracy with the Cholesky factor of T^T * T. Q is implicitly
computed from the formula Q = T * inv(R^T R) * R, i.e., for ill
conditioned problems this factor is of very limited value.
2
The algorithm requires 0(K*L *M*N) floating point operations.

```
```  None
```
Example

Program Text

```*     MB02JD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          KMAX, LMAX, MMAX, NMAX
PARAMETER        ( KMAX = 10, LMAX = 10, MMAX = 20, NMAX = 20 )
INTEGER          LDR, LDQ, LDTC, LDTR, LDWORK
PARAMETER        ( LDR  = NMAX*LMAX, LDQ  = MMAX*KMAX,
\$                   LDTC = MMAX*KMAX, LDTR = KMAX,
\$                   LDWORK = ( MMAX*KMAX + NMAX*LMAX )
\$                            *( LMAX + 2*KMAX ) + 6*LMAX
\$                            + MMAX*KMAX + NMAX*LMAX )
*     .. Local Scalars ..
INTEGER          I, INFO, J, K, L, M, N, S
CHARACTER        JOB
*     .. Local Arrays ..
DOUBLE PRECISION DWORK(LDWORK), Q(LDQ,NMAX*LMAX),
\$                 R(LDR,NMAX*LMAX), TC(LDTC,LMAX),
\$                 TR(LDTR,NMAX*LMAX)
*     .. External Functions ..
LOGICAL          LSAME
EXTERNAL         LSAME
*     .. External Subroutines ..
EXTERNAL         MB02JD
*     .. Intrinsic Functions ..
INTRINSIC        MIN
*
*     .. Executable Statements ..
WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) K, L, M, N, JOB
IF( K.LE.0 .OR. K.GT.KMAX ) THEN
WRITE ( NOUT, FMT = 99994 ) K
ELSE IF( L.LE.0 .OR. L.GT.LMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) L
ELSE IF( M.LE.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99992 ) M
ELSE IF( N.LE.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99991 ) N
ELSE
READ ( NIN, FMT = * ) ( ( TC(I,J), J = 1,L ), I = 1,M*K )
READ ( NIN, FMT = * ) ( ( TR(I,J), J = 1,( N - 1 )*L ),
\$                             I = 1,K )
S = ( MIN( M*K, N*L ) + L - 1 ) / L
*        Compute the required part of the QR factorization.
CALL MB02JD( JOB, K, L, M, N, 0, S, TC, LDTC, TR, LDTR, Q, LDQ,
\$                R, LDR, DWORK, LDWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
IF ( LSAME( JOB, 'Q' ) ) THEN
WRITE ( NOUT, FMT = 99997 )
DO 10  I = 1, M*K
WRITE ( NOUT, FMT = 99995 )
\$                  ( Q(I,J), J = 1, MIN( N*L, M*K ) )
10          CONTINUE
END IF
WRITE ( NOUT, FMT = 99996 )
DO 20  I = 1, N*L
WRITE ( NOUT, FMT = 99995 )
\$               ( R(I,J), J = 1, MIN( N*L, M*K ) )
20       CONTINUE
END IF
END IF
*
STOP
*
99999 FORMAT (' MB02JD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB02JD = ',I2)
99997 FORMAT (/' The factor Q is ')
99996 FORMAT (/' The factor R is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (/' K is out of range.',/' K = ',I5)
99993 FORMAT (/' L is out of range.',/' L = ',I5)
99992 FORMAT (/' M is out of range.',/' M = ',I5)
99991 FORMAT (/' N is out of range.',/' N = ',I5)
END
```
Program Data
```MB02JD EXAMPLE PROGRAM DATA
2   3    4    3    Q
1.0     4.0     0.0
4.0     1.0     2.0
4.0     2.0     2.0
5.0     3.0     2.0
2.0     4.0     4.0
5.0     3.0     4.0
2.0     2.0     5.0
4.0     2.0     3.0
3.0     4.0     2.0     5.0     0.0     4.0
5.0     1.0     1.0     2.0     4.0     1.0
```
Program Results
``` MB02JD EXAMPLE PROGRAM RESULTS

The factor Q is
-0.0967   0.7166  -0.4651   0.1272   0.4357   0.0435   0.2201   0.0673
-0.3867  -0.3108  -0.0534   0.5251   0.0963  -0.3894   0.1466   0.5412
-0.3867  -0.0990  -0.1443  -0.7021   0.3056  -0.3367  -0.3233   0.1249
-0.4834  -0.0178  -0.3368  -0.1763  -0.5446   0.5100   0.1503   0.2054
-0.1933   0.5859   0.3214   0.1156  -0.4670  -0.3199  -0.4185   0.0842
-0.4834  -0.0178   0.1072   0.0357  -0.0575  -0.2859   0.4339  -0.6928
-0.1933   0.1623   0.7251  -0.1966   0.2736   0.3058   0.3398   0.2968
-0.3867  -0.0990   0.0777   0.3615   0.3386   0.4421  -0.5693  -0.2641

The factor R is
-10.3441   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000
-6.3805   4.7212   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000
-7.3472   1.9320   4.5040   0.0000   0.0000   0.0000   0.0000   0.0000
-10.0541   2.5101   0.5065   3.6550   0.0000   0.0000   0.0000   0.0000
-6.5738   3.6127   1.2702  -1.3146   3.5202   0.0000   0.0000   0.0000
-5.2204   2.4764   2.4113   1.3890   1.2780   2.4976   0.0000   0.0000
-9.6674   3.2445  -0.5099  -0.0224   2.6548   2.9491   1.0049   0.0000
-6.3805   0.6968   1.9483   0.3050   0.7002  -2.0220  -2.8246   2.3147
-4.1570   2.4309  -0.7190  -0.1455   3.0149   0.5454   0.9394  -0.0548
```