! Compute value of a Convertible Bond
      use feynman_kac_int
      use hqsval_int
      use mp_types
      use umach_int

      implicit none
! The face value
      real(kind(1e0)) :: ks = 1.0e0
! The sigma  or volatility value
      real(kind(1e0)) :: sigma = 0.25e0

! Time values for the options
      integer, parameter :: nt = 2
      real(kind(1e0)) :: time(nt)=(/0.5e0, 1.0e0/)
! Values of the underlying where evaluation are made
      integer, parameter :: nv = 13
      integer, parameter :: nlbc = 3, nrbc = 3, ndeg = 6
      integer :: i
      real(kind(1e0)) :: xs(nv) = (/((i-1)*0.25e0,i=1,nv)/)
! Value of the interest rate, continuous dividend and factor
      real(kind(1e0)) :: r = 0.1e0, dividend = 0.02e0,&
                         factor =1.125e0
! Values of the min and max underlying values modeled
      real(kind(1e0)) :: x_min = 0.0e0, x_max = 4.0e0
! Define parameters for the integration step.  
      integer, parameter :: nx = 61, nint = nx-1, n = 3*nx
      real(kind(1e0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt),&
                         dx, f(nv,0:nt)
! Relative and absolute error tolerances
      real(kind(1e0)) :: atol(1), rtol(1)
      type(s_fcn_data) fcn_data
      real(kind(1e0)), external :: fkcoef_cbond, fkinitcond_cbond
      external fkbc_cbond, fkforce_cbond, fkinit_cbond

      integer :: nout

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

! Use a pure absolute error tolerance for the integration
! The default values require too much integation time.
      atol(1) = 1.0e-3
      rtol(1) = 0.0e0

! Pass the data for evaluation  
      fcn_data % rdata = (/ks,x_max,sigma,r,dividend,factor,&
                           atol(1)/)
      fcn_data % idata = (/ndeg/)
          
! Compute value of convertible bond                 
      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_cbond,&
                     fkinitcond_cbond, fkbc_cbond, y, yprime,&
                     ATOL=atol, RTOL=rtol, FKINIT = fkinit_cbond,&
                     FKFORCE = fkforce_cbond, FCN_DATA = fcn_data)
   
! Evaluate and display solutions at vector of points XS(:), at each
! time value prior to expiration.                        
      do i=0,nt
        f(:,i) = hqsval (xs, xgrid, y(:,i))
      end do

      write(nout,'(T05,A)')&
      "Convertible Bond Value, 0+, 6 and 12 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,/'//&
         'T08,"Interest Rate =",F5.2,", Dividend =",F5.2,'//&
         '", Factor = ",F5.3,//T08,"Underlying", T26,"Bond Value",'//&
         '/(T10,4F8.4))') ks,sigma,r,dividend,factor,&
           (xs(i), f(i,0:nt),i=1,nv)
      end

! These routines define the coefficients, payoff, boundary 
! conditions and forcing term.
      function fkcoef_cbond(x, tx, iflag, fcn_data) 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
      real(kind(1e0)) :: value

      real(kind(1e0)) :: sigma, 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			
          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_cbond

      function fkinitcond_cbond(x, fcn_data) result(value)
      use mp_types
      implicit none
      real(kind(1e0)), intent(in) :: x
      type (s_fcn_data), optional :: fcn_data
      real(kind(1e0)) :: value

      real(kind(1e0)) :: strike_price, factor

      strike_price = fcn_data % rdata(1)
      factor = fcn_data % rdata(6)
      value = max(factor * x, strike_price)
      return
      end function fkinitcond_cbond

      subroutine fkbc_cbond (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

      real(kind(1e0)) :: interest_rate, strike_price, dp,&
                         factor, x_max

      select case (iflag)
        case (1)
          strike_price = fcn_data % rdata(1)
          interest_rate = fcn_data % rdata(4)
          dp = strike_price * exp(tx*interest_rate)
          bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0,    dp/)
          bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0,  0.0e0/)
          bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0,  0.0e0/)
          return
        case (2)
          x_max = fcn_data % rdata(2)
          factor = fcn_data % rdata(6)
          bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, factor * x_max/)
          bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, factor/)
          bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
      end select
! Note no time dependence                  
      iflag = 0
      return
      end subroutine fkbc_cbond

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

      integer :: ndeg
      real(kind(1.e0)) :: yl(local), bf(local)
      real(kind(1.e0)) :: value, strike_price, interest_rate,&
                          zero = 0.0e0, one = 1.0e0, rt, mu, factor
      type(s_fcn_data), optional :: fcn_data

      yl = y(3*interval-2:3*interval+3)
      phi = zero
      dphi = zero
      value = fcn_data % rdata(7)
      strike_price = fcn_data % rdata(1)
      interest_rate = fcn_data % rdata(4)
      factor = fcn_data % rdata(6)
      ndeg = fcn_data % idata(1)
      mu = 2
! This is the local definition of the forcing term
! It "forces" the constraint f >= factor*x.        
      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 - factor * xlocal(l))
          phi(j) = phi(j) + qw(l) * bf(j) * rt**mu
        end do
      end do
      phi = -phi * hx * factor * strike_price
! This is the local derivative matrix for the forcing term -        
      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 - factor * xlocal(l))
            dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j)&
                     * (value * rt)**mu * rt
          end do
        end do
      end do
      dphi = -mu * dphi * hx * factor * strike_price
      return                 
      end subroutine fkforce_cbond
    	
      subroutine fkinit_cbond(xgrid,tgrid,t,yprime,y,atol,rtol,&
                              fcn_data)
      use mp_types
      implicit none
      real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)
      real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), yprime(:),&
                                     t
      type(s_fcn_data), optional :: fcn_data

      integer :: i
      if (t == 0.0e0) then
! Set initial data precisely.
        do i=1,size(Xgrid)
          if (xgrid(i)*fcn_data % rdata(6) <&
              fcn_data % rdata(1)) then
            y(3*i-2) = fcn_data % rdata(1)
            y(3*i-1) = 0.0e0
            y(3*i  ) = 0.0e0
          else
            y(3*i-2) = xgrid(i) * fcn_data % rdata(6)
            y(3*i-1) = fcn_data % rdata(6)
            y(3*i  ) = 0.0e0
          end if
        end do
      end if
      end subroutine fkinit_cbond
