! Compute Bermudan-Style Option Premium for Vanilla Put
      use feynman_kac_int
      use hqsval_int
      use mp_types
      use umach_int

      implicit none
      integer :: nout
! The strike price
      real(kind(1e0)) :: ks = 10.0e0
! The sigma value
      real(kind(1e0)) :: sigma = 0.4e0
! Program working stores
      real(kind(1e0)) :: week
! Time values for the options
      integer, parameter :: nt = 26
      integer, parameter :: ndeg = 6
      real(kind(1e0)) :: time(nt), time_end = 0.5e0
! Values of the underlying where evaluation are made
      integer, parameter :: nv = 9, nlbc = 2, nrbc = 3
      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

! Define parameters for the integration step.  
      integer, parameter :: nx = 61, nint = nx-1, n = 3*nx
      real(kind(1e0)) :: xgrid(nx), yb(n,0:nt), ybprime(n,0:nt),&
                         ya(n,0:nt), yaprime(n,0:nt),&
                         ytemp(n,0:1), ytempprime(n,0:1),&
                         dx, fb(nv,nt), fa(nv,nt)
      real(kind(1e0)) :: atol
      type(s_fcn_data) fcn_data_amer, fcn_data_berm
      real(kind(1e0)), external :: fkcoef_put, fkinitcond_put
      external fkbc_put, fkforce_put, fkinit_amer_put, fkinit_berm_put

      call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
      allocate(fcn_data_amer % rdata (6), fcn_data_amer % 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

! Compute time values where American option is computed
      week = time_end/real(nt,kind(week))
      time(1) = week
      do i=2,nt-1
        time(i) = time(i-1) + week
      end do
      time(nt) = time_end

      atol = 1.0e-3
      fcn_data_amer % rdata = (/ks,x_max,sigma,r,dividend,atol/)
      fcn_data_amer % idata = (/ndeg/)

! Compute American Put Option Values at the weekly grid.
      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&
                        fkinitcond_put, fkbc_put, ya, yaprime,&
                        FKINIT = fkinit_amer_put,&
                        FKFORCE = fkforce_put,&
                        FCN_DATA = fcn_data_amer)
! Integrate once again over the weekly grid, using the American
! Option values as initial data for a piece-wise European option
!integration.

! Allocate space to hold coefficient data and initial values.
      allocate(fcn_data_berm % rdata(5+n))
      fcn_data_berm % rdata(1:5) = fcn_data_amer % rdata(1:5)
! Copy initial data so the payoff value is the same for
! American and Bermudan option values.        
      yb(1:n,0) = ya(1:n,0)
      ybprime(1:n,0) = ya(1:n,0)

      do i=0,nt-1
! Move American Option values into place as initial conditions,
! but now integrating with European style over each period of
! the weekly grid.       
        fcn_data_berm % rdata(6:) = ya(1:n,i)
        if (i .eq. 0) then
          call feynman_kac (xgrid, (/time(1)/), nlbc, nrbc,&
                            fkcoef_put, fkinitcond_put, fkbc_put,&
                            ytemp(:,0:1), ytempprime(:,0:1),&
                            FKINIT = fkinit_berm_put,&
                            FCN_DATA = fcn_data_berm)
        else
          call feynman_kac (xgrid, (/time(i+1)-time(i)/),&
                            nlbc, nrbc, fkcoef_put,&
                            fkinitcond_put, fkbc_put,&
                            ytemp(:,0:1), ytempprime(:,0:1),&
                            FKINIT = fkinit_berm_put,&
                            FCN_DATA = fcn_data_berm)
        end if
! Record values of the Bermudan option at the end of each integration.
        yb(1:n,i+1) = ytemp(1:n,1)
        ybprime(1:n,i+1) = ytempprime(1:n,1) 
      end do                                            
