Tridiagonal linear system of equations example using gttrf from LAPACK (Intel MKL)

test_gttrf.f90

! LINK="${MKLROOT}/lib/intel64/libmkl_lapack95_lp64.a -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a ${MKLROOT}/lib/intel64/libmkl_sequential.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl"
! INC="-I${MKLROOT}/include/intel64/lp64/ -I${MKLROOT}/include/"
! ifort -warn all -o test_gttrf $INC test_gttrf.f90 $LINK 

! ./test_gttrf > test_gttrf.out
! gnuplot
! gnuplot> Th = 1
! gnuplot> plot "test_gttrf.out" u 1:2 w p pt 7, cosh((1-x)*Th)/cosh(Th)

module rd_setup

    use lapack95, only: gttrf,gttrs
    use f95_precision, only: wp => dp
    implicit none
    private

    public :: setup, Thiele, wp

    real(wp), parameter :: length = 1.0_wp
    real(wp), parameter :: diff = 1.0_wp
    real(wp), parameter :: reac = 1.0_wp
    real(wp), parameter :: Thiele = length*sqrt(reac/diff)

contains

    subroutine setup(n,dl,d,du,b,x)
        integer, intent(in) :: n
        real(wp), intent(out) :: dl(n-1), d(n), du(n-1), b(n)
        real(wp), intent(out) :: x(n)

        real(wp) :: h, twodiff,rfactor
        integer :: i

        h = length/real(n-1,wp) ! step size
        x = [((i-1)*h,i=1,n)]   ! x coordinates

        twodiff = 2._wp*diff
        rfactor = h**2*reac

        d(1) = 1.0_wp; du(1) = 0.0_wp ! Dirichlet boundary conditon
        do i = 2, n-1
            dl(i-1) = diff             ! D
            d(i) =  -twodiff - rfactor ! -2*D - h**2*k
            du(i) = diff               ! D
            b(i) = 0._wp
        end do
        d(n) = -twodiff - rfactor; dl(n-1) = twodiff ! zero-flux (ghost-node trick)


        b(1) = 1.0_wp ! left Dirichlet boundary condition
        b(n) = 0.0_wp ! right zero-flux boundary condition

    end subroutine

end module

program test_tridiagonal

    use lapack95, only: gttrf, gttrs ! factorize and solve general system with tridiagonal coefficient matrix
    use rd_setup, only: wp, setup, Thiele
    implicit none
    
    integer, parameter :: n = 11
    real(wp), allocatable :: dl(:), d(:), du(:), du2(:), b(:), x(:)
    integer, allocatable :: ipiv(:)
    integer :: info, i

    allocate(dl(n-1),d(n),du(n-1),du2(n-2),ipiv(n),b(n),x(n))

    ! Populate tridiagonal matrix
    call setup(n,dl,d,du,b,x)
    print *, "# Thiele = ", Thiele

    ! Factorize matrix
    call gttrf(dl,d,du,du2,ipiv,info=info)
    print *, "# factorize info = ", info
    
    ! Solve tridiagonal matrix
    call gttrs(dl,d,du,du2,b,ipiv,info=info) ! trans = 'N'
    print *, "# solve info = ", info

    ! Print solution
    do i = 1, n
        print *, x(i), b(i)
    end do
end program