! Compute American Option Premium for Vanilla Put
      use feynman_kac_int
      use hqsval_int
      use mp_types
      use umach_int
      implicit none
! The strike price
      real(kind(1e0)) :: ks = 10.0e0
! The sigma value
      real(kind(1e0)) :: sigma = 0.4e0

! Time values for the options
      integer, parameter :: nt = 2
      real(kind(1e0)) :: time(nt)=(/0.25e0, 0.5e0/)
! Values of the underlying where evaluations are made
      integer, parameter :: nv = 9
      integer, parameter :: nlbc = 2, nrbc = 3, ndeg = 6

      integer :: i
      real(kind(1e0)) :: xs(nv) = (/((i-1)*2.0e0,i=1,nv)/)
! Value of the interest rate and continuous dividend
      real(kind(1e0)) :: r = 0.1e0, dividend = 0.0e0
! Values of the min and max underlying values modeled
      real(kind(1e0)) :: x_min = 0.0e0, x_max = 30.0e0
      real(kind(1e0)) :: atol(1), rtol(1)
 
! Define parameters for the integration step.  
      integer, parameter :: nx = 121, nint = nx-1, n = 3*nx
      real(kind(1e0)) :: xgrid(nx), ye(n,0:nt), yeprime(n,0:nt),&
                         ya(n,0:nt), yaprime(n,0:nt),&
                         dx, fe(nv,nt), fa(nv,nt)
      type(s_fcn_data) fcn_data
      integer :: nout
      real(kind(1e0)), external :: fkcoef_put, fkinitcond_put
      external fkbc_put, fkinit_put, fkforce_put
      
      call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
      allocate(fcn_data % rdata (6), fcn_data % idata (1))
! Define an equally-spaced grid of points for the underlying price
      dx = (x_max-x_min)/real(nint)
      xgrid(1) = x_min
      xgrid(nx) = x_max
      do i=2,nx-1
        xgrid(i) = xgrid(i-1) + dx
      end do

! Place a breakpoint at the strike price.
      do i = 1,nx
        if (xgrid(i) > ks) then
          xgrid(i-1) = ks
          exit
        end if
      end do
! Request less accuracy than the default values provide.
      atol(1) = 0.5e-2
      rtol(1) = 0.5e-2
      fcn_data % rdata = (/ks,x_max,sigma,r,dividend,atol(1)/)
      fcn_data % idata = (/ndeg/)

! Compute European then American Put Option Values.
      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&
                        fkinitcond_put, fkbc_put, ye, yeprime,&
                        FKINIT=fkinit_put, ATOL=atol,RTOL=rtol,&
                        FCN_DATA = fcn_data)
      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&
                        fkinitcond_put, fkbc_put, ya, yaprime,&
                        FKINIT=fkinit_put, ATOL=atol, RTOL=rtol,&
                        FKFORCE=fkforce_put, FCN_DATA = fcn_data)

! Evaluate solutions at vector of points XS(:), at each time value
! prior to expiration.                        
      do i=1,nt
        fe(:,i) = hqsval (xs, xgrid, ye(:,I))
        fa(:,I) = hqsval (xs, xgrid, ya(:,I))
      end do 
      write(nout,'(T05,A,/,T05,A)')&
         "American Option Premium for Vanilla Put, 3 and 6 Months "//&
         "Prior to", "Expiry"
      write(nout,'(T08,"Number of equally spaced spline knots ",I4,'//&
                  '/T08,"Number of unknowns ",I4)') nx,n          
      write(nout,'(T08,"Strike= ",F5.2,", Sigma=", F5.2,", Interest'//&
               ' Rate=",F5.2,/T08,"Underlying", T26,"European",'//&
               'T42,"American",/(T10,5F8.4))') ks,sigma,r,&
               (xs(i), fe(i,1:nt), fa(i,1:nt),i=1,nv)
      end
! These routines define the coefficients, payoff, boundary 
! conditions, forcing term and initial conditions for American and
! European Options.
      function fkcoef_put(x, tx, iflag, fcn_data)
        use mp_types
        implicit none
        integer, intent(inout) :: iflag
        real(kind(1e0)), intent(in) :: x, tx
        type(s_fcn_data), optional :: fcn_data
        real(kind(1e0)) :: fkcoef_put

        real(kind(1e0)) :: sigma, strike_price, interest_rate, &
                           dividend, zero=0.e0
        sigma = fcn_data % rdata(3)
        interest_rate = fcn_data % rdata(4)
        dividend = fcn_data % rdata(5)
        select case (iflag)
          case (1)
! The coefficient derivative d(sigma)/dx
            fkcoef_put = sigma
! The coefficient sigma(x)
          case (2)
            fkcoef_put = sigma*x 
          case (3)
! The coefficient mu(x)
            fkcoef_put = (interest_rate - dividend)*x
          case (4)
