`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.