program matrixexp
c      Matrix exponential
c      ************************************************************
c      n: Order of the matrix
c      tdelta: time step (decrease tdelta for higher accuracy)
c      tend: time $T$ upto which time stepping is carried out
c      A(i,j): Input matrix (a(i,j) in the program)
c      program finds Exp(AT)
c      Redirect output to file: Example  matrixexp.e >> matrixexp.out
c      ************************************************************
       integer i,j,k,m,n,nsteps,ierr
       real*8 tdelta,tend,a(100,100),stiff(100,100),b(100,100),
     & asq(100,100),ident(100,100),zn(100,100),znp1(100,100)
       open(unit=5,file='matrixexp.dat',status='old')
       read(5,*) n,tdelta,tend
       do 10 i=1,n
         read(5,*) (a(i,j),j=1,n)
10     continue
       nsteps=nint(tend/tdelta)
       ident=0.0d0
       asq=0.0d0
       do 20 i=1,n
        ident(i,i)=1.0d0
        zn(i,i)=1.0d0
20     continue
       do 40 i=1,n
        do 30 j=1,n
         asq(i,j)=0.0d0
         do 25 k=1,n
          asq(i,j)=asq(i,j)+a(i,k)*a(k,j)
25       continue
30      continue
40     continue
       stiff(1:n,1:n)=6.0d0*(2.0d0*ident(1:n,1:n)-tdelta*a(1:n,1:n))+
     & tdelta**2*asq(1:n,1:n)
       b(1:n,1:n)=6.0d0*(2.0d0*ident(1:n,1:n)+tdelta*a(1:n,1:n))+
     & tdelta**2*asq(1:n,1:n)
      call gaussj(stiff,n,100,b,n,100,ierr)
      do 100 m=1,nsteps
       do 70 i=1,n
        do 60 j=1,n
          znp1(i,j)=0.0d0
         do 50 k=1,n
          znp1(i,j)=znp1(i,j)+b(i,k)*zn(k,j)
50       continue
60      continue
70     continue
       zn(1:n,1:n)=znp1(1:n,1:n)
100   continue
      do 200 i=1,n
       print*, (znp1(i,j),j=1,n)
200   continue
      end
      subroutine gaussj(a,n,np,b,m,mp,ierr)

c  Purpose: Solution of the system of linear equations AX = B by
c     Gauss-Jordan elimination, where A is a matrix of order N and B is
c     an N x M matrix.  On output A is replaced by its matrix inverse
c     and B is preplaced by the corresponding set of solution vectors.

c  Source: W.H. Press et al, "Numerical Recipes," 1989, p. 28.

c  Modifications: 
c     1. Double  precision.
c     2. Error parameter IERR included.  0 = no error. 1 = singular 
c        matrix encountered; no inverse is returned.

c  Prepared by J. Applequist, 8/17/91.

      integer i,j,k,l,ll,m,mp,n,np,nmax,ipiv,indxr,indxc,irow,icol,ierr
      real*8 a,b,big,dum,pivinv

c        Set largest anticipated value of N.

      parameter (nmax=500)
      dimension a(np,np),b(np,mp),ipiv(nmax),indxr(nmax),indxc(nmax)
      ierr=0
      irow=0
      icol=0
      do 11 j=1,n
      ipiv(j)=0
 11   continue
      do 22 i=1,n
      big=0.d0
      do 13 j=1,n
      if (ipiv(j).ne.1) then
      do 12 k=1,n
      if (ipiv(k).eq.0) then
      if (dabs(a(j,k)).ge.big) then
      big=dabs(a(j,k))
      irow=j
      icol=k
      endif
      else if (ipiv(k).gt.1) then
      ierr=1
      return
      endif
 12   continue
      endif
 13   continue
      ipiv(icol)=ipiv(icol)+1
      if (irow.ne.icol) then
      do 14 l=1,n
      dum=a(irow,l)
      a(irow,l)=a(icol,l)
      a(icol,l)=dum
 14   continue
      do 15 l=1,m
      dum=b(irow,l)
      b(irow,l)=b(icol,l)
      b(icol,l)=dum
 15   continue
      endif
      indxr(i)=irow
      indxc(i)=icol
      if (a(icol,icol).eq.0.d0) then
      ierr=1
      return
      endif
      pivinv=1.d0/a(icol,icol)
      a(icol,icol)=1.d0
      do 16 l=1,n
      a(icol,l)=a(icol,l)*pivinv
 16   continue
      do 17 l=1,m
      b(icol,l)=b(icol,l)*pivinv
 17   continue
      do 21 ll=1,n
      if (ll.ne.icol) then
      dum=a(ll,icol)
      a(ll,icol)=0.d0
      do 18 l=1,n
      a(ll,l)=a(ll,l)-a(icol,l)*dum
 18   continue
      do 19 l=1,m
      b(ll,l)=b(ll,l)-b(icol,l)*dum
 19   continue
      endif
 21   continue
 22   continue
      do 24 l=n,1,-1
      if (indxr(l).ne.indxc(l)) then
      do 23 k=1,n
      dum=a(k,indxr(l))
      a(k,indxr(l))=a(k,indxc(l))
      a(k,indxc(l))=dum
 23   continue
      endif
 24   continue
      return
      end