In cubic spline interpolation, we use the set of knots and function values $(x_i,y_i),i=1,...,n$ to construct a (tridiagonal) system of equations for the unknowns $\sigma_i$:
$$ h_{i-1}\sigma_{i-1} + 2(h_{i-1} + h_i)\sigma_i + h_i \sigma_{i+1} = \Delta_i - \Delta_{i-1} $$
where $h_i = x_{i+1} - x_i$ are the distances between knots, and $\Delta_i = (y_{i+1} - y_i)/(x_{i+1} - x_i)$ are approximations to the first derivatives. Upon solving the system of equations (together with suitable boundary equations), the values $\sigma_i = s''(x_i)/6$ (related to the second derivative of the cubic spline interpolant), together with the function values $y_i$ are used to form a piecewise cubic interpolant of the form: $$ s(x)=w y_{i+1} + \bar{w}y_i + h_i^2 \left[(w^3 - w)\sigma_{i+1} + (\bar{w}^3 - \bar{w})\sigma_i\right] $$ and $w = (x - x_i)/h_i$, $\bar{w} = 1-w$. For any choice of $x_1 < x_2 < ... < x_n$ the matrix of this system of equations will be nonsingular and diagonally dominant, meaning we can calculate a solution using Gaussian elimination without scaling or pivoting, and still obtain accurate results.
For general tridiagonal matrices LAPACK offers the routine dgtsv using LU factorization with pivoting. Intel MKL also offers the routine dtsvb for which does not perform pivoting (and therefore increases performance).
For certain types of boundary conditions (e.g. natural, or not-a-knot), on top of the diagonal dominance, the system of equations will also be symmetric. A specialized solver valid for symmetric, tridiagonal, and positive-definite systems is also available in LAPACK under dptsv (it computes a $LDL^T$ factorization). My primary source,
Forsythe, Malcolm, and Moler. Computer methods for Mathematical Computations. Prentice-Hall, 1977. (pg. 70-76)
does not indicate if the cubic spline system of equations is positive definite or not. I have found some references which state that for the B-spline basis, the collocation matrix is positive definite (see Novaković, 2018), however this appears to be a fundamentally different spline representation.
Is the positive-definite solver dptsv suitable for cubic spline interpolation in the form given above, or should I stick to the non-symmetric routines based upon LU factorization?
Edit: I have done the test suggested by Charlie S. It can be found in my gist.
The "not-a-knot" spline matrix is not positive definite: $$ \begin{pmatrix} -h_1 & h_1 & & & & 0 \\ h_1 & 2(h_1 + h_2) & h_2 \\ 0 & h_2 & 2(h_2 + h_3) & h_3 \\ & & \ddots & \ddots & \ddots \\ & & & h_{n-2} & 2(h_{n-2} + h_{n-1}) & h_{n-1} \\ 0 & & & & h_{n-1} & -h_{n-1} \end{pmatrix} \begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \vdots \\ \sigma_{n-1} \\ \sigma_n \end{pmatrix} = \begin{pmatrix} h_1^2\Delta_1^{(3)}\\ \Delta_2 - \Delta_1 \\ \Delta_3 - \Delta_2 \\ \vdots \\ \Delta_{n-1} - \Delta_{n-2} \\ -h_{n-1}^2 \Delta_{n-3}^{(3)} \end{pmatrix} $$
The spline with natural boundary conditions $s''(x) = 0$ at the ends, meaning $\sigma_1 = \sigma_n = 0$, is positive definite (in agreement with the answer of Federico Poloni):
$$ \begin{pmatrix} 2(h_1 + h_2) & h_2 & & & 0\\ h_2 & 2(h_2 + h_3) & h_3 \\ & \ddots & \ddots & \ddots \\ & & h_{n-3} & 2(h_{n-3} + h_{n-2}) & h_{n-2} \\ 0 & & & h_{n-2} & 2(h_{n-2} + h_{n-1}) \end{pmatrix} \begin{pmatrix} \sigma_2 \\ \sigma_3 \\ \vdots \\ \sigma_{n-2} \\ \sigma_{n-1} \end{pmatrix} = \begin{pmatrix} \Delta_2 - \Delta_1 \\ \Delta_3 - \Delta_2 \\ \vdots \\ \Delta_{n-2} - \Delta_{n-3} \\ \Delta_{n-1} - \Delta_{n-2} \end{pmatrix} $$
For the natural cubic spline it is possible to use the symmetric routine dptsv.
Based on the answer of Federico Poloni, I suspect the periodic cubic spline will have a positive definite matrix too (but not tridiagonal).
