Example for scicomp.stackexchange thread "Factorization of cubic spline interpolation matrix"

spline_test_driver.f90

! For an explanation see:
!   https://scicomp.stackexchange.com/questions/36261/factorization-of-cubic-spline-interpolation-matrix
!
! compile with:
!   gfortran -Wall spline_test_driver.f90 -o spline_test_driver -llapack
!
module spline_test

  implicit none

  integer, parameter :: wp = kind(1.0d0)

  interface
    subroutine dgtsv(n,nrhs,dl,d,du,b,ldb,info)
       INTEGER            INFO, LDB, N, NRHS
       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * )
    end subroutine
    subroutine dptsv(n,nrhs,d,e,b,ldb,info)
       INTEGER            INFO, LDB, N, NRHS
       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * )
    end subroutine
  end interface

contains

  subroutine fill_diagonals_and_rhs(n,x,y,b,d,c)
    integer, intent(in) :: n ! must be larger than 3
    real(wp), intent(in) :: x(n), y(n)
    real(wp), intent(out) :: b(n), c(n), d(n)

    integer :: i

    ! b - diagonal
    ! d - off-diagonal
    ! c - right hand side

    d(1) = x(2) - x(1)
    c(2) = (y(2) - y(1))/d(1)
    
    do i = 2, n - 1
      d(i) = x(i+1) - x(i)
      b(i) = 2.0_wp*(d(i-1) + d(i))
      c(i+1) = (y(i+1) - y(i))/d(i)
      c(i) = c(i+1) - c(i)
    end do

    b(1) = -d(1)
    b(n) = -d(n-1)

    c(1) = c(3)/(x(4) - x(2)) - c(2)/(x(3) - x(1))
    c(n) = c(n-1)/(x(n) - x(n-2)) - c(n-2)/(x(n-1) - x(n-3))

    c(1) = c(1)*d(1)**2/(x(4)-x(1))
    c(n) = -c(n)*d(n-1)**2/(x(n) - x(n-3))

  end subroutine

end module

program spline_test_driver

  use spline_test
  implicit none

  integer, parameter :: n = 51
  real(wp) :: x(n), y(n)
  real(wp) :: b(n), c(n), d(n), dl(n)
  real(wp) :: dx
  integer :: i, info

  dx = 2*4.0_wp*atan(1.0_wp)/real(n-1,wp)

  x = [((i-1)*dx,i=1,n)]
  y = sin(x)

  call fill_diagonals_and_rhs(n,x,y,b,d,c)
  ! copy lower diagonal
  dl = d 

  ! Not-a-knot boundary conditions
  call dgtsv(n,1,dl,b,d,c,n,info)
  print *, "info = ", info
  
  do i = 1, n
    print *, i, c(i)*6.0_wp, -sin(x(i))
  end do

  ! Refill matrix diagonals
  call fill_diagonals_and_rhs(n,x,y,b,d,c)

  ! Natural boundary conditions
  call dptsv(n-2,1,b(2:n-1),d(2:n-1),c(2:n-1),n-2,info)
  print *, "info = ", info
  c(1) = 0.0_wp
  c(n) = 0.0_wp


  do i = 1, n
    print *, i, c(i)*6.0_wp, -sin(x(i))
  end do

end program