## MB02FD

### Incomplete Cholesky factor of a positive definite block Toeplitz matrix

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

Purpose

```  To compute the incomplete Cholesky (ICC) factor of a symmetric
positive definite (s.p.d.) block Toeplitz matrix T, defined by
either its first block row, or its first block column, depending
on the routine parameter TYPET.

By subsequent calls of this routine, further rows / columns of
the Cholesky factor can be added.
Furthermore, the generator of the Schur complement of the leading
(P+S)*K-by-(P+S)*K block in T is available, which can be used,
e.g., for measuring the quality of the ICC factorization.

```
Specification
```      SUBROUTINE MB02FD( TYPET, K, N, P, S, T, LDT, R, LDR, DWORK,
\$                   LDWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER         TYPET
INTEGER           INFO, K, LDR, LDT, LDWORK, N, P, S
C     .. Array Arguments ..
DOUBLE PRECISION  DWORK(*), R(LDR,*), T(LDT,*)

```
Arguments

Mode Parameters

```  TYPET   CHARACTER*1
Specifies the type of T, as follows:
= 'R':  T contains the first block row of an s.p.d. block
Toeplitz matrix; the ICC factor R is upper
trapezoidal;
= 'C':  T contains the first block column of an s.p.d.
block Toeplitz matrix; the ICC factor R is lower
trapezoidal; this choice leads to better
localized memory references and hence a faster
algorithm.
Note:   in the sequel, the notation x / y means that
x corresponds to TYPET = 'R' and y corresponds to
TYPET = 'C'.

```
Input/Output Parameters
```  K       (input)  INTEGER
The number of rows / columns in T, which should be equal
to the blocksize.  K >= 0.

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

P       (input)  INTEGER
The number of previously computed block rows / columns
of R.  0 <= P <= N.

S       (input)  INTEGER
The number of block rows / columns of R to compute.
0 <= S <= N-P.

T       (input/output)  DOUBLE PRECISION array, dimension
(LDT,(N-P)*K) / (LDT,K)
On entry, if P = 0, then the leading K-by-N*K / N*K-by-K
part of this array must contain the first block row /
column of an s.p.d. block Toeplitz matrix.
If P > 0, the leading K-by-(N-P)*K / (N-P)*K-by-K must
contain the negative generator of the Schur complement of
the leading P*K-by-P*K part in T, computed from previous
calls of this routine.
On exit, if INFO = 0, then the leading K-by-(N-P)*K /
(N-P)*K-by-K part of this array contains, in the first
K-by-K block, the upper / lower Cholesky factor of
T(1:K,1:K), in the following S-1 K-by-K blocks, the
Householder transformations applied during the process,
and in the remaining part, the negative generator of the
Schur complement of the leading (P+S)*K-by(P+S)*K part
in T.

LDT     INTEGER
The leading dimension of the array T.
LDT >= MAX(1,K),        if TYPET = 'R';
LDT >= MAX(1,(N-P)*K),  if TYPET = 'C'.

R       (input/output)  DOUBLE PRECISION array, dimension
(LDR, N*K)       / (LDR, S*K )     if P = 0;
(LDR, (N-P+1)*K) / (LDR, (S+1)*K ) if P > 0.
On entry, if P > 0, then the leading K-by-(N-P+1)*K /
(N-P+1)*K-by-K part of this array must contain the
nonzero blocks of the last block row / column in the
ICC factor from a previous call of this routine. Note that
this part is identical with the positive generator of
the Schur complement of the leading P*K-by-P*K part in T.
If P = 0, then R is only an output parameter.
On exit, if INFO = 0 and P = 0, then the leading
S*K-by-N*K / N*K-by-S*K part of this array contains the
upper / lower trapezoidal ICC factor.
On exit, if INFO = 0 and P > 0, then the leading
(S+1)*K-by-(N-P+1)*K / (N-P+1)*K-by-(S+1)*K part of this
array contains the upper / lower trapezoidal part of the
P-th to (P+S)-th block rows / columns of the ICC factor.
The elements in the strictly lower / upper trapezoidal
part are not referenced.

LDR     INTEGER
The leading dimension of the array R.
LDR >= MAX(1, S*K ),        if TYPET = 'R' and P = 0;
LDR >= MAX(1, (S+1)*K ),    if TYPET = 'R' and P > 0;
LDR >= MAX(1, N*K ),        if TYPET = 'C' and P = 0;
LDR >= MAX(1, (N-P+1)*K ),  if TYPET = 'C' and P > 0.

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0,  DWORK(1)  returns the optimal
value of LDWORK.
On exit, if  INFO = -11,  DWORK(1)  returns the minimum
value of LDWORK.

LDWORK  INTEGER
The length of the array DWORK.
LDWORK >= MAX(1,(N+1)*K,4*K),   if P = 0;
LDWORK >= MAX(1,(N-P+2)*K,4*K), if P > 0.
For optimum performance LDWORK should be larger.

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  the reduction algorithm failed; the Toeplitz matrix
associated with T is not (numerically) positive
definite in its leading (P+S)*K-by-(P+S)*K part.

```
Method
```  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 is numerically stable.
3
The algorithm requires 0(K S (N-P)) floating point operations.

```
```  None
```
Example