! Evaluate solutions at vector of points XS(:), at each time value
! prior to expiration.                        
      do i=1,nt
        fa(:,i) = hqsval (xs, xgrid, ya(:,i))
        fb(:,i) = hqsval (xs, xgrid, yb(:,i))
      end do

      write(nout,'(T05,A)')&
        "American Option Premium for Vanilla Put, 6 Months "//&
        "Prior to Expiry"
      write(nout,'(T05,A)')&
        "Exercise Opportunities At Weekly Intervals"
      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",'//&
                   'T20,"Bermudan Style", T42,"American",'//&
                   '/(T10,F8.4, T26, F8.4, T42, F8.4))')&
                   KS,SIGMA,R,&
                  (xs(i), fb(i,nt:nt), fa(i,nt:nt),i=1,nv)
      end

! These subprograms set the coefficients, payoff, boundary
! conditions and forcing term for American and European Options.
      function fkcoef_put(x, tx, iflag, fcn_data_amer)&
                          result(value)
      use mp_types
      implicit none
      integer, intent(inout) :: iflag
      real(kind(1e0)), intent(in) :: x, tx
      type(s_fcn_data), optional :: fcn_data_amer
      real(kind(1e0)) :: value

      real(kind(1e0)) :: sigma, interest_rate, dividend, zero=0.0e0

      sigma = fcn_data_amer % rdata(3)
      interest_rate = fcn_data_amer % rdata(4)
      dividend = fcn_data_amer % rdata(5)
      select case (iflag)
        case (1)
! The coefficient derivative d(sigma)/dx
          value = sigma
! The coefficient sigma(x)
        case (2)
          value = sigma * x 
        case (3)
! The coefficient mu(x)
          value = (interest_rate - dividend) * x
        case (4)
! The coefficient kappa(x)
          value = interest_rate
      end select
! Note that there is no time dependence
      iflag = 0
      return
      end function fkcoef_put

      function fkinitcond_put(x, fcn_data_amer) result(value)
      use mp_types
      implicit none

      real(kind(1e0)), intent(in) :: x
      type (s_fcn_data), optional :: fcn_data_amer

      real(kind(1e0)) :: value
      real(kind(1e0)) :: strike_price, zero = 0.0e0

      strike_price = fcn_data_amer % rdata(1)
! The payoff function
      value = 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_amer)
      use mp_types
      implicit none
      integer, parameter :: local = 6
      integer :: i, j, l, ndeg
      integer, intent(in) :: interval
      real(kind(1.e0)), intent(in) :: y(:), t, hx, qw(:),&
                                      xlocal(:), u(:,:)
      real(kind(1.e0)), intent(out) :: phi(:), dphi(:,:)
      type (s_fcn_data), optional :: fcn_data_amer

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

      yl = y(3*interval-2:3*interval+3)
      phi = zero
      value = fcn_data_amer % rdata(6)
      strike_price = fcn_data_amer % rdata(1)
      interest_rate = fcn_data_amer % rdata(4)
      ndeg = fcn_data_amer % 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
      end subroutine fkforce_put

      subroutine fkinit_amer_put(xgrid,tgrid,t,yprime,y,atol,rtol,&
                                 fcn_data_amer)
      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_amer
      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_amer % rdata(1)) then
            y(3*i-2) = fcn_data_amer % rdata(1) - xgrid(i)
            y(3*i-1) = -1.0e0
            y(3*i) = 0.0e0
          else if (xgrid(i) == fcn_data_amer % 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_amer_put

      subroutine fkinit_berm_put(xgrid,tgrid,t,yprime,y,atol,rtol,&
                                 fcn_data_berm)
      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_berm
      integer :: i

      if (t == 0.0e0) then
! Set initial data for each week at the previously computed
! American Option values.  These coefficients are passed 
! in the derived type fcn_data_berm.
        do i=1,size(xgrid)
          y(3*i-2) = fcn_data_berm % rdata(3+3*i)
          y(3*i-1) = fcn_data_berm % rdata(4+3*i)
          y(3*i  ) = fcn_data_berm % rdata(5+3*i)
        end do
      end if
      end subroutine fkinit_berm_put
