program f90_spl
    implicit none
    integer,parameter::NDmax=100
    integer::nd,md,i,io
    real(8)::xdat(1:NDmax),ydat(1:NDmax)
    real(8),allocatable::datax(:),datay(:)
    character(len=50)::fnameR
    character(len=20)::dummy

    call getarg(1,fnameR) !Input data file
    call getarg(2,dummy)
    read(dummy,*) md

    open(11,file=fnameR,status='old')
    do i=1,NDmax
            read(11,*,iostat=io) xdat(i),ydat(i)
            if(io<0)exit
    end do
    close(11)
    nd=i-1
    allocate(datax(1:nd),datay(1:nd))

    do i=1,nd
        datax(i)=xdat(i)
        datay(i)=ydat(i)
    end do
    call SPL3(nd,datax,datay,md)

stop

contains

    subroutine SPL3(nd,datax,datay,md)
        integer,intent(in)::nd,md
        real(8),intent(in)::datax(1:nd),datay(1:nd)

        integer::i,j,nn,mm
        real(8)::a(1:nd),b(1:nd),c(1:nd),d(1:nd)
        real(8)::h0,h1,hj
        real(8)::dx,xx,sj
        real(8),allocatable::ary(:,:),vec(:)

        nn=nd-2
        mm=2
        allocate(ary(1:nn,1:mm),vec(1:nn))
        do i=1,nn
            do j=1,mm
                ary(i,j)=0.0D0
            end do
            vec(i)=0.0D0
        end do
        do i=1,nd
            a(i)=datay(i)
            b(i)=0.0D0
            c(i)=0.0D0
            d(i)=0.0D0
        end do
        do i=2,nd-1
            h0=datax(i)-datax(i-1)
            h1=datax(i+1)-datax(i)
            ary(i-1,1)=2.0D0*(h0+h1)
            ary(i-1,2)=h1
            vec(i-1)=3.0D0*(a(i+1)-a(i))/h1-3.0D0*(a(i)-a(i-1))/h0
        end do

        call CHBAND(nn,mm,ary,vec)

        c(1)=0.0D0
        do i=2,nd-1
            c(i)=vec(i-1)
        end do
        c(nd)=0.0D0
        do i=1,nd-1
            hj=datax(i+1)-datax(i)
            d(i)=(c(i+1)-c(i))/(3.0D0*hj)
            b(i)=(a(i+1)-a(i))/hj-(2.0D0*c(i)+c(i+1))*hj/3.0D0
        end do

        do i=1,nd-1
            dx=(datax(i+1)-datax(i))/dble(md)
            do j=1,md
                xx=dble(j-1)*dx
                sj=a(i)+b(i)*xx+c(i)*xx**2.0D0+d(i)*xx**3.0D0
                write(6,*) datax(i)+xx,sj
            end do
        end do
        i=nd-1
        j=md+1
        xx=dble(j-1)*dx
        sj=a(i)+b(i)*xx+c(i)*xx**2.0D0+d(i)*xx**3.0D0
        write(6,*) datax(i)+xx,sj

    end subroutine SPL3


    subroutine CHBAND(n1,m1,array,vector)
        integer,intent(in)::n1,m1
        real(8),intent(inout)::array(1:n1,1:m1)
        real(8),intent(inout)::vector(1:n1)
        integer :: i,j,k,mm,mn
        real(8) :: z

        array(1,1)=sqrt(array(1,1))
        do j=2,m1
            array(1,j)=array(1,j)/array(1,1)
        end do
        do i=2,n1
            if(n1-i+1>m1)then
                mm=m1
            else
                mm=n1-i+1
            end if
            do j=1,mm
                z=array(i,j)
                if(j/=m1)then
                    if(m1-j+1<i)then
                        mn=m1-j+1
                    else
                        mn=i
                    end if
                    do k=2,mn
                        z=z-array(i-k+1,k)*array(i-k+1,j+k-1)
                    end do
                    if(j>1)then
                        array(i,j)=z/array(i,1)
                    else
                        array(i,1)=sqrt(z)
                    end if
                else
                    array(i,j)=z/array(i,1)
                end if
            end do
        end do

        vector(1)=vector(1)/array(1,1)
        do i=2,n1
            z=vector(i)
            if(i>m1)then
                mm=m1
            else
                mm=i
            end if
            do j=2,mm
                z=z-array(i-j+1,j)*vector(i-j+1)
            end do
            vector(i)=z/array(i,1)
        end do

        vector(n1)=vector(n1)/array(n1,1)
        do i=n1-1,1,-1
            z=vector(i)
            do j=2,m1
                if(i+j-1>n1)exit
                z=z-array(i,j)*vector(i+j-1)
            end do
            vector(i)=z/array(i,1)
        end do
    end subroutine CHBAND

end program f90_spl
