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.