**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.

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(*)

**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).

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.

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.

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.

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.

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.

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

3 The algorithm requires 0(N ) operations.

None

**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

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

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