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 [1], [2].

References
  [1] Kailath, T. and Sayed, A.
      Fast Reliable Algorithms for Matrices with Structure.
      SIAM Publications, Philadelphia, 1999.

  [2] 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.

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

Return to index