Numerics extras

Scientific work defines VECMAT_USE_F64 before the include. -DVECMAT_USE_F64=ON does that for this build's tests and benchmarks; it does not drop name32.

GEMM

BLAS-style dense multiply C = alpha * op(A) * op(B) + beta * C. Row-major and column-major layouts.

Packed 8×8 ukernel, TLS pack workspace MC=NC=KC=128, persistent thread pool. Cap threads with vm_gemm_set_threads(n) or VECMAT_GEMM_THREADS. fp16 / bf16 are not in this release. This is not a BLAS.

Contract: vm_gemm matches vm_gemm_ref for trans / no-trans, row / col, odd sizes, beta ≠ 0, and batch shared-B, on a size band of about 32–512.

Column-major C = A * B for square panels:

vm_gemm(C, n, A, n, B, n, n, n, n,
        1.0f, 0.0f, false, false, VM_LAYOUT_COL_MAJOR);

Dense linear algebra

Heap vm_mat (M×N, column-major). LU with partial pivoting, Householder QR, thin one-sided Jacobi SVD, Cholesky. Scale-aware cutoffs tol = n * VECMAT_EPSILON * max|A|. Rank-deficient work returns false (an ok flag), not a silent NaN. Integer matNi_inverse is truncated, not modular inverse — the declaration says so.

vm_mat A = vm_mat_alloc(3, 3);
vm_mat_set(&A, 0, 0, 2.0);
vm_mat_set(&A, 1, 1, 2.0);
vm_mat_set(&A, 2, 2, 2.0);
vm_mat inv;
bool ok = vm_mat_inverse(&inv, &A);
vm_mat_free(&A);
vm_mat_free(&inv);

3×3 factorizations

mat3_inverse_ok returns false on a singular matrix (scale-aware cutoff). mat3_sym_eigen is the analytical symmetric 3×3 eigensolve, eigenvalues sorted ascending, eigenvectors on request. mat3_polar is the polar decomposition F = R U (Newton on C = Fᵀ F). A reflection is a failure.

Dense spectral primitives

include/vecmat/spectral.h:

  • vm_matvec, vm_symv — dense matvec (symmetric reads the lower triangle)
  • vm_syev — Jacobi symmetric eigensolve, eigenvalues descending, A overwritten
  • vm_laplacian — unnormalized graph Laplacian L = D − A
  • vm_syev_topk — thin wrapper for the leading k eigenpairs

Sparse systems

  • vm_spmat — square CSR from triplets
  • vm_spmv, vm_cg (SPD), vm_bicgstab (nonsymmetric)
  • Left preconditioners: Jacobi, SSOR (ω = 1), IC(0) (falls back to Jacobi on pivot breakdown)
  • vm_ksp_info reports iters, rel_res, ok
  • Relative residual is ||r|| / max(||b||, ε)

Time integration and grid primitives

Semi-implicit Euler on one particle. x, v, and a are length n.

vm_float_t x[3] = {0.0f, 1.0f, 0.0f};
vm_float_t v[3] = {1.0f, 0.0f, 0.0f};
const vm_float_t a[3] = {0.0f, -9.81f, 0.0f};
vm_euler_semi(x, v, a, 3, 1.0f / 60.0f);

Computer graphics (core)

Clip-space helpers

Build projection and view matrices for different graphics APIs and depth conventions.

Perspective projections — camera frustum matrices (radians; _deg if FOV is in degrees). The unsuffixed mat4_perspective / mat4_perspective_fov / mat4_perspective_infinite helpers also take radians. Use mat4_perspective_deg (and friends) for degrees:

Orthographic projections — parallel projection matrices from frustum bounds:

Look-at view matrices — world-to-view transforms from eye, target, and up:

Look-from-direction view matrices — the same basis as look-at, but the camera aims along a direction (FPS / fly camera, no target point):

Infinite / reverse-Z projections — infinite far plane, optionally with reversed depth (near → 1, infinity → 0 on ZO):

Viewport, world ↔ window — NDC to a pixel box and back. Geometric vec3_project (onto a direction) is unchanged:

Affine inverse and normal matrices — skip the 4×4 adjugate when the transform is [A t; 0 1]:

Clip conventions — handedness + depth range selectors used by the *_clip helpers:

A Vulkan-style camera is a projection times a view. Angles are radians; VM_DEG converts a degree literal.

const vector3 eye = vec3(0.0f, 1.6f, 4.0f);
const vector3 at  = vec3(0.0f, 1.0f, 0.0f);
const vector3 up  = vec3(0.0f, 1.0f, 0.0f);
const matrix4 view = mat4_look_at_rh(eye, at, up);
const matrix4 proj = mat4_perspective_rh_zo(VM_DEG(60.0), 16.0f / 9.0f, 0.1f, 200.0f);
const matrix4 view_proj = mat4_mul(proj, view);

Rotation helpers

Build 4×4 rotation matrices from axis angles in radians (mat4_rotation / mat4_rotation_x / mat4_rotation_y / mat4_rotation_z). Use *_deg or VM_DEG(...) when the angle is in degrees: