Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Example for scicomp.stackexchange thread "Factorization of cubic spline interpolation matrix"
! 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment