Output.

• v the (k)-dimensional vector of spline interpolant values, v() = S(x()).

The interface of this routine should be

  subroutine spline(t,y,x,v)
    implicit none
    double precision, intent(in),  dimension(:) :: t,y,x
    double precision, intent(out), dimension(:) :: v
  end subroutine spline

You will need to compute the parameters first. Use natural splines. For the solution of the resulting tridiagonal system. use LAPACK's function dgtsv. Be very careful on the data vectors needed (we do not use the full matrix, only the 3 diagonals as vectors).

Go to http://www.netlib.org/lapack/double/dgtsv.f for more information on how to use this function.

      SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDB, N, NRHS
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * )
*     ..
*
*  Purpose
*  =======
*
*  DGTSV  solves the equation
*
*     A*X = B,
*
*  where A is an n by n tridiagonal matrix, by Gaussian elimination with
*  partial pivoting.
*
*  Note that the equation  A'*X = B  may be solved by interchanging the
*  order of the arguments DU and DL.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix A.  N >= 0.
*
*  NRHS    (input) INTEGER
*          The number of right hand sides, i.e., the number of columns
*          of the matrix B.  NRHS >= 0.
*
*  DL      (input/output) DOUBLE PRECISION array, dimension (N-1)
*          On entry, DL must contain the (n-1) sub-diagonal elements of
*          A.
*
*          On exit, DL is overwritten by the (n-2) elements of the
*          second super-diagonal of the upper triangular matrix U from
*          the LU factorization of A, in DL(1), ..., DL(n-2).
*
*  D       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, D must contain the diagonal elements of A.
*
*          On exit, D is overwritten by the n diagonal elements of U.
*
*  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)
*          On entry, DU must contain the (n-1) super-diagonal elements
*          of A.
*
*          On exit, DU is overwritten by the (n-1) elements of the first
*          super-diagonal of U.
*
*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
*          On entry, the N by NRHS matrix of right hand side matrix B.
*          On exit, if INFO = 0, the N by NRHS solution matrix X.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B.  LDB >= max(1,N).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -i, the i-th argument had an illegal value
*          > 0: if INFO = i, U(i,i) is exactly zero, and the solution
*               has not been computed.  The factorization has not been
*               completed unless i = N.
*

Recall that the compilation is done with

After you found 's you can start to compute the interpolating values v.

• For every x() find the interval where it lies, , and compute the value of the respective polynomial v()= (x());
• If x() then the point is outside left the interpolating interval and we simply return for the value of the interpolating polynomial at that point.
• If x() then the point is outside right the interpolating interval and we simply return for the value of the interpolating polynomial at that point.