nalgebra-sparse-linalg 0.1.10

Sparse linear algebra library for Rust using nalgebra including linear solvers and SVD.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

use nalgebra_sparse::{na::RealField, CooMatrix};
use super::{super::*, coarsen::Mark};

pub(crate) fn build_p<N>(
    a: &CsrMatrix<N>,
    marks: &[Mark],
    coarse_of: &[usize],
    s: &[Vec<usize>],
) -> CsrMatrix<N>
where
    N: RealField + Copy
{
    let n = a.nrows();
    let n_coarse = coarse_of.iter().filter(|&&c| c != usize::MAX).count();

    // Better capacity estimation based on connectivity
    let f_points = marks.iter().filter(|&&m| matches!(m, Mark::F)).count();
    let estimated_nnz = n_coarse + f_points * 4; // More accurate estimate
    
    let mut trip = CooMatrix::new(n, n_coarse);
    trip.reserve(estimated_nnz);

    // Pre-allocate reusable vectors
    let mut c_neighbors = Vec::with_capacity(8);
    let mut weights = Vec::with_capacity(8);

    for i in 0..n {
        match marks[i] {
            Mark::C => {
                trip.push(i, coarse_of[i], N::one());
            }
            Mark::F => {
                c_neighbors.clear();
                // Collect coarse neighbors
                for &nbr in &s[i] {
                    if matches!(marks[nbr], Mark::C) {
                        c_neighbors.push(nbr);
                    }
                }

                // If no coarse neighbors, add a zero entry
                if c_neighbors.is_empty() {
                    if n_coarse > 0 {
                        trip.push(i, 0, N::one());
                    }
                    continue;
                }

                // Early exit if no diagonal entry
                let Some(diag_entry) = a.get_entry(i, i) else {
                    let equal_weight = N::one() / N::from_usize(c_neighbors.len()).unwrap();
                    for &nbr in &c_neighbors {
                        trip.push(i, coarse_of[nbr], equal_weight);
                    }
                    continue;
                };

                let diag = diag_entry.into_value();
                let mut weight_sum = N::zero();
                
                weights.clear();
                weights.reserve(c_neighbors.len());
                
                // Compute weights in single pass
                for &nbr in &c_neighbors {
                    if let Some(a_ij) = a.get_entry(i, nbr) {
                        let w = -a_ij.into_value() / diag;
                        weights.push((coarse_of[nbr], w));
                        weight_sum += w;
                    }
                }
                
                // Add entries with normalized weights
                if weight_sum != N::zero() {
                    let inv_sum = N::one() / weight_sum;
                    for (col, w) in weights.drain(..) {
                        trip.push(i, col, w * inv_sum);
                    }
                } else {
                    let equal_weight = N::one() / N::from_usize(c_neighbors.len()).unwrap();
                    for &nbr in &c_neighbors {
                        trip.push(i, coarse_of[nbr], equal_weight);
                    }
                }
            }
            Mark::Unmarked => unreachable!("All points should be marked during coarsening"),
        }
    }
    
    CsrMatrix::from(&trip)
}