MB05OD

Matrix exponential for a real matrix, with accuracy estimate

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

Purpose

  To compute exp(A*delta) where A is a real N-by-N matrix and delta
  is a scalar value. The routine also returns the minimal number of
  accurate digits in the 1-norm of exp(A*delta) and the number of
  accurate digits in the 1-norm of exp(A*delta) at 95% confidence
  level.

Specification
      SUBROUTINE MB05OD( BALANC, N, NDIAG, DELTA, A, LDA, MDIG, IDIG,
     $                   IWORK, DWORK, LDWORK, IWARN, INFO )
C     .. Scalar Arguments ..
      CHARACTER         BALANC
      INTEGER           IDIG, INFO, IWARN, LDA, LDWORK, MDIG, N,
     $                  NDIAG
      DOUBLE PRECISION  DELTA
C     .. Array Arguments ..
      INTEGER           IWORK(*)
      DOUBLE PRECISION  A(LDA,*), DWORK(*)

Arguments

Mode Parameters

  BALANC  CHARACTER*1
          Specifies whether or not a balancing transformation (done
          by SLICOT Library routine MB04MD) is required, as follows:
          = 'N', do not use balancing;
          = 'S', use balancing (scaling).

Input/Output Parameters
  N       (input) INTEGER
          The order of the matrix A.  N >= 0.

  NDIAG   (input) INTEGER
          The specified order of the diagonal Pade approximant.
          In the absence of further information NDIAG should
          be set to 9.  NDIAG should not exceed 15.  NDIAG >= 1.

  DELTA   (input) DOUBLE PRECISION
          The scalar value delta of the problem.

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On input, the leading N-by-N part of this array must
          contain the matrix A of the problem. (This is not needed
          if DELTA = 0.)
          On exit, if INFO = 0, the leading N-by-N part of this
          array contains the solution matrix exp(A*delta).

  LDA     INTEGER
          The leading dimension of array A.  LDA >= MAX(1,N).

  MDIG    (output) INTEGER
          The minimal number of accurate digits in the 1-norm of
          exp(A*delta).

  IDIG    (output) INTEGER
          The number of accurate digits in the 1-norm of
          exp(A*delta) at 95% confidence level.

Workspace
  IWORK   INTEGER array, dimension (N)

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK >= N*(2*N+NDIAG+1)+NDIAG, if N >  1.
          LDWORK >= 1,                     if N <= 1.

Warning Indicator
  IWARN   INTEGER
          = 0:  no warning;
          = 1:  if MDIG = 0 and IDIG > 0, warning for possible
                inaccuracy (the exponential has been computed);
          = 2:  if MDIG = 0 and IDIG = 0, warning for severe
                inaccuracy (the exponential has been computed);
          = 3:  if balancing has been requested, but it failed to
                reduce the matrix norm and was not actually used.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if the norm of matrix A*delta (after a possible
                balancing) is too large to obtain an accurate
                result;
          = 2:  if the coefficient matrix (the denominator of the
                Pade approximant) is exactly singular; try a
                different value of NDIAG;
          = 3:  if the solution exponential would overflow, possibly
                due to a too large value DELTA; the calculations
                stopped prematurely. This error is not likely to
                appear.

Method
  The exponential of the matrix A is evaluated from a diagonal Pade
  approximant. This routine is a modification of the subroutine
  PADE, described in reference [1]. The routine implements an
  algorithm which exploits the identity

      (exp[(2**-m)*A]) ** (2**m) = exp(A),

  where m is an integer determined by the algorithm, to improve the
  accuracy for matrices with large norms.

References
  [1] Ward, R.C.
      Numerical computation of the matrix exponential with accuracy
      estimate.
      SIAM J. Numer. Anal., 14, pp. 600-610, 1977.

Numerical Aspects
                            3
  The algorithm requires 0(N ) operations.

Further Comments
  None
Example

Program Text

*     MB05OD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX
      PARAMETER        ( NMAX = 20 )
      INTEGER          LDA
      PARAMETER        ( LDA = NMAX )
      INTEGER          NDIAG
      PARAMETER        ( NDIAG = 9 )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = NMAX*( 2*NMAX+NDIAG+1 )+NDIAG )
*     .. Local Scalars ..
      DOUBLE PRECISION DELTA
      INTEGER          I, IDIG, INFO, IWARN, J, MDIG, N
      CHARACTER*1      BALANC
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK)
      INTEGER          IWORK(NMAX)
*     .. External Subroutines ..
      EXTERNAL         MB05OD
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, DELTA, BALANC
      IF ( N.LE.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99994 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
*        Find the exponential of the real defective matrix A*DELTA.
         CALL MB05OD( BALANC, N, NDIAG, DELTA, A, LDA, MDIG, IDIG,
     $                IWORK, DWORK, LDWORK, IWARN, INFO )
*
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            IF ( IWARN.NE.0 )
     $         WRITE ( NOUT, FMT = 99993 ) IWARN
            WRITE ( NOUT, FMT = 99997 )
            DO 20 I = 1, N
               WRITE ( NOUT, FMT = 99996 ) ( A(I,J), J = 1,N )
   20       CONTINUE
            WRITE ( NOUT, FMT = 99995 ) MDIG, IDIG
         END IF
      END IF
      STOP
*
99999 FORMAT (' MB05OD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB05OD = ',I2)
99997 FORMAT (' The solution matrix E = exp(A*DELTA) is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' Minimal number of accurate digits in the norm of E =',
     $       I4,/' Number of accurate digits in the norm of E',/'     ',
     $       '            at 95 per cent confidence interval =',I4)
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (' IWARN on exit from MB05OD = ',I2)
      END
Program Data
 MB05OD EXAMPLE PROGRAM DATA
   3     1.0     S
   2.0   1.0   1.0
   0.0   3.0   2.0
   1.0   0.0   4.0
Program Results
 MB05OD EXAMPLE PROGRAM RESULTS

 The solution matrix E = exp(A*DELTA) is 
  22.5984  17.2073  53.8144
  24.4047  27.6033  83.2241
  29.4097  12.2024  81.4177

 Minimal number of accurate digits in the norm of E =  12
 Number of accurate digits in the norm of E
                 at 95 per cent confidence interval =  15

Return to index