outram_foam_basic_lib/krylov/
vecops.rs

1// SPDX-License-Identifier: GPL-3.0-only
2// Copyright (C) 2026 OUTRAM PARK contributors
3//
4// This file is part of OUTRAM PARK.
5//
6// OUTRAM PARK is free software: you can redistribute it and/or modify it
7// under the terms of the GNU General Public License as published by the
8// Free Software Foundation, either version 3 of the License, or (at your
9// option) any later version.
10//
11// OUTRAM PARK is distributed in the hope that it will be useful, but
12// WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14// General Public License for more details.
15//
16// You should have received a copy of the GNU General Public License along
17// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.
18
19//! Dense BLAS-1 vector primitives for the Krylov solvers.
20//!
21//! Pure-Rust replacements for the handful of level-1 BLAS operations the
22//! iterative solvers need. All operands are dimensionless `&[f64]` slices whose
23//! length equals the number of unknowns (mesh cells); no `uom` typing is applied
24//! here because a Krylov subspace mixes residuals, search directions, and
25//! solution increments that share no single physical dimension.
26//!
27//! Every function is `O(n)` in the slice length and allocation-free (results are
28//! either scalars or written in place), so they are safe to call inside the
29//! innermost solver loops.
30
31/// Euclidean inner product `Σ_i a_i · b_i` (dimensionless).
32///
33/// Both slices must have the same length; a mismatch panics via
34/// `debug_assert`. Valid for any finite inputs; returns `0.0` for empty slices.
35pub fn dot(a: &[f64], b: &[f64]) -> f64 {
36    debug_assert_eq!(a.len(), b.len(), "dot: length mismatch");
37    a.iter().zip(b.iter()).map(|(ai, bi)| ai * bi).sum()
38}
39
40/// Euclidean 2-norm `sqrt(Σ_i x_i²)` (dimensionless, always `>= 0`).
41///
42/// Computed as `sqrt(dot(x, x))`. For very large magnitudes this can overflow to
43/// `+inf`; inputs are expected to be within normal `f64` range, which holds for
44/// well-scaled linear systems.
45pub fn nrm2(x: &[f64]) -> f64 {
46    dot(x, x).sqrt()
47}
48
49/// AXPY update `y := alpha · x + y`, in place.
50///
51/// `alpha` is a dimensionless scalar; `x` and `y` must have equal length (a
52/// mismatch panics via `debug_assert`). `y` is overwritten with the result.
53pub fn axpy(alpha: f64, x: &[f64], y: &mut [f64]) {
54    debug_assert_eq!(x.len(), y.len(), "axpy: length mismatch");
55    for (yi, xi) in y.iter_mut().zip(x.iter()) {
56        *yi += alpha * xi;
57    }
58}
59
60/// Scale `x := alpha · x`, in place.
61///
62/// `alpha` is a dimensionless scalar. Every element of `x` is multiplied by
63/// `alpha`.
64pub fn scal(alpha: f64, x: &mut [f64]) {
65    for xi in x.iter_mut() {
66        *xi *= alpha;
67    }
68}
69
70#[cfg(test)]
71mod tests {
72    use super::*;
73
74    #[test]
75    fn dot_known() {
76        let a = [1.0, 2.0, 3.0];
77        let b = [4.0, -5.0, 6.0];
78        // 1*4 + 2*-5 + 3*6 = 4 - 10 + 18 = 12
79        assert!((dot(&a, &b) - 12.0).abs() < 1e-12);
80    }
81
82    #[test]
83    fn nrm2_known() {
84        let x = [3.0, 4.0];
85        assert!((nrm2(&x) - 5.0).abs() < 1e-12);
86        assert!(nrm2(&[]).abs() < 1e-12);
87    }
88
89    #[test]
90    fn axpy_known() {
91        let x = [1.0, 1.0, 1.0];
92        let mut y = [10.0, 20.0, 30.0];
93        axpy(2.0, &x, &mut y);
94        assert_eq!(y, [12.0, 22.0, 32.0]);
95    }
96
97    #[test]
98    fn scal_known() {
99        let mut x = [1.0, -2.0, 3.0];
100        scal(-2.0, &mut x);
101        assert_eq!(x, [-2.0, 4.0, -6.0]);
102    }
103}
outram_foam_basic_lib/krylov/vecops.rs

outram_foam_basic_lib/krylov/
vecops.rs
