program main

c*********************************************************************72
c
cc MAIN is a test program for exponential propagator using Arnoldi approach.
c
c  Discussion:
c
c    This main program is a very simple test using diagonal matrices
c    (Krylov subspace methods are blind to the structure of the matrix
c    except for symmetry). This provides a good way of testing the
c    accuracy of the method as well as the error estimates.
c
c  Modified:
c
c    02 July 2005
c
      implicit none

      integer nmax
      parameter ( nmax = 150 )
      integer ih0
      parameter ( ih0 = 60 )
      integer ndmx
      parameter ( ndmx = 20 )

      integer nzmax
      parameter ( nzmax = 7 * nmax )

      double precision a(nzmax)
      double precision :: a0 = 0.0
      double precision :: b0 = 1.0
      double precision ddot
      double precision eps
      double precision h
      integer ioff(10)
      integer j
      integer k
      integer m
      integer n
      integer ndiag
      double precision t
      double precision tn
      double precision u(ih0*nmax)
      double precision w(nmax)
      double precision w1(nmax)
      double precision x(nmax)
      double precision y(nmax)

      call timestamp ( )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SPARSEKIT_PRB12'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '  Test EXPPRO, which computes the matrix exponential.'
c
c  Set dimension of matrix.
c
      n = 100
c
c  Define matrix.
c  A is a single diagonal matrix (ndiag = 1 and ioff(1) = 0 )
c
      ndiag = 1
      ioff(1) = 0
c
c  entries in the diagonal are uniformly distributed.
c
      h = 1.0 / dble ( n + 1 )
      do j = 1, n
        a(j) = dble ( j + 1 ) / dble ( n + 1 )
      end do
c
c  Set a value for TN
c
      tn = 2.0

      eps = 0.0001

      m = 5
      write ( *, '(a,i6)' ) '  Dimension of Krylov subspace M = ', m
c
c  Define initial conditions: chosen so that solution = (1,1,1,1..1)^T
c
      do j = 1, n
        w(j) = exp ( a(j) * tn )
        w1(j) = w(j)
      end do

      call expprod ( n, m, eps, tn, u, w, x, y, a, ioff, ndiag )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  First 10 components of final answer:'
      write ( *, '(a)' ) ' '
      do k = 1, 10
        write ( *, * ) w(k)
      end do
      write ( *, '(a)' ) ' '

      do k = 1, n
        w1(k) = exp ( -a(k) * tn ) * w1(k)
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  First 10 components of exact solution '
      write ( *, '(a)' ) ' '
      do k = 1, 10
        write ( *, * ) w1(k)
      end do
c
c  Compute actual 2-norm of error.
c
      t = 0.0
      do k = 1, n
        t = t + ( w1(k) - w(k) )**2
      end do
      t = sqrt ( t / ddot ( n, w, 1, w, 1 ) )

      write ( *, '(a)' ) ' '
      write ( *, * ) '  RMS error (approx-exact)=', t

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SPARSEKIT_PRB12'
      write ( *, '(a)' ) '  Normal end of execution.'

      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      function afun ( x, y, z )

c*********************************************************************72
c
cc AFUN
c
      implicit none

      double precision afun
      double precision x
      double precision y
      double precision z

      afun = -1.0D+00

      return
      end
      function bfun ( x, y, z )

c*********************************************************************72
c
cc BFUN
c
      implicit none

      double precision bfun
      double precision x
      double precision y
      double precision z

      bfun = -1.0D+00

      return
      end
      function cfun ( x, y, z )

c*********************************************************************72
c
cc CFUN
c
      implicit none

      double precision cfun
      double precision x
      double precision y
      double precision z

      cfun = -1.0D+00

      return
      end
      function dfun ( x, y, z )

c*********************************************************************72
c
cc DFUN
c
      implicit none

      double precision dfun
      double precision x
      double precision y
      double precision z

      dfun = 0.0D+00

      return
      end
      function efun ( x, y, z )

c*********************************************************************72
c
cc EFUN
c
      implicit none

      double precision efun
      double precision x
      double precision y
      double precision z

      efun = 0.0D+00

      return
      end
      function ffun ( x, y, z )

c*********************************************************************72
c
cc FFUN
c
      implicit none

      double precision ffun
      double precision x
      double precision y
      double precision z

      ffun = 0.0D+00

      return
      end
      function gfun ( x, y, z )

c*********************************************************************72
c
cc GFUN
c
      implicit none

      double precision gfun
      double precision x
      double precision y
      double precision z

      gfun = 0.0D+00

      return
      end
      subroutine afunbl ( nfree, x, y, z, coeff )

c*********************************************************************72
c
cc AFUNBL ???
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
        coeff((j-1)*nfree+j) = -1.0
      end do

      return
      end
      subroutine bfunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
        coeff((j-1)*nfree+j) = -1.0
      end do

      return
      end
      subroutine cfunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
        coeff((j-1)*nfree+j) = -1.0
      end do

      return
      end
      subroutine dfunbl (nfree,x,y,z,coeff)
 
c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end
      subroutine efunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end
      subroutine ffunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end
      subroutine gfunbl (nfree,x,y,z,coeff)

c*********************************************************************72
c
      implicit none

      integer i
      integer j
      integer nfree
      double precision x, y, z, coeff(100)

      do j=1, nfree
        do i=1, nfree
          coeff((j-1)*nfree+i) = 0.0
        end do
      end do

      return
      end