Struct CsrMatrix 

pub struct CsrMatrix {
    pub values: Vec<f64>,
    pub col_indices: Vec<usize>,
    pub row_ptr: Vec<usize>,
    pub size: usize,
    pub diag_indices: Vec<usize>,
}

Compressed Sparse Row Matrix

Examples

use solver::linalg::csr::CsrMatrix;

// Create a 3x3 identity matrix
let mut mat = CsrMatrix::new(3, 3);

// Row 0
mat.add_entry(0, 1.0);
mat.finish_row();

// Row 1
mat.add_entry(1, 1.0);
mat.finish_row();

// Row 2
mat.add_entry(2, 1.0);
mat.finish_row();

assert_eq!(mat.values, vec![1.0, 1.0, 1.0]);
assert_eq!(mat.col_indices, vec![0, 1, 2]);
assert_eq!(mat.row_ptr, vec![0, 1, 2, 3]);

Fields

values: Vec<f64>

col_indices: Vec<usize>

row_ptr: Vec<usize>

size: usize

diag_indices: Vec<usize>

Implementations

impl CsrMatrix

pub fn new(size: usize, estimated_entries: usize) -> Self

pub fn add_entry(&mut self, col: usize, val: f64)

pub fn finish_row(&mut self)

pub fn solve_tridiagonal_system(&self, b: &[f64]) -> Vec<f64>

Solves the linear system Ax = b using LAPACK’s DGTSV (General Tridiagonal Solve). Assumes the matrix is tridiagonal.

Examples

use solver::linalg::csr::CsrMatrix;

// Solve A x = b where A is:
// [ 2 -1  0 ]
// [-1  2 -1 ]
// [ 0 -1  2 ]
// and b = [1, 0, 1]
// Solution should be x = [1, 1, 1]

let mut mat = CsrMatrix::new(3, 9);

// Row 0
mat.add_entry(0, 2.0);
mat.add_entry(1, -1.0);
mat.finish_row();

// Row 1
mat.add_entry(0, -1.0);
mat.add_entry(1, 2.0);
mat.add_entry(2, -1.0);
mat.finish_row();

// Row 2
mat.add_entry(1, -1.0);
mat.add_entry(2, 2.0);
mat.finish_row();

let b = vec![1.0, 0.0, 1.0];
let x = mat.solve_tridiagonal_system(&b);

assert!((x[0] - 1.0).abs() < 1e-6);
assert!((x[1] - 1.0).abs() < 1e-6);
assert!((x[2] - 1.0).abs() < 1e-6);

impl CsrMatrix

pub fn solve_sor( &self, b: &[f64], x: &mut [f64], tol: f64, max_iter: usize, omega: f64, ) -> usize

Successive Over-Relaxation (SOR)

The Successive Over-Relaxation (SOR) method is an iterative technique for solving a linear system $Ax = b$. It is a variant of the Gauss-Seidel method that uses a relaxation factor $\omega$ to accelerate convergence.

Mathematical Formulation

For a square matrix $A$ with elements $a_{ij}$, the solution $x$ is updated in place, element by element. The update rule for the $i$-th component at iteration $k+1$ is:

$$ x_i^{(k+1)} = (1-\omega),x_i^{(k)} + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j < i} a_{ij} x_j^{(k+1)} - \sum_{j > i} a_{ij} x_j^{(k)} \right) $$

Here $a_{ii} \neq 0$ is required; the implementation skips rows with a zero diagonal entry. Because the method overwrites $x_i$ immediately, the most recently updated values are always used on the right-hand side.

Defining the Gauss-Seidel iterate

$$ x_i^{(GS)} = \frac{1}{a_{ii}} \left(b_i - \sum_{j \neq i} a_{ij} x_j \right), $$

where $x_j$ takes the most recent available value, the SOR update is the weighted combination

$$ x_i^{(k+1)} = (1-\omega),x_i^{(k)} + \omega,x_i^{(GS)}. $$

Convergence

SOR converges for an arbitrary initial guess when $A$ is symmetric positive definite and $0 < \omega < 2$. The iteration matrix is $\omega$-dependent; no value outside this range is guaranteed to converge.

• $\omega = 1$: reduces to Gauss-Seidel.
• $0 < \omega < 1$: under-relaxation; used to stabilize systems that oscillate or diverge under Gauss-Seidel.
• $1 < \omega < 2$: over-relaxation; used to accelerate convergence for systems that already converge.
• $\omega \ge 2$: outside the guaranteed convergence interval for symmetric positive definite systems.

The HJB discretizations in this project typically yield symmetric M-matrices, which are symmetric positive definite; when the system is merely diagonally dominant, convergence for $0 < \omega < 2$ is not guaranteed and the value should be verified numerically. For these systems $\omega \in (1, 2)$ is the usual choice. The optimal $\omega$ depends on the spectral radius of the Gauss-Seidel iteration matrix and is not derived here.

Validation targets are the linear systems in solver/tests/infra/linalg.rs (test_sor_identity, test_sor_diagonal_system, test_sor_tridiagonal_system, and test_sor_3x3_nontrivial), each checked against a known exact solution.

Examples

use solver::linalg::csr::CsrMatrix;

// Solve A x = b where A is:
// [ 4 -1 ]
// [-1  4 ]
// and b = [3, 3]
// Solution should be x = [1, 1]

let mut mat = CsrMatrix::new(2, 4);

// Row 0
mat.add_entry(0, 4.0);
mat.add_entry(1, -1.0);
mat.finish_row();

// Row 1
mat.add_entry(0, -1.0);
mat.add_entry(1, 4.0);
mat.finish_row();

let b = vec![3.0, 3.0];
let mut x = vec![0.0, 0.0]; // Initial guess

let iters = mat.solve_sor(&b, &mut x, 1e-6, 100, 1.0);

assert!((x[0] - 1.0).abs() < 1e-5);
assert!((x[1] - 1.0).abs() < 1e-5);