oxiphysics_core/sparse/
functions_2.rs

1//! Auto-generated module
2//!
3//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
4
5#![allow(clippy::needless_range_loop)]
6#[allow(unused_imports)]
7use super::functions::*;
8use super::types::CsrMatrix;
9
10/// Extract the upper triangular part of a sparse matrix (including diagonal).
11#[allow(dead_code)]
12pub fn upper_triangular(a: &CsrMatrix) -> CsrMatrix {
13    let mut rows = Vec::new();
14    let mut cols = Vec::new();
15    let mut vals = Vec::new();
16    for i in 0..a.nrows {
17        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
18            let j = a.col_idx[k];
19            if j >= i {
20                rows.push(i);
21                cols.push(j);
22                vals.push(a.values[k]);
23            }
24        }
25    }
26    CsrMatrix::from_triplets(a.nrows, a.ncols, &rows, &cols, &vals)
27}
28/// Symmetrise a sparse matrix: A_sym = (A + A^T) / 2.
29#[allow(dead_code)]
30pub fn symmetrise_sparse(a: &CsrMatrix) -> CsrMatrix {
31    let at = a.transpose();
32    let sum = a.add(&at);
33    sum.scale(0.5)
34}
35/// Compute the 1-norm (max column sum of absolute values) of a sparse matrix.
36#[allow(dead_code)]
37pub fn sparse_one_norm(a: &CsrMatrix) -> f64 {
38    let mut col_sums = vec![0.0f64; a.ncols];
39    for k in 0..a.values.len() {
40        col_sums[a.col_idx[k]] += a.values[k].abs();
41    }
42    col_sums.iter().cloned().fold(0.0f64, f64::max)
43}
44/// Compute the infinity-norm (max row sum of absolute values) of a sparse matrix.
45#[allow(dead_code)]
46pub fn sparse_inf_norm(a: &CsrMatrix) -> f64 {
47    (0..a.nrows)
48        .map(|i| {
49            (a.row_ptr[i]..a.row_ptr[i + 1])
50                .map(|k| a.values[k].abs())
51                .sum::<f64>()
52        })
53        .fold(0.0f64, f64::max)
54}
55/// Compute the Frobenius norm of a sparse matrix.
56#[allow(dead_code)]
57pub fn sparse_frobenius_norm(a: &CsrMatrix) -> f64 {
58    a.values.iter().map(|x| x * x).sum::<f64>().sqrt()
59}
60/// Set the diagonal of a CSR matrix (in-place). Only existing diagonal entries are modified.
61#[allow(dead_code)]
62pub fn set_diagonal(a: &mut CsrMatrix, d: &[f64]) {
63    assert_eq!(d.len(), a.nrows.min(a.ncols));
64    for i in 0..d.len() {
65        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
66            if a.col_idx[k] == i {
67                a.values[k] = d[i];
68                break;
69            }
70        }
71    }
72}
73/// Add a scalar multiple of the identity to a CSR matrix: A + alpha * I.
74/// Requires the diagonal entries to already exist in the sparsity pattern.
75#[allow(dead_code)]
76pub fn add_identity_scaled(a: &CsrMatrix, alpha: f64) -> CsrMatrix {
77    let n = a.nrows.min(a.ncols);
78    let mut new_vals = a.values.clone();
79    for i in 0..n {
80        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
81            if a.col_idx[k] == i {
82                new_vals[k] += alpha;
83                break;
84            }
85        }
86    }
87    CsrMatrix {
88        nrows: a.nrows,
89        ncols: a.ncols,
90        row_ptr: a.row_ptr.clone(),
91        col_idx: a.col_idx.clone(),
92        values: new_vals,
93    }
94}
95#[cfg(test)]
96mod tests_extended {
97    use super::*;
98    use crate::sparse::BlockJacobi;
99    fn approx(a: f64, b: f64, tol: f64) -> bool {
100        (a - b).abs() < tol
101    }
102    #[test]
103    fn test_bicgstab_identity_system() {
104        let a = identity_csr(5);
105        let b = vec![1.0, 2.0, 3.0, 4.0, 5.0];
106        let (x, conv, _) = bicgstab_solve(&a, &b, None, 1e-10, 100);
107        assert!(conv, "should converge");
108        for i in 0..5 {
109            assert!(approx(x[i], b[i], 1e-6));
110        }
111    }
112    #[test]
113    fn test_bicgstab_laplacian() {
114        let a = laplacian_1d(5, 1.0);
115        let x_true = vec![1.0, 2.0, 3.0, 2.0, 1.0];
116        let b = a.matvec(&x_true);
117        let (x, conv, _) = bicgstab_solve(&a, &b, None, 1e-8, 200);
118        assert!(conv, "BICGStab should converge on 1D Laplacian");
119        for i in 0..5 {
120            assert!(approx(x[i], x_true[i], 1e-4), "x[{i}]={}", x[i]);
121        }
122    }
123    #[test]
124    fn test_bicgstab_zero_rhs() {
125        let a = laplacian_1d(4, 1.0);
126        let b = vec![0.0; 4];
127        let (x, conv, _) = bicgstab_solve(&a, &b, None, 1e-10, 50);
128        assert!(conv);
129        for i in 0..4 {
130            assert!(approx(x[i], 0.0, 1e-10));
131        }
132    }
133    #[test]
134    fn test_bicgstab_initial_guess() {
135        let a = identity_csr(3);
136        let b = vec![4.0, 5.0, 6.0];
137        let x0 = vec![3.9, 4.9, 5.9];
138        let (x, conv, iters) = bicgstab_solve(&a, &b, Some(&x0), 1e-10, 50);
139        assert!(conv, "should converge from near-exact initial guess");
140        for i in 0..3 {
141            assert!(approx(x[i], b[i], 1e-6));
142        }
143        assert!(iters < 10, "iters={iters}");
144    }
145    #[test]
146    fn test_gauss_seidel_diagonal_system() {
147        let a = identity_csr(4);
148        let b = vec![1.0, 2.0, 3.0, 4.0];
149        let (x, iters) = gauss_seidel_solve(&a, &b, 1e-10, 100);
150        assert!(iters <= 3, "iters={iters}");
151        for i in 0..4 {
152            assert!(approx(x[i], b[i], 1e-10));
153        }
154    }
155    #[test]
156    fn test_gauss_seidel_converges_laplacian() {
157        let a = laplacian_1d(5, 1.0);
158        let a2 = add_identity_scaled(&a, 4.0);
159        let x_true = vec![1.0, 2.0, 3.0, 4.0, 5.0];
160        let b = a2.matvec(&x_true);
161        let (x, _iters) = gauss_seidel_solve(&a2, &b, 1e-8, 500);
162        for i in 0..5 {
163            assert!(approx(x[i], x_true[i], 1e-4), "x[{i}]={}", x[i]);
164        }
165    }
166    #[test]
167    fn test_forward_substitution_identity() {
168        let l = identity_csr(4);
169        let b = vec![1.0, 2.0, 3.0, 4.0];
170        let x = forward_substitution(&l, &b);
171        for i in 0..4 {
172            assert!(approx(x[i], b[i], 1e-12));
173        }
174    }
175    #[test]
176    fn test_backward_substitution_identity() {
177        let u = identity_csr(4);
178        let b = vec![4.0, 3.0, 2.0, 1.0];
179        let x = backward_substitution(&u, &b);
180        for i in 0..4 {
181            assert!(approx(x[i], b[i], 1e-12));
182        }
183    }
184    #[test]
185    fn test_forward_substitution_lower_triangular() {
186        let l = CsrMatrix::from_triplets(2, 2, &[0, 1, 1], &[0, 0, 1], &[2.0, 3.0, 4.0]);
187        let b = vec![4.0, 11.0];
188        let x = forward_substitution(&l, &b);
189        assert!(approx(x[0], 2.0, 1e-10));
190        assert!(approx(x[1], 1.25, 1e-10));
191    }
192    #[test]
193    fn test_backward_substitution_upper_triangular() {
194        let u = CsrMatrix::from_triplets(2, 2, &[0, 0, 1], &[0, 1, 1], &[2.0, 3.0, 4.0]);
195        let b = vec![11.0, 8.0];
196        let x = backward_substitution(&u, &b);
197        assert!(approx(x[1], 2.0, 1e-10));
198        assert!(approx(x[0], 2.5, 1e-10));
199    }
200    #[test]
201    fn test_banded_matrix_tridiagonal() {
202        let a = banded_matrix(4, &[-1, 0, 1], &[-1.0, 2.0, -1.0]);
203        assert!(approx(a.get(0, 0), 2.0, 1e-12));
204        assert!(approx(a.get(0, 1), -1.0, 1e-12));
205        assert!(approx(a.get(1, 0), -1.0, 1e-12));
206        assert!(approx(a.get(3, 3), 2.0, 1e-12));
207        assert!(approx(a.get(0, 3), 0.0, 1e-12));
208    }
209    #[test]
210    fn test_banded_matrix_diagonal_only() {
211        let a = banded_matrix(3, &[0], &[5.0]);
212        for i in 0..3 {
213            assert!(approx(a.get(i, i), 5.0, 1e-12));
214        }
215        assert!(approx(a.get(0, 1), 0.0, 1e-12));
216    }
217    #[test]
218    fn test_row_scale_identity() {
219        let a = identity_csr(3);
220        let s = vec![2.0, 3.0, 5.0];
221        let b = row_scale(&a, &s);
222        for i in 0..3 {
223            assert!(approx(b.get(i, i), s[i], 1e-12));
224        }
225    }
226    #[test]
227    fn test_col_scale_identity() {
228        let a = identity_csr(3);
229        let s = vec![7.0, 11.0, 13.0];
230        let b = col_scale(&a, &s);
231        for i in 0..3 {
232            assert!(approx(b.get(i, i), s[i], 1e-12));
233        }
234    }
235    #[test]
236    fn test_equilibration_scales_identity() {
237        let a = identity_csr(4);
238        let (rs, cs) = equilibration_scales(&a);
239        for i in 0..4 {
240            assert!(approx(rs[i], 1.0, 1e-10));
241        }
242        for j in 0..4 {
243            assert!(approx(cs[j], 1.0, 1e-10));
244        }
245    }
246    #[test]
247    fn test_multi_rhs_cg_identity() {
248        let a = identity_csr(3);
249        let b1 = vec![1.0, 2.0, 3.0];
250        let b2 = vec![4.0, 5.0, 6.0];
251        let xs = multi_rhs_cg(&a, &[b1.clone(), b2.clone()], 1e-10, 50);
252        for i in 0..3 {
253            assert!(approx(xs[0][i], b1[i], 1e-6));
254        }
255        for i in 0..3 {
256            assert!(approx(xs[1][i], b2[i], 1e-6));
257        }
258    }
259    #[test]
260    fn test_lower_triangular_extraction() {
261        let a = laplacian_1d(4, 1.0);
262        let l = lower_triangular(&a);
263        for i in 0..4 {
264            for j in (i + 1)..4 {
265                assert!(
266                    approx(l.get(i, j), 0.0, 1e-12),
267                    "L[{i},{j}]={}",
268                    l.get(i, j)
269                );
270            }
271        }
272    }
273    #[test]
274    fn test_upper_triangular_extraction() {
275        let a = laplacian_1d(4, 1.0);
276        let u = upper_triangular(&a);
277        for i in 1..4 {
278            for j in 0..i {
279                assert!(
280                    approx(u.get(i, j), 0.0, 1e-12),
281                    "U[{i},{j}]={}",
282                    u.get(i, j)
283                );
284            }
285        }
286    }
287    #[test]
288    fn test_symmetrise_sparse_already_symmetric() {
289        let a = laplacian_1d(4, 1.0);
290        let asym = symmetrise_sparse(&a);
291        for i in 0..4 {
292            for j in 0..4 {
293                assert!(approx(asym.get(i, j), a.get(i, j), 1e-10));
294            }
295        }
296    }
297    #[test]
298    fn test_symmetrise_sparse_asymmetric_input() {
299        let a = CsrMatrix::from_triplets(2, 2, &[0, 0, 1], &[0, 1, 0], &[1.0, 2.0, 0.0]);
300        let asym = symmetrise_sparse(&a);
301        assert!(approx(asym.get(0, 1), 1.0, 1e-10));
302        assert!(approx(asym.get(1, 0), 1.0, 1e-10));
303    }
304    #[test]
305    fn test_sparse_one_norm_identity() {
306        let a = identity_csr(5);
307        assert!(approx(sparse_one_norm(&a), 1.0, 1e-12));
308    }
309    #[test]
310    fn test_sparse_inf_norm_identity() {
311        let a = identity_csr(5);
312        assert!(approx(sparse_inf_norm(&a), 1.0, 1e-12));
313    }
314    #[test]
315    fn test_sparse_frobenius_norm_identity() {
316        let a = identity_csr(5);
317        assert!(approx(sparse_frobenius_norm(&a), (5.0f64).sqrt(), 1e-12));
318    }
319    #[test]
320    fn test_sparse_one_norm_laplacian() {
321        let a = laplacian_1d(4, 1.0);
322        assert!(approx(sparse_one_norm(&a), 4.0, 1e-10));
323    }
324    #[test]
325    fn test_add_identity_scaled_zero() {
326        let a = laplacian_1d(3, 1.0);
327        let b = add_identity_scaled(&a, 0.0);
328        for i in 0..3 {
329            for j in 0..3 {
330                assert!(approx(b.get(i, j), a.get(i, j), 1e-12));
331            }
332        }
333    }
334    #[test]
335    fn test_add_identity_scaled_positive() {
336        let a = identity_csr(3);
337        let b = add_identity_scaled(&a, 2.0);
338        for i in 0..3 {
339            assert!(approx(b.get(i, i), 3.0, 1e-12));
340        }
341    }
342    #[test]
343    fn test_block_jacobi_identity() {
344        let a = identity_csr(4);
345        let bj = BlockJacobi::new(&a, 2);
346        let b = vec![1.0, 2.0, 3.0, 4.0];
347        let x = bj.apply(&b);
348        for i in 0..4 {
349            assert!(approx(x[i], b[i], 1e-10));
350        }
351    }
352    #[test]
353    fn test_block_jacobi_diagonal_matrix() {
354        let a = CsrMatrix::from_triplets(4, 4, &[0, 1, 2, 3], &[0, 1, 2, 3], &[2.0, 4.0, 3.0, 6.0]);
355        let bj = BlockJacobi::new(&a, 2);
356        let b = vec![2.0, 4.0, 3.0, 6.0];
357        let x = bj.apply(&b);
358        for i in 0..4 {
359            assert!(approx(x[i], 1.0, 1e-10));
360        }
361    }
362    #[test]
363    fn test_lanczos_identity_eigenvalue() {
364        let a = identity_csr(5);
365        let b0 = vec![1.0, 0.0, 0.0, 0.0, 0.0];
366        let (alphas, _betas, _q) = lanczos(&a, &b0, 3);
367        assert!(approx(alphas[0], 1.0, 1e-10));
368    }
369    #[test]
370    fn test_lanczos_returns_k_alphas() {
371        let a = laplacian_1d(6, 1.0);
372        let b0 = vec![1.0; 6];
373        let (alphas, betas, q) = lanczos(&a, &b0, 4);
374        assert!(alphas.len() <= 4);
375        assert!(betas.len() <= 4);
376        assert!(!q.is_empty());
377    }
378    #[test]
379    fn test_lanczos_vectors_orthogonal() {
380        let a = laplacian_1d(6, 1.0);
381        let b0: Vec<f64> = (0..6).map(|i| (i + 1) as f64).collect();
382        let (_alphas, _betas, q) = lanczos(&a, &b0, 4);
383        for i in 0..q.len() {
384            for j in (i + 1)..q.len() {
385                let dot: f64 = q[i].iter().zip(q[j].iter()).map(|(a, b)| a * b).sum();
386                assert!(dot.abs() < 1e-8, "q[{i}]·q[{j}] = {dot}");
387            }
388        }
389    }
390    #[test]
391    fn test_ritz_values_1x1() {
392        let alphas = vec![3.125];
393        let betas: Vec<f64> = vec![];
394        let rv = ritz_values(&alphas, &betas);
395        assert_eq!(rv.len(), 1);
396        assert!(approx(rv[0], 3.125, 1e-6));
397    }
398    #[test]
399    fn test_ritz_values_count() {
400        let alphas = vec![2.0, 5.0, 7.0];
401        let betas = vec![0.5, 0.3];
402        let rv = ritz_values(&alphas, &betas);
403        assert_eq!(rv.len(), 3, "should return 3 Ritz values");
404        for v in &rv {
405            assert!(v.is_finite(), "Ritz value should be finite: {v}");
406        }
407    }
408    #[test]
409    fn test_minres_identity() {
410        let a = identity_csr(4);
411        let b = vec![1.0, 2.0, 3.0, 4.0];
412        let (x, conv, _) = minres_solve(&a, &b, 1e-8, 100);
413        assert!(
414            conv || {
415                let ax = a.matvec(&x);
416                let res: Vec<f64> = b.iter().zip(ax.iter()).map(|(bi, ai)| bi - ai).collect();
417                res.iter().map(|x| x * x).sum::<f64>().sqrt() < 1e-4
418            }
419        );
420    }
421    #[test]
422    fn test_minres_zero_rhs() {
423        let a = identity_csr(3);
424        let b = vec![0.0, 0.0, 0.0];
425        let (x, _conv, _) = minres_solve(&a, &b, 1e-10, 50);
426        for i in 0..3 {
427            assert!(x[i].abs() < 1e-10);
428        }
429    }
430}
oxiphysics_core/sparse/functions_2.rs

oxiphysics_core/sparse/
functions_2.rs
