sciforge_lib/maths/sparse/
ops.rs

1use super::csr::SparseMatrix;
2
3pub fn sparse_kronecker(a: &SparseMatrix, b: &SparseMatrix) -> SparseMatrix {
4    let mut triplets = Vec::new();
5    for ia in 0..a.rows {
6        for ka in a.row_ptr[ia]..a.row_ptr[ia + 1] {
7            let ja = a.col_idx[ka];
8            let va = a.values[ka];
9            for ib in 0..b.rows {
10                for kb in b.row_ptr[ib]..b.row_ptr[ib + 1] {
11                    let jb = b.col_idx[kb];
12                    let vb = b.values[kb];
13                    triplets.push((ia * b.rows + ib, ja * b.cols + jb, va * vb));
14                }
15            }
16        }
17    }
18    SparseMatrix::from_triplets(a.rows * b.rows, a.cols * b.cols, &mut triplets)
19}
20
21pub fn sparse_from_diagonals(n: usize, diags: &[(i64, Vec<f64>)]) -> SparseMatrix {
22    let mut triplets = Vec::new();
23    for &(offset, ref vals) in diags {
24        for (k, &v) in vals.iter().enumerate() {
25            let i = if offset >= 0 {
26                k
27            } else {
28                ((-offset) as usize) + k
29            };
30            let j = if offset >= 0 { k + offset as usize } else { k };
31            if i < n && j < n {
32                triplets.push((i, j, v));
33            }
34        }
35    }
36    SparseMatrix::from_triplets(n, n, &mut triplets)
37}
38
39pub fn sparse_add(a: &SparseMatrix, b: &SparseMatrix) -> SparseMatrix {
40    let mut triplets = Vec::new();
41    for i in 0..a.rows {
42        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
43            triplets.push((i, a.col_idx[k], a.values[k]));
44        }
45    }
46    for i in 0..b.rows {
47        for k in b.row_ptr[i]..b.row_ptr[i + 1] {
48            triplets.push((i, b.col_idx[k], b.values[k]));
49        }
50    }
51    SparseMatrix::from_triplets(a.rows, a.cols, &mut triplets)
52}
53
54pub fn sparse_scale(a: &SparseMatrix, scalar: f64) -> SparseMatrix {
55    let mut result = a.clone();
56    for v in &mut result.values {
57        *v *= scalar;
58    }
59    result
60}
61
62pub fn sparse_mat_mul(a: &SparseMatrix, b: &SparseMatrix) -> SparseMatrix {
63    let mut triplets = Vec::new();
64    for i in 0..a.rows {
65        for ka in a.row_ptr[i]..a.row_ptr[i + 1] {
66            let ja = a.col_idx[ka];
67            let va = a.values[ka];
68            for kb in b.row_ptr[ja]..b.row_ptr[ja + 1] {
69                triplets.push((i, b.col_idx[kb], va * b.values[kb]));
70            }
71        }
72    }
73    SparseMatrix::from_triplets(a.rows, b.cols, &mut triplets)
74}
75
76pub fn sparse_trace(a: &SparseMatrix) -> f64 {
77    let mut trace = 0.0;
78    for i in 0..a.rows.min(a.cols) {
79        trace += a.get(i, i);
80    }
81    trace
82}
83
84pub fn sparse_frobenius_norm(a: &SparseMatrix) -> f64 {
85    a.values.iter().map(|v| v * v).sum::<f64>().sqrt()
86}
87
88pub fn sparse_infinity_norm(a: &SparseMatrix) -> f64 {
89    let mut max_row_sum = 0.0_f64;
90    for i in 0..a.rows {
91        let row_sum: f64 = (a.row_ptr[i]..a.row_ptr[i + 1])
92            .map(|k| a.values[k].abs())
93            .sum();
94        max_row_sum = max_row_sum.max(row_sum);
95    }
96    max_row_sum
97}
98
99pub fn sparse_one_norm(a: &SparseMatrix) -> f64 {
100    let mut col_sums = vec![0.0; a.cols];
101    for i in 0..a.rows {
102        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
103            col_sums[a.col_idx[k]] += a.values[k].abs();
104        }
105    }
106    col_sums.iter().cloned().fold(0.0, f64::max)
107}
108
109pub fn sparse_diagonal(a: &SparseMatrix) -> Vec<f64> {
110    let n = a.rows.min(a.cols);
111    (0..n).map(|i| a.get(i, i)).collect()
112}
113
114pub fn sparse_lower_triangular(a: &SparseMatrix) -> SparseMatrix {
115    let mut triplets = Vec::new();
116    for i in 0..a.rows {
117        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
118            if a.col_idx[k] <= i {
119                triplets.push((i, a.col_idx[k], a.values[k]));
120            }
121        }
122    }
123    SparseMatrix::from_triplets(a.rows, a.cols, &mut triplets)
124}
125
126pub fn sparse_upper_triangular(a: &SparseMatrix) -> SparseMatrix {
127    let mut triplets = Vec::new();
128    for i in 0..a.rows {
129        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
130            if a.col_idx[k] >= i {
131                triplets.push((i, a.col_idx[k], a.values[k]));
132            }
133        }
134    }
135    SparseMatrix::from_triplets(a.rows, a.cols, &mut triplets)
136}
137
138pub fn sparse_row_sum(a: &SparseMatrix) -> Vec<f64> {
139    (0..a.rows)
140        .map(|i| (a.row_ptr[i]..a.row_ptr[i + 1]).map(|k| a.values[k]).sum())
141        .collect()
142}
143
144pub fn sparse_col_sum(a: &SparseMatrix) -> Vec<f64> {
145    let mut sums = vec![0.0; a.cols];
146    for i in 0..a.rows {
147        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
148            sums[a.col_idx[k]] += a.values[k];
149        }
150    }
151    sums
152}
153
154pub fn sparse_is_symmetric(a: &SparseMatrix, tol: f64) -> bool {
155    if a.rows != a.cols {
156        return false;
157    }
158    for i in 0..a.rows {
159        for k in a.row_ptr[i]..a.row_ptr[i + 1] {
160            let j = a.col_idx[k];
161            if (a.values[k] - a.get(j, i)).abs() > tol {
162                return false;
163            }
164        }
165    }
166    true
167}
sciforge_lib/maths/sparse/ops.rs

sciforge_lib/maths/sparse/
ops.rs