Program Text

```*     MB02FD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          ITMAX, KMAX, NMAX
PARAMETER        ( ITMAX = 10, KMAX = 20, NMAX = 20 )
INTEGER          LDR, LDT, LDWORK
PARAMETER        ( LDR = NMAX*KMAX, LDT = KMAX,
\$                   LDWORK = ( NMAX + 1 )*KMAX )
*     .. Local Scalars ..
INTEGER          I, INFO, IT, J, K, LEN, M, N, P, PIT, POS, POSR,
\$                 S1, SCIT
CHARACTER        TYPET
DOUBLE PRECISION NNRM
*     .. Local Arrays .. (Dimensioned for TYPET = 'R'.)
INTEGER          S(ITMAX)
DOUBLE PRECISION DWORK(LDWORK), R(LDR, NMAX*KMAX),
\$                 T(LDT, NMAX*KMAX), V(NMAX*KMAX), W(NMAX*KMAX),
\$                 Z(NMAX*KMAX)
*     .. External Functions ..
LOGICAL          LSAME
DOUBLE PRECISION DNRM2
EXTERNAL         DNRM2, LSAME
*     .. External Subroutines ..
EXTERNAL         DAXPY, DCOPY, DGEMV, DLASET, DSCAL, DTRMV, MB02FD
*
*     .. Executable Statements ..
WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, K, IT
TYPET = 'R'
M = N*K
IF( N.LE.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) N
ELSE IF( K.LE.0 .OR. K.GT.KMAX ) THEN
WRITE ( NOUT, FMT = 99992 ) K
ELSE IF( IT.LE.0 .OR. IT.GT.ITMAX ) THEN
WRITE ( NOUT, FMT = 99991 ) IT
ELSE
READ ( NIN, FMT = * ) ( S(I), I = 1, IT )
READ ( NIN, FMT = * ) ( ( T(I,J), J = 1,M ), I = 1,K )
P   = 0
POS = 1
WRITE ( NOUT, FMT = 99997 )
DO 90  SCIT = 1, IT
CALL MB02FD( TYPET, K, N, P, S(SCIT), T(1,POS), LDT,
\$                   R(POS,POS), LDR, DWORK, LDWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
STOP
END IF
S1 = S(SCIT) + P
IF ( S1.EQ.0 ) THEN
*              Estimate the 2-norm of the Toeplitz matrix with 5 power
*              iterations.
LEN = N*K
CALL DLASET( 'All', LEN, 1, ONE, ONE, V, 1 )
DO 30  PIT = 1, 5
DO 10  I = 1, N
CALL DGEMV( 'NoTranspose', K, LEN-(I-1)*K, ONE, T,
\$                           LDT, V((I-1)*K+1), 1, ZERO,
\$                           W((I-1)*K+1), 1 )
10             CONTINUE
DO 20 I = 1, N-1
CALL DGEMV( 'Transpose', K, (N-I)*K, ONE,
\$                           T(1,K+1), LDT, V((I-1)*K+1), 1,
\$                           ONE, W(I*K+1), 1 )
20             CONTINUE
CALL DCOPY( LEN, W, 1, V, 1 )
NNRM = DNRM2( LEN, V, 1 )
CALL DSCAL( LEN, ONE/NNRM, V, 1 )
30          CONTINUE
ELSE
*              Estimate the 2-norm of the Schur complement with 5 power
*              iterations.
LEN = ( N - S1 )*K
CALL DLASET( 'All', LEN, 1, ONE, ONE, V, 1 )
DO 80  PIT = 1, 5
POSR = ( S1 - 1 )*K + 1
DO 40  I = 1, N - S1
CALL DGEMV( 'NoTranspose', K, LEN-(I-1)*K, ONE,
\$                           T(1,POSR+K), LDT, V((I-1)*K+1), 1,
\$                           ZERO, W((I-1)*K+1), 1 )
40             CONTINUE
DO 50  I = 1, N - S1
CALL DTRMV( 'Upper', 'NoTranspose', 'NonUnit', K,
\$                           R(POSR,POSR), LDR, V((I-1)*K+1), 1 )
CALL DGEMV( 'NoTranspose', K, LEN-I*K, ONE,
\$                           R(POSR,POSR+K), LDR, V(I*K+1), 1, ONE,
\$                           V((I-1)*K+1), 1 )
50             CONTINUE
CALL DLASET( 'All', LEN, 1, ZERO, ZERO, Z, 1 )
DO 60  I = 1, N - S1
CALL DGEMV( 'Transpose', K, LEN-I*K, ONE,
\$                           R(POSR,POSR+K), LDR, V((I-1)*K+1), 1,
\$                           ONE, Z(I*K+1), 1 )
CALL DTRMV( 'Upper', 'Transpose', 'NonUnit', K,
\$                           R(POSR,POSR), LDR, V((I-1)*K+1), 1 )
CALL DAXPY( K, ONE, V((I-1)*K+1), 1, Z((I-1)*K+1),
\$                           1 )
60             CONTINUE
CALL DLASET( 'All', LEN, 1, ZERO, ZERO, V, 1 )
DO 70  I = 1, N - S1
CALL DGEMV( 'Transpose', K, LEN-(I-1)*K, ONE,
\$                           T(1,POSR+K), LDT, W((I-1)*K+1), 1,
\$                           ONE, V((I-1)*K+1), 1 )
70             CONTINUE
CALL DAXPY( LEN, -ONE, Z, 1, V, 1 )
NNRM = DNRM2( LEN, V, 1 )
CALL DSCAL( LEN, -ONE/NNRM, V, 1 )
80          CONTINUE
POS = ( S1 - 1 )*K + 1
P   = S1
END IF
WRITE ( NOUT, FMT = 99995 ) P*K, NNRM
90    CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 100  I = 1, P*K
WRITE ( NOUT, FMT = 99994 ) ( R(I,J), J = 1, M )
100    CONTINUE
END IF
STOP
*
99999 FORMAT (' MB02FD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB02FD = ',I2)
99997 FORMAT ('   Incomplete Cholesky factorization ',
\$         //'   rows    norm(Schur complement)',/)
99996 FORMAT (/' The upper ICC factor of the block Toeplitz matrix is '
\$       )
99995 FORMAT (I4,5X,F8.4)
99994 FORMAT (20(1X,F8.4))
99993 FORMAT (/' N is out of range.',/' N = ',I5)
99992 FORMAT (/' K is out of range.',/' K = ',I5)
99991 FORMAT (/' IT is out of range.',/' IT = ',I5)
END
```
Program Data
```MB02FD EXAMPLE
4 2 3
0 1 1
3.0000    1.0000    0.1000    0.1000    0.2000    0.0500   0.2000   0.3000
1.0000    4.0000    0.4000    0.1000    0.0400    0.2000   0.1000   0.2000
```
Program Results
``` MB02FD EXAMPLE PROGRAM RESULTS

Incomplete Cholesky factorization

rows    norm(Schur complement)

0       5.5509
2       5.1590
4       4.8766

The upper ICC factor of the block Toeplitz matrix is
1.7321   0.5774   0.0577   0.0577   0.1155   0.0289   0.1155   0.1732
0.0000   1.9149   0.1915   0.0348  -0.0139   0.0957   0.0174   0.0522
0.0000   0.0000   1.7205   0.5754   0.0558   0.0465   0.1104   0.0174
0.0000   0.0000   0.0000   1.9142   0.1890   0.0357  -0.0161   0.0931
```