Spline interpolation matrices in FORTRAN

spline.f90

! ----------------------------------------------------------------------
! SPDX-License-Identifier: CC0-1.0
! Authored by Fereydoun Memarzanjany
!
! To the maximum extent possible under law, the Author waives all
! copyright and related or neighboring rights to this code.
!
! You should have received a copy of the CC0 legalcode along with this
! work; if not, see <http://creativecommons.org/publicdomain/zero/1.0/>
! ----------------------------------------------------------------------

program spline_coefficients
    implicit none ! I want types to be explicit, not Hungarian notation.
    integer, parameter :: n = 6 ! Number of data points
    real               :: xj(n), fj(n) ! Points and function values
    real               :: hj(n-1), aj(n), bj(n-1), cj(n), dj(n-1)
    real               :: alpha(n-1), l(n), mu(n), z(n)
    integer            :: i ! Counter for loop indices (1-indexed)

    ! Input data points
    xj = (/ 0.0, 1.0, 2.0, 10.0, 50.0, 100.0 /)
    fj = (/ 100.0, 45.0, 39.0, 22.0, 5.0, 0.05 /)

    ! Compute hj (i.e., heights or distances between intervals)
    do i = 1, n-1
        hj(i) = xj(i+1) - xj(i)
    end do

    ! ------------------------------------------------------------------
    ! Solve the Sparse Matrix
    ! ------------------------------------------------------------------
        ! Compute alpha (temporary variable)
        alpha(1) = 0.0
        do i = 2, n-1
            alpha(i) = (3.0/hj(i)  ) * (fj(i+1) - fj(i)  ) &
                     - (3.0/hj(i-1)) * (fj(i)   - fj(i-1))
        end do

        ! Tridiagonal system initialization
        l(1)  = 1.0 ! Diagonal Element
        mu(1) = 0.0 ! Multipliers (Greek \Mu letter)
        z(1)  = 0.0 ! Temporary holders

        ! Forward elimination
        do i = 2, n-1
            l(i)  = 2.0 * (xj(i+1) - xj(i-1)) - hj(i-1) * mu(i-1)
            mu(i) = hj(i) / l(i)
            z(i)  = (alpha(i) - hj(i-1) * z(i-1)) / l(i)
        end do

        ! Boundary conditions
        l(n) = 1.0
        z(n) = 0.0
        cj(n) = 0.0

        ! Back substitution
        do i = n-1, 1, -1
            cj(i) = z(i) - mu(i) * cj(i+1)
            bj(i) = (fj(i+1) - fj(i)) / hj(i) &
                  - (hj(i)   * (cj(i+1) + 2.0 * cj(i))) / 3.0
            dj(i) = (cj(i+1) - cj(i)) / (3.0 * hj(i))
        end do
        cj(1) = 0.0
    ! ------------------------------------------------------------------

    ! Output the coefficients
    print '(A)', 'j    cj           bj           dj'
    do i = 1, n-1
        print '(I2,2X,F10.6,2X,F10.6,2X,F12.6)', &
                i-1,  cj(i),   bj(i),   dj(i)
    end do
    ! Print the last cj alone.
    print '(I2,2X,F10.6)', &
            n-1,  cj(n)

end program spline_coefficients

! ----------------------------------------------------------------------
! END OF FILE: spline.f90
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Sample Output:
! ----------------------------------------------------------------------
! j    cj           bj           dj
!  0    0.000000  -67.375244     12.375247
!  1   37.125740  -30.249504    -12.876235
!  2   -1.502965    5.373271      0.070710
!  3    0.194078   -5.097822     -0.001931
!  4   -0.037695    1.157507      0.000251
!  5    0.000000
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Sample Output:
! ----------------------------------------------------------------------
! j    cj           bj           dj
!  0    0.000000  -67.375244     12.375247
!  1   37.125740  -30.249504    -12.876235
!  2   -1.502965    5.373271      0.070710
!  3    0.194078   -5.097822     -0.001931
!  4   -0.037695    1.157507      0.000251
!  5    0.000000
! ----------------------------------------------------------------------