Function thomas 

pub fn thomas(
    n: usize,
    lower: &[f64],
    diag: &[f64],
    upper: &[f64],
    rhs: &[f64],
    x: &mut [f64],
)

Thomas algorithm (tridiagonal matrix algorithm) for solving a 1-D tridiagonal system.

Solves: lower[i]*x[i-1] + diag[i]*x[i] + upper[i]*x[i+1] = rhs[i]

Boundary conventions:

• lower[0] is unused (no element below the first row).
• upper[n-1] is unused (no element above the last row).

Arguments

• n - System size.
• lower - Sub-diagonal coefficients (length n; lower[0] ignored).
• diag - Main diagonal coefficients (length n).
• upper - Super-diagonal coefficients (length n; upper[n-1] ignored).
• rhs - Right-hand side vector (length n).
• x - Output solution vector (length n).

Panics

Panics in debug mode if n == 0.

Examples

use solver::linalg::adi::thomas;

// Solve the system:
//  [ 2 -1  0 ] [x0]   [1]
//  [-1  2 -1 ] [x1] = [0]
//  [ 0 -1  2 ] [x2]   [1]
// Solution: x = [1, 1, 1]
let lo  = [0.0, -1.0, -1.0];
let di  = [2.0,  2.0,  2.0];
let up  = [-1.0, -1.0, 0.0];
let rhs = [1.0,  0.0,  1.0];
let mut x = [0.0f64; 3];
thomas(3, &lo, &di, &up, &rhs, &mut x);
assert!((x[0] - 1.0).abs() < 1e-12);
assert!((x[1] - 1.0).abs() < 1e-12);
assert!((x[2] - 1.0).abs() < 1e-12);