src/features/sparse.c file

Classes

struct vm_trip
Sparse matrix triplet (row, column, value).
struct vm_ic0
Incomplete Cholesky IC(0) factor for preconditioning.
struct vm_prec
Left preconditioner for iterative solvers (CG / BiCGSTAB).

Functions

void vm_spmat_init(vm_spmat* A)
Initialize a sparse matrix to empty.
void vm_spmat_free(vm_spmat* A)
Free sparse matrix storage and reset it.
static auto vm_trip_cmp(const void* pa, const void* pb) -> int
Compare triplets by (row, col) for qsort.
static auto vm_vec_dot(const vm_float_t* a, const vm_float_t* b, const int n) -> vm_float_t
Dot product of two vectors.
static auto vm_vec_norm2(const vm_float_t* a, const int n) -> vm_float_t
Euclidean (L2) norm of a vector.
static void vm_vec_axpy(vm_float_t* y, const vm_float_t a, const vm_float_t* x, const int n)
y := y + a * x (AXPY).
static void vm_vec_copy(vm_float_t* dst, const vm_float_t* src, const int n)
Copy a vector.
static void vm_vec_zero(vm_float_t* y, const int n)
Set a vector to zero.
auto vm_spmat_from_triplets(vm_spmat* A, const int n, const int nnz, const int* row, const int* col, const vm_float_t* val) -> bool
Build a square CSR matrix from unsorted (row, col, val) triplets.
void vm_spmv(vm_float_t* y, const vm_spmat* A, const vm_float_t* x)
Sparse matrix–vector product y = A x.
auto vm_spmat_diag(const vm_spmat* A, vm_float_t* d) -> bool
Extract the main diagonal of A into d.
static void vm_jacobi_apply(vm_float_t* z, const vm_float_t* d, const vm_float_t* r, const int n)
Apply Jacobi preconditioner: z_i = r_i / d_i.
static void vm_ssor_apply(vm_float_t* z, const vm_spmat* A, const vm_float_t* d, const vm_float_t* r, vm_float_t* tmp)
Apply symmetric Gauss–Seidel (SSOR, ω = 1): (D+L) D⁻¹ (D+U) z = r.
static void vm_ic0_free(vm_ic0* ic)
Free IC(0) factor data.
static auto vm_ic0_factor(vm_ic0* ic, const vm_spmat* A) -> bool
Build incomplete Cholesky IC(0) factor of A (lower triangle).
static void vm_ic0_apply(vm_float_t* z, const vm_ic0* ic, const vm_float_t* r, vm_float_t* tmp)
Apply IC(0): solve L Lᵀ z = r.
static void vm_prec_free(vm_prec* p)
Free preconditioner data.
static auto vm_prec_setup(vm_prec* p, const vm_spmat* A, const vm_ksp_prec_t kind) -> bool
Set up a left preconditioner for A.
static void vm_prec_apply(vm_float_t* z, const vm_prec* p, const vm_float_t* r)
Apply preconditioner: z ≈ M⁻¹ r.
static void vm_ksp_set(vm_ksp_info* info, const int iters, const vm_float_t relres, const bool ok)
Fill KSP result info.
auto vm_cg(const vm_spmat* A, const vm_float_t* b, vm_float_t* x, vm_float_t tol, int max_iter, const vm_ksp_prec_t pre_cond, vm_ksp_info* info) -> bool
Conjugate gradient for SPD A x = b.
auto vm_bicgstab(const vm_spmat* A, const vm_float_t* b, vm_float_t* x, vm_float_t tol, int max_iter, const vm_ksp_prec_t pre_cond, vm_ksp_info* info) -> bool
BiCGSTAB for general (possibly nonsymmetric) A x = b.
static auto vm_mat_at(const vm_mat* m, const int r, const int c) -> vm_float_t*
Pointer to entry A(r, c) in column-major dense storage.
auto vm_chol_factor(vm_mat* A) -> bool
In-place dense Cholesky A = L Lᵀ (lower triangle overwritten).
auto vm_chol_solve(const vm_mat* L, const vm_float_t* b, vm_float_t* x) -> bool
Solve L Lᵀ x = b after vm_chol_factor.

Function documentation

void vm_spmat_init(vm_spmat* A)

Initialize a sparse matrix to empty.

Parameters
A Sparse matrix (may be NULL).

void vm_spmat_free(vm_spmat* A)

Free sparse matrix storage and reset it.

Parameters
A Sparse matrix (may be NULL).

static int vm_trip_cmp(const void* pa, const void* pb)

Compare triplets by (row, col) for qsort.

Parameters
pa First triplet.
pb Second triplet.
Returns -1, 0, or 1.

static vm_float_t vm_vec_dot(const vm_float_t* a, const vm_float_t* b, const int n)

Dot product of two vectors.

Parameters
a First vector.
b Second vector.
n Length.
Returns a·b.

static vm_float_t vm_vec_norm2(const vm_float_t* a, const int n)

Euclidean (L2) norm of a vector.

Parameters
a Vector.
n Length.
Returns ||a||₂.

static void vm_vec_axpy(vm_float_t* y, const vm_float_t a, const vm_float_t* x, const int n)

y := y + a * x (AXPY).

Parameters
y Input/output vector.
a Scalar.
x Input vector.
n Length.

static void vm_vec_copy(vm_float_t* dst, const vm_float_t* src, const int n)

Copy a vector.

Parameters
dst Destination.
src Source.
n Length.

static void vm_vec_zero(vm_float_t* y, const int n)

Set a vector to zero.

Parameters
y Vector.
n Length.

bool vm_spmat_from_triplets(vm_spmat* A, const int n, const int nnz, const int* row, const int* col, const vm_float_t* val)

Build a square CSR matrix from unsorted (row, col, val) triplets.

Parameters
A Output sparse matrix.
n Matrix size (n×n).
nnz Number of input triplets.
row Row indices (length nnz).
col Column indices (length nnz).
val Values (length nnz).
Returns true on success, false on error.

Duplicate (i, j) entries are summed. Out-of-range indices are skipped. On success, existing storage in A is freed and replaced.

void vm_spmv(vm_float_t* y, const vm_spmat* A, const vm_float_t* x)

Sparse matrix–vector product y = A x.

Parameters
y Output vector (length A->n).
A CSR matrix.
x Input vector (length A->n).

bool vm_spmat_diag(const vm_spmat* A, vm_float_t* d)

Extract the main diagonal of A into d.

Parameters
A CSR matrix.
d Output diagonal (length A->n).
Returns true on success, false on error.

Missing diagonal entries are set to 0.

static void vm_jacobi_apply(vm_float_t* z, const vm_float_t* d, const vm_float_t* r, const int n)

Apply Jacobi preconditioner: z_i = r_i / d_i.

Parameters
z Output vector.
d Diagonal.
r Residual.
n Length.

static void vm_ssor_apply(vm_float_t* z, const vm_spmat* A, const vm_float_t* d, const vm_float_t* r, vm_float_t* tmp)

Apply symmetric Gauss–Seidel (SSOR, ω = 1): (D+L) D⁻¹ (D+U) z = r.

Parameters
z Output vector.
A CSR matrix.
d Diagonal of A.
r Residual.
tmp Workspace (length n).

static void vm_ic0_free(vm_ic0* ic)

Free IC(0) factor data.

Parameters
ic IC(0) handle (may be NULL).

static bool vm_ic0_factor(vm_ic0* ic, const vm_spmat* A)

Build incomplete Cholesky IC(0) factor of A (lower triangle).

Parameters
ic Output IC(0) handle.
A SPD-like CSR matrix with a full diagonal.
Returns true on success, false on breakdown or error.

static void vm_ic0_apply(vm_float_t* z, const vm_ic0* ic, const vm_float_t* r, vm_float_t* tmp)

Apply IC(0): solve L Lᵀ z = r.

Parameters
z Output vector.
ic IC(0) factor.
r Residual.
tmp Workspace (length n).

static void vm_prec_free(vm_prec* p)

Free preconditioner data.

Parameters
p Preconditioner (may be NULL).

static bool vm_prec_setup(vm_prec* p, const vm_spmat* A, const vm_ksp_prec_t kind)

Set up a left preconditioner for A.

Parameters
p Output preconditioner.
A CSR matrix.
kind NONE, JACOBI, SSOR, or IC0.
Returns true on success, false on allocation/setup failure.

IC(0) breakdown falls back to Jacobi.

static void vm_prec_apply(vm_float_t* z, const vm_prec* p, const vm_float_t* r)

Apply preconditioner: z ≈ M⁻¹ r.

Parameters
z Output vector.
p Preconditioner.
r Residual.

static void vm_ksp_set(vm_ksp_info* info, const int iters, const vm_float_t relres, const bool ok)

Fill KSP result info.

Parameters
info Output info (may be NULL).
iters Iteration count.
relres Relative residual.
ok Convergence flag.

bool vm_cg(const vm_spmat* A, const vm_float_t* b, vm_float_t* x, vm_float_t tol, int max_iter, const vm_ksp_prec_t pre_cond, vm_ksp_info* info)

Conjugate gradient for SPD A x = b.

Parameters
A SPD CSR matrix.
b Right-hand side (length A->n).
x Initial guess / solution (length A->n).
tol Relative residual tolerance (≤0 picks a default).
max_iter Max iterations (≤0 defaults to n).
pre_cond NONE, JACOBI, SSOR, or IC0.
info Optional solver stats (may be NULL).
Returns true if converged, false otherwise.

x is the initial guess and the solution. tol is relative residual ||r|| / max(||b||, ε).

bool vm_bicgstab(const vm_spmat* A, const vm_float_t* b, vm_float_t* x, vm_float_t tol, int max_iter, const vm_ksp_prec_t pre_cond, vm_ksp_info* info)

BiCGSTAB for general (possibly nonsymmetric) A x = b.

Parameters
A CSR matrix.
b Right-hand side (length A->n).
x Initial guess / solution (length A->n).
tol Relative residual tolerance (≤0 picks a default).
max_iter Max iterations (≤0 defaults to 2n).
pre_cond NONE, JACOBI, SSOR, or IC0.
info Optional solver stats (may be NULL).
Returns true if converged, false otherwise.

Same calling convention as vm_cg. Jacobi / SSOR / IC0 are left preconditioners; IC0 still expects an SPD-like diagonal.

static vm_float_t* vm_mat_at(const vm_mat* m, const int r, const int c)

Pointer to entry A(r, c) in column-major dense storage.

Parameters
m Dense matrix.
r Row index.
c Column index.
Returns Pointer to the entry.

bool vm_chol_factor(vm_mat* A)

In-place dense Cholesky A = L Lᵀ (lower triangle overwritten).

Parameters
A Square dense matrix (destroyed / factored in place).
Returns true on success, false if not SPD or invalid input.

A must be square SPD. The strict upper triangle is left untouched.

bool vm_chol_solve(const vm_mat* L, const vm_float_t* b, vm_float_t* x)

Solve L Lᵀ x = b after vm_chol_factor.

Parameters
L Factored lower triangle from vm_chol_factor.
b Right-hand side (length n).
x Solution (length n).
Returns true on success, false on error.