Function solve_tridiagonal 

pub fn solve_tridiagonal(
    lower: &[f64],
    diag: &[f64],
    upper: &[f64],
    rhs: &[C64],
) -> Vec<C64> 

The solution of the tridiagonal system lower[i] u[i-1] + diag[i] u[i] + upper[i] u[i+1] = rhs[i] (lower[0] and the last upper unused), by the Thomas algorithm without pivoting: for diagonally dominant systems, as those of cubic splines.

Examples

use comnum_math::{C64, linalg::solve_tridiagonal};
// [2 1 0; 1 2 1; 0 1 2] u = [3, 4, 3] has u = [1, 1, 1]
let r = [3.0, 4.0, 3.0].map(|v| C64::new(v, 0.0));
let u = solve_tridiagonal(&[0.0, 1.0, 1.0], &[2.0; 3], &[1.0, 1.0, 0.0], &r);
assert!(u.iter().all(|v| (v - C64::new(1.0, 0.0)).norm() < 1e-15));