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 of equations $Ax = b$. It is a variant of the Gauss-Seidel method that uses a relaxation factor $\omega$ to speed up convergence.

Mathematical Formulation

For a square matrix $A$ with elements $a_{ij}$, the solution $x$ is updated element-wise. 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) $$

This can be rewritten as a weighted average of the previous iterate and the Gauss-Seidel iterate:

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

where $x_i^{(GS)}$ is the value computed by the Gauss-Seidel method:

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

(using the most recent values for $x_j$).

Convergence Properties

The convergence of the SOR method depends on the choice of the relaxation factor $\omega$:

• $\omega = 1$: The method reduces to the Gauss-Seidel method.
• $0 < \omega < 1$: This is under-relaxation. It can help converge systems that would otherwise diverge or oscillate.
• $1 < \omega < 2$: This is over-relaxation. It is typically used to accelerate convergence for systems that are already convergent (like those arising from elliptic PDEs).
• $\omega \ge 2$: The method generally fails to converge.

For the HJB equation discretization (which often results in an M-matrix or diagonally dominant matrix), a value of $\omega \in (1, 2)$ is usually optimal.

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);