326 lines
10 KiB
Rust
326 lines
10 KiB
Rust
//! Shared test helpers for the ruvector-solver integration test suite.
|
|
//!
|
|
//! Provides deterministic random matrix generators, dense reference solvers,
|
|
//! and floating-point comparison utilities used across all test modules.
|
|
|
|
use ruvector_solver::types::CsrMatrix;
|
|
|
|
// ---------------------------------------------------------------------------
|
|
// Random number generator (simple LCG for deterministic reproducibility)
|
|
// ---------------------------------------------------------------------------
|
|
|
|
/// A minimal linear congruential generator for deterministic test data.
|
|
///
|
|
/// Uses the Numerical Recipes LCG parameters. Not cryptographically secure,
|
|
/// but perfectly adequate for generating reproducible test matrices.
|
|
pub struct Lcg {
|
|
state: u64,
|
|
}
|
|
|
|
impl Lcg {
|
|
/// Create a new LCG with the given seed.
|
|
pub fn new(seed: u64) -> Self {
|
|
Self { state: seed }
|
|
}
|
|
|
|
/// Generate the next u64 value.
|
|
pub fn next_u64(&mut self) -> u64 {
|
|
self.state = self
|
|
.state
|
|
.wrapping_mul(6364136223846793005)
|
|
.wrapping_add(1442695040888963407);
|
|
self.state
|
|
}
|
|
|
|
/// Generate a uniform f64 in [0, 1).
|
|
pub fn next_f64(&mut self) -> f64 {
|
|
(self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
|
|
}
|
|
|
|
/// Generate a uniform f64 in [lo, hi).
|
|
pub fn next_f64_range(&mut self, lo: f64, hi: f64) -> f64 {
|
|
lo + (hi - lo) * self.next_f64()
|
|
}
|
|
}
|
|
|
|
// ---------------------------------------------------------------------------
|
|
// Matrix generators
|
|
// ---------------------------------------------------------------------------
|
|
|
|
/// Generate a random diagonally dominant CSR matrix of dimension `n`.
|
|
///
|
|
/// Each row has approximately `density * n` non-zero off-diagonal entries
|
|
/// (at least 1). The diagonal entry is set to `1 + sum_of_abs_off_diag`
|
|
/// to guarantee strict diagonal dominance.
|
|
///
|
|
/// The resulting matrix is suitable for Neumann and Jacobi solvers.
|
|
pub fn random_diag_dominant_csr(n: usize, density: f64, seed: u64) -> CsrMatrix<f64> {
|
|
let mut rng = Lcg::new(seed);
|
|
let mut entries: Vec<(usize, usize, f64)> = Vec::new();
|
|
|
|
for i in 0..n {
|
|
let mut off_diag_sum = 0.0f64;
|
|
|
|
for j in 0..n {
|
|
if i == j {
|
|
continue;
|
|
}
|
|
if rng.next_f64() < density {
|
|
let val = rng.next_f64_range(-1.0, 1.0);
|
|
entries.push((i, j, val));
|
|
off_diag_sum += val.abs();
|
|
}
|
|
}
|
|
|
|
// Ensure at least one off-diagonal entry per row for non-trivial testing.
|
|
if off_diag_sum == 0.0 && n > 1 {
|
|
let j = (i + 1) % n;
|
|
let val = rng.next_f64_range(0.1, 0.5);
|
|
entries.push((i, j, val));
|
|
off_diag_sum = val;
|
|
}
|
|
|
|
// Diagonal: strictly dominant.
|
|
let diag_val = off_diag_sum + 1.0 + rng.next_f64();
|
|
entries.push((i, i, diag_val));
|
|
}
|
|
|
|
CsrMatrix::<f64>::from_coo(n, n, entries)
|
|
}
|
|
|
|
/// Generate a random graph Laplacian CSR matrix of dimension `n`.
|
|
///
|
|
/// A graph Laplacian `L = D - A` where:
|
|
/// - `A` is the adjacency matrix of a random undirected graph.
|
|
/// - `D` is the degree matrix.
|
|
/// - Each row sums to zero (L * ones = 0).
|
|
///
|
|
/// The resulting matrix is symmetric positive semi-definite.
|
|
pub fn random_laplacian_csr(n: usize, density: f64, seed: u64) -> CsrMatrix<f64> {
|
|
let mut rng = Lcg::new(seed);
|
|
|
|
// Build symmetric adjacency: for i < j, randomly include edge (i,j).
|
|
let mut adj = vec![vec![0.0f64; n]; n];
|
|
|
|
for i in 0..n {
|
|
for j in (i + 1)..n {
|
|
if rng.next_f64() < density {
|
|
let weight = rng.next_f64_range(0.1, 2.0);
|
|
adj[i][j] = weight;
|
|
adj[j][i] = weight;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Ensure the graph is connected: add a path 0-1-2-...-n-1.
|
|
for i in 0..n.saturating_sub(1) {
|
|
if adj[i][i + 1] == 0.0 {
|
|
let weight = rng.next_f64_range(0.1, 1.0);
|
|
adj[i][i + 1] = weight;
|
|
adj[i + 1][i] = weight;
|
|
}
|
|
}
|
|
|
|
// Build Laplacian: L[i][j] = -A[i][j] for i != j, L[i][i] = sum_j A[i][j].
|
|
let mut entries: Vec<(usize, usize, f64)> = Vec::new();
|
|
|
|
for i in 0..n {
|
|
let mut degree = 0.0f64;
|
|
for j in 0..n {
|
|
if i != j && adj[i][j] != 0.0 {
|
|
entries.push((i, j, -adj[i][j]));
|
|
degree += adj[i][j];
|
|
}
|
|
}
|
|
entries.push((i, i, degree));
|
|
}
|
|
|
|
CsrMatrix::<f64>::from_coo(n, n, entries)
|
|
}
|
|
|
|
/// Generate a random SPD (symmetric positive definite) matrix.
|
|
///
|
|
/// Constructs `A = B^T B + epsilon * I` where `B` has random entries,
|
|
/// guaranteeing positive definiteness.
|
|
pub fn random_spd_csr(n: usize, density: f64, seed: u64) -> CsrMatrix<f64> {
|
|
let mut rng = Lcg::new(seed);
|
|
|
|
// Build a random dense matrix B, then compute A = B^T B + eps * I.
|
|
// For efficiency with CSR, we do this differently: build a sparse
|
|
// symmetric matrix and add a diagonal shift.
|
|
let mut dense = vec![vec![0.0f64; n]; n];
|
|
|
|
for i in 0..n {
|
|
for j in i..n {
|
|
if i == j || rng.next_f64() < density {
|
|
let val = rng.next_f64_range(-1.0, 1.0);
|
|
dense[i][j] += val;
|
|
if i != j {
|
|
dense[j][i] += val;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Compute A = M^T M where M = dense (makes it PSD).
|
|
let mut a = vec![vec![0.0f64; n]; n];
|
|
for i in 0..n {
|
|
for j in 0..n {
|
|
let mut sum = 0.0;
|
|
for k in 0..n {
|
|
sum += dense[k][i] * dense[k][j];
|
|
}
|
|
a[i][j] = sum;
|
|
}
|
|
}
|
|
|
|
// Add diagonal shift to ensure positive definiteness.
|
|
for i in 0..n {
|
|
a[i][i] += 1.0;
|
|
}
|
|
|
|
// Convert to COO.
|
|
let mut entries: Vec<(usize, usize, f64)> = Vec::new();
|
|
for i in 0..n {
|
|
for j in 0..n {
|
|
if a[i][j].abs() > 1e-15 {
|
|
entries.push((i, j, a[i][j]));
|
|
}
|
|
}
|
|
}
|
|
|
|
CsrMatrix::<f64>::from_coo(n, n, entries)
|
|
}
|
|
|
|
/// Generate a deterministic random vector of length `n`.
|
|
pub fn random_vector(n: usize, seed: u64) -> Vec<f64> {
|
|
let mut rng = Lcg::new(seed);
|
|
(0..n).map(|_| rng.next_f64_range(-1.0, 1.0)).collect()
|
|
}
|
|
|
|
/// Build a simple undirected graph as a CSR adjacency matrix.
|
|
///
|
|
/// Each entry `(u, v)` in `edges` creates entries `A[u][v] = 1` and
|
|
/// `A[v][u] = 1`.
|
|
pub fn adjacency_from_edges(n: usize, edges: &[(usize, usize)]) -> CsrMatrix<f64> {
|
|
let mut entries: Vec<(usize, usize, f64)> = Vec::new();
|
|
for &(u, v) in edges {
|
|
entries.push((u, v, 1.0));
|
|
if u != v {
|
|
entries.push((v, u, 1.0));
|
|
}
|
|
}
|
|
CsrMatrix::<f64>::from_coo(n, n, entries)
|
|
}
|
|
|
|
// ---------------------------------------------------------------------------
|
|
// Dense reference solver
|
|
// ---------------------------------------------------------------------------
|
|
|
|
/// Solve `Ax = b` using dense Gaussian elimination with partial pivoting.
|
|
///
|
|
/// This is an O(n^3) reference solver used only for small test problems
|
|
/// to verify iterative solver accuracy.
|
|
///
|
|
/// # Panics
|
|
///
|
|
/// Panics if the matrix is singular or dimensions are inconsistent.
|
|
pub fn dense_solve(matrix: &CsrMatrix<f64>, rhs: &[f64]) -> Vec<f64> {
|
|
let n = matrix.rows;
|
|
assert_eq!(n, matrix.cols, "dense_solve requires a square matrix");
|
|
assert_eq!(rhs.len(), n, "rhs length must match matrix dimension");
|
|
|
|
// Convert CSR to dense augmented matrix [A | b].
|
|
let mut aug = vec![vec![0.0f64; n + 1]; n];
|
|
for i in 0..n {
|
|
aug[i][n] = rhs[i];
|
|
let start = matrix.row_ptr[i];
|
|
let end = matrix.row_ptr[i + 1];
|
|
for idx in start..end {
|
|
let j = matrix.col_indices[idx];
|
|
aug[i][j] = matrix.values[idx];
|
|
}
|
|
}
|
|
|
|
// Forward elimination with partial pivoting.
|
|
for col in 0..n {
|
|
// Find pivot.
|
|
let mut max_row = col;
|
|
let mut max_val = aug[col][col].abs();
|
|
for row in (col + 1)..n {
|
|
if aug[row][col].abs() > max_val {
|
|
max_val = aug[row][col].abs();
|
|
max_row = row;
|
|
}
|
|
}
|
|
assert!(max_val > 1e-15, "matrix is singular or near-singular");
|
|
aug.swap(col, max_row);
|
|
|
|
// Eliminate.
|
|
let pivot = aug[col][col];
|
|
for row in (col + 1)..n {
|
|
let factor = aug[row][col] / pivot;
|
|
for j in col..=n {
|
|
aug[row][j] -= factor * aug[col][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Back substitution.
|
|
let mut x = vec![0.0f64; n];
|
|
for i in (0..n).rev() {
|
|
let mut sum = aug[i][n];
|
|
for j in (i + 1)..n {
|
|
sum -= aug[i][j] * x[j];
|
|
}
|
|
x[i] = sum / aug[i][i];
|
|
}
|
|
|
|
x
|
|
}
|
|
|
|
// ---------------------------------------------------------------------------
|
|
// Floating-point comparison utilities
|
|
// ---------------------------------------------------------------------------
|
|
|
|
/// Compute the L2 norm of a vector.
|
|
pub fn l2_norm(v: &[f64]) -> f64 {
|
|
v.iter().map(|&x| x * x).sum::<f64>().sqrt()
|
|
}
|
|
|
|
/// Compute the L2 distance between two vectors.
|
|
pub fn l2_distance(a: &[f64], b: &[f64]) -> f64 {
|
|
assert_eq!(a.len(), b.len(), "vectors must have same length");
|
|
a.iter()
|
|
.zip(b.iter())
|
|
.map(|(&ai, &bi)| (ai - bi) * (ai - bi))
|
|
.sum::<f64>()
|
|
.sqrt()
|
|
}
|
|
|
|
/// Compute the relative error ||approx - exact|| / ||exact||.
|
|
///
|
|
/// Returns absolute error if the exact solution has zero norm.
|
|
pub fn relative_error(approx: &[f64], exact: &[f64]) -> f64 {
|
|
let exact_norm = l2_norm(exact);
|
|
let error = l2_distance(approx, exact);
|
|
if exact_norm > 1e-15 {
|
|
error / exact_norm
|
|
} else {
|
|
error
|
|
}
|
|
}
|
|
|
|
/// Compute the residual `b - A*x` for a sparse system.
|
|
pub fn compute_residual(matrix: &CsrMatrix<f64>, x: &[f64], rhs: &[f64]) -> Vec<f64> {
|
|
let n = matrix.rows;
|
|
let mut ax = vec![0.0f64; n];
|
|
matrix.spmv(x, &mut ax);
|
|
(0..n).map(|i| rhs[i] - ax[i]).collect()
|
|
}
|
|
|
|
/// Convert an f32 solution vector to f64 for comparison.
|
|
pub fn f32_to_f64(v: &[f32]) -> Vec<f64> {
|
|
v.iter().map(|&x| x as f64).collect()
|
|
}
|