! The coefficient kappa(x)
            fkcoef_put = interest_rate
        end select
! Note that there is no time dependence		
        iflag = 0
        return
      end function fkcoef_put

      function fkinitcond_put(x, fcn_data)
        use mp_types
        implicit none
        real(kind(1e0)), intent(in) :: x
        type (s_fcn_data), optional :: fcn_data
        real(kind(1e0)) :: fkinitcond_put
        real(kind(1e0)) :: zero = 0.0e0
        real(kind(1e0)) :: strike_price

        strike_price = fcn_data % rdata(1)
! The payoff function			 
        fkinitcond_put = max(strike_price - x, zero)
        return
      end function fkinitcond_put

      subroutine fkbc_put (tx, iflag, bccoefs, fcn_data)
        use mp_types
        implicit none
        real(kind(1e0)), intent(in) :: tx
        integer, intent(inout) :: iflag
        real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
        type (s_fcn_data), optional :: fcn_data

        select case (iflag)
          case (1)       
            bccoefs(1,1:4) = ((/0.0e0, 1.0e0, 0.0e0, -1.0e0/))
            bccoefs(2,1:4) = ((/0.0e0, 0.0e0, 1.0e0,  0.0e0/))
          case (2)
            bccoefs(1,1:4) = ((/1.0e0, 0.0e0, 0.0e0, 0.0e0/))
            bccoefs(2,1:4) = ((/0.0e0, 1.0e0, 0.0e0, 0.0e0/))
            bccoefs(3,1:4) = ((/0.0e0, 0.0e0, 1.0e0, 0.0e0/))
        end select
! Note no time dependence
        iflag = 0		
      end subroutine fkbc_put

      subroutine fkforce_put (interval, t, hx, y, xlocal, qw, u,&
                              phi, dphi, fcn_data)
        use mp_types
        implicit none
        integer, parameter :: local = 6
        integer :: i, j, l, ndeg
        integer, intent(in) :: interval
        real(kind(1e0)), intent(in) :: y(:), t, hx, qw(:),&
                                       xlocal(:), u(:,:)
        real(kind(1e0)), intent(out) :: phi(:), dphi(:,:)
        type (s_fcn_data), optional :: fcn_data

        real(kind(1e0)) :: yl(local), bf(local)
        real(kind(1e0)) :: value, strike_price, interest_rate,&
                           zero=0.0e0, one=1.0e0, rt, mu

        yl = y(3*interval-2:3*interval+3)
        phi = zero

        value = fcn_data % rdata(6)
        strike_price = fcn_data % rdata(1)
        interest_rate = fcn_data % rdata(4)
        ndeg = fcn_data % idata(1)

        mu = 2
! This is the local definition of the forcing term
        do j=1,local
          do l=1,ndeg
            bf(1:3) = u(l,1:3)
            bf(4:6) = u(l,7:9)
            rt = dot_product(yl,bf)
            rt = value/(rt + value - (strike_price - xlocal(l)))
            phi(j) = phi(j) + qw(l) * bf(j) * rt**mu
          end do
        end do
        phi = -phi*hx*interest_rate*strike_price

! This is the local derivative matrix for the forcing term

        dphi = zero
        do j =1,local
          do i = 1,local
            do l=1,ndeg
              bf(1:3) = u(l,1:3)
              bf(4:6) = u(l,7:9)
              rt = dot_product(yl,bf)
              rt = one/(rt + value - (strike_price - xlocal(l)))
              dphi(i,j) = dphi(i,j) + qw(l) * bf(I) * bf(j) *&
                          rt**(mu+1)
            end do
          end do
        end do
        dphi = mu*dphi*hx*value**mu*interest_rate*strike_price
        return
      end subroutine fkforce_put

      subroutine fkinit_put(xgrid,tgrid,t,yprime,y,atol,rtol,&
                            fcn_data)
        use mp_types
        implicit none
        real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), t,&
                                       yprime(:)
        real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)
        type (s_fcn_data), optional :: fcn_data
        integer :: i

        if (t == 0.0e0) then
! Set initial data precisely.  The strike price is a breakpoint.
! Average the derivative limit values from either side.
          do i=1,size(xgrid)
            if (xgrid(i) < fcn_data % rdata(1)) then
              y(3*i-2) = fcn_data % rdata(1) - xgrid(i)
              y(3*i-1) = -1.0e0
              y(3*i)= 0.0e0
            else if (xgrid(i) ==  fcn_data % rdata(1)) then
              y(3*i-2) = 0.0e0
              y(3*i-1) = -0.5e0
              y(3*i) = 0.0e0
            else
              y(3*i-2) = 0.0e0
              y(3*i-1) = 0.0e0
              y(3*i) = 0.0e0
            end if
          end do
        end if
      end subroutine fkinit_put
