MB02KD

Computation of the product C = alpha op( T ) B + beta C, with T a block Toeplitz matrix

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

Purpose

  To compute the matrix product

            C = alpha*op( T )*B + beta*C,

  where alpha and beta are scalars and T is a block Toeplitz matrix
  specified by its first block column TC and first block row TR;
  B and C are general matrices of appropriate dimensions.

Specification
      SUBROUTINE MB02KD( LDBLK, TRANS, K, L, M, N, R, ALPHA, BETA,
     $                   TC, LDTC, TR, LDTR, B, LDB, C, LDC, DWORK,
     $                   LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER         LDBLK, TRANS
      INTEGER           INFO, K, L, LDB, LDC, LDTC, LDTR, LDWORK, M, N,
     $                  R
      DOUBLE PRECISION  ALPHA, BETA
C     .. Array Arguments ..
      DOUBLE PRECISION  B(LDB,*), C(LDC,*), DWORK(*), TC(LDTC,*),
     $                  TR(LDTR,*)

Arguments

Mode Parameters

  LDBLK   CHARACTER*1
          Specifies where the (1,1)-block of T is stored, as
          follows:
          = 'C':  in the first block of TC;
          = 'R':  in the first block of TR.

  TRANS   CHARACTER*1
          Specifies the form of op( T ) to be used in the matrix
          multiplication as follows:
          = 'N':  op( T ) = T;
          = 'T':  op( T ) = T';
          = 'C':  op( T ) = T'.

Input/Output Parameters
  K       (input) INTEGER
          The number of rows in the blocks of T.  K >= 0.

  L       (input) INTEGER
          The number of columns in the blocks of T.  L >= 0.

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

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

  R       (input) INTEGER
          The number of columns in B and C.  R >= 0.

  ALPHA   (input) DOUBLE PRECISION
          The scalar alpha. When alpha is zero then TC, TR and B
          are not referenced.

  BETA    (input) DOUBLE PRECISION
          The scalar beta. When beta is zero then C need not be set
          before entry.

  TC      (input)  DOUBLE PRECISION array, dimension (LDTC,L)
          On entry with LDBLK = 'C', the leading M*K-by-L part of
          this array must contain the first block column of T.
          On entry with LDBLK = 'R', the leading (M-1)*K-by-L part
          of this array must contain the 2nd to the M-th blocks of
          the first block column of T.

  LDTC    INTEGER
          The leading dimension of the array TC.
          LDTC >= MAX(1,M*K),      if LDBLK = 'C';
          LDTC >= MAX(1,(M-1)*K),  if LDBLK = 'R'.

  TR      (input)  DOUBLE PRECISION array, dimension (LDTR,k)
          where k is (N-1)*L when LDBLK = 'C' and is N*L when
          LDBLK = 'R'.
          On entry with LDBLK = 'C', the leading K-by-(N-1)*L part
          of this array must contain the 2nd to the N-th blocks of
          the first block row of T.
          On entry with LDBLK = 'R', the leading K-by-N*L part of
          this array must contain the first block row of T.

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

  B       (input)  DOUBLE PRECISION array, dimension (LDB,R)
          On entry with TRANS = 'N', the leading N*L-by-R part of
          this array must contain the matrix B.
          On entry with TRANS = 'T' or TRANS = 'C', the leading
          M*K-by-R part of this array must contain the matrix B.

  LDB     INTEGER
          The leading dimension of the array B.
          LDB >= MAX(1,N*L),  if TRANS = 'N';
          LDB >= MAX(1,M*K),  if TRANS = 'T' or TRANS = 'C'.

  C       (input/output)  DOUBLE PRECISION array, dimension (LDC,R)
          On entry with TRANS = 'N', the leading M*K-by-R part of
          this array must contain the matrix C.
          On entry with TRANS = 'T' or TRANS = 'C', the leading
          N*L-by-R part of this array must contain the matrix C.
          On exit with TRANS = 'N', the leading M*K-by-R part of
          this array contains the updated matrix C.
          On exit with TRANS = 'T' or TRANS = 'C', the leading
          N*L-by-R part of this array contains the updated matrix C.

  LDC     INTEGER
          The leading dimension of the array C.
          LDC >= MAX(1,M*K),  if TRANS = 'N';
          LDC >= MAX(1,N*L),  if TRANS = 'T' or TRANS = 'C'.

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

  LDWORK  INTEGER
          The length of the array DWORK.  LDWORK >= 1.
          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.

Method
  For point Toeplitz matrices or sufficiently large block Toeplitz
  matrices, this algorithm uses convolution algorithms based on
  the fast Hartley transforms [1]. Otherwise, TC is copied in
  reversed order into the workspace such that C can be computed from
  barely M matrix-by-matrix multiplications.

References
  [1] Van Loan, Charles.
      Computational frameworks for the fast Fourier transform.
      SIAM, 1992.

Numerical Aspects
  The algorithm requires O( (K*L+R*L+K*R)*(N+M)*log(N+M) + K*L*R )
  floating point operations.

Further Comments
  None
Example

Program Text

*     MB02KD EXAMPLE PROGRAM TEXT
*     RELEASE 5.0, SLICOT COPYRIGHT 2001.
*
*     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          KMAX, LMAX, MMAX, NMAX, RMAX
      PARAMETER        ( KMAX = 20, LMAX = 20, MMAX = 20, NMAX = 20,
     $                   RMAX = 20 )
      INTEGER          LDB, LDC, LDTC, LDTR, LDWORK
      PARAMETER        ( LDB  = LMAX*NMAX, LDC  = KMAX*MMAX,
     $                   LDTC = MMAX*KMAX, LDTR = KMAX,
     $                   LDWORK = 2*( KMAX*LMAX + KMAX*RMAX
     $                            + LMAX*RMAX + 1 )*( MMAX + NMAX ) )
*     .. Local Scalars .. 
      INTEGER          I, INFO, J, K, L, M, N, R
      CHARACTER        LDBLK, TRANS
      DOUBLE PRECISION ALPHA, BETA
*     .. Local Arrays .. (Dimensioned for TRANS = 'N'.)
      DOUBLE PRECISION B(LDB,RMAX), C(LDC,RMAX), DWORK(LDWORK),
     $                 TC(LDTC,LMAX), TR(LDTR,NMAX*LMAX)
*     .. External Functions ..
      LOGICAL          LSAME 
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         MB02KD
*     .. 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 = * )  K, L, M, N, R, LDBLK, TRANS
      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 IF( R.LE.0 .OR. R.GT.RMAX ) THEN
         WRITE ( NOUT, FMT = 99990 ) N
      ELSE
         IF ( LSAME( LDBLK, 'R' ) ) THEN
            READ ( NIN, FMT = * ) ( ( TC(I,J), J = 1,L ),
     $                              I = 1,(M-1)*K )
            READ ( NIN, FMT = * ) ( ( TR(I,J), J = 1,N*L ), I = 1,K )      
         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 )
         END IF
         IF ( LSAME( TRANS, 'N' ) ) THEN
            READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,R ), I = 1,N*L )
         ELSE
            READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,R ), I = 1,M*K )
         END IF
         ALPHA = ONE
         BETA  = ZERO
         CALL MB02KD( LDBLK, TRANS, K, L, M, N, R, ALPHA, BETA, TC,
     $                LDTC, TR, LDTR, B, LDB, C, LDC, DWORK, LDWORK,
     $                INFO )
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            IF ( LSAME( TRANS, 'N' ) ) THEN      
               WRITE ( NOUT, FMT = 99997 )
               DO 10  I = 1, M*K
                  WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,R )
   10          CONTINUE      
            ELSE
               WRITE ( NOUT, FMT = 99996 )
               DO 20  I = 1, N*L
                  WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,R )
   20          CONTINUE
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' MB02KD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB02KD = ',I2)
99997 FORMAT (' The product C = T * B is ')
99996 FORMAT (' The product C = T^T * B 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)
99990 FORMAT (/' R is out of range.',/' R = ',I5)
      END
Program Data
MB02KD EXAMPLE PROGRAM DATA
   3    2   4    5    1    C    N
     4.0     1.0
     3.0     5.0
     2.0     1.0
     4.0     1.0
     3.0     4.0
     2.0     4.0
     3.0     1.0
     3.0     0.0
     4.0     4.0
     5.0     1.0
     3.0     1.0
     4.0     3.0
     5.0     2.0     2.0     2.0     2.0     1.0     1.0     3.0
     4.0     1.0     5.0     4.0     5.0     4.0     1.0     2.0
     2.0     3.0     4.0     1.0     3.0     3.0     3.0     3.0
     0.0
     2.0
     2.0
     2.0
     1.0
     3.0
     3.0
     4.0
     2.0
     3.0
Program Results
 MB02KD EXAMPLE PROGRAM RESULTS

 The product C = T * B is 
  45.0000
  76.0000
  55.0000
  44.0000
  84.0000
  56.0000
  52.0000
  70.0000
  54.0000
  49.0000
  63.0000
  59.0000

Return to index