peroxide/structure/
sparse.rs1use crate::structure::matrix::Matrix;
6use crate::traits::math::LinearOp;
7#[allow(unused_imports)]
8use crate::traits::matrix::{Form, LinearAlgebra, SolveKind, PQLU, QR, SVD, UPLO, WAZD};
9use crate::util::non_macro::zeros;
11use std::ops::Mul;
12
13#[derive(Debug, Clone)]
14pub struct SPMatrix {
15 pub row: usize,
16 pub col: usize,
17 pub nnz: usize,
18 pub col_ptr: Vec<usize>,
19 pub row_ics: Vec<usize>,
20 pub data: Vec<f64>,
21}
22
23impl SPMatrix {
24 pub fn new(row: usize, col: usize, nnz: usize) -> Self {
25 SPMatrix {
26 row,
27 col,
28 nnz,
29 col_ptr: vec![0usize; col + 1],
30 row_ics: vec![0usize, nnz],
31 data: vec![0f64; nnz],
32 }
33 }
34
35 pub fn from_dense(m: &Matrix) -> Self {
36 let mut data: Vec<f64> = Vec::new();
37 let mut row_ics: Vec<usize> = Vec::new();
38 let mut col_ptr: Vec<usize> = vec![0usize; m.col + 1];
39 let mut k = 0usize;
40 for j in 0..m.col {
41 for i in 0..m.row {
42 let val = m[(i, j)];
43 if val != 0f64 {
44 data.push(val);
45 row_ics.push(i);
46 k += 1;
47 }
48 }
49 col_ptr[j + 1] = k;
50 }
51
52 SPMatrix {
53 row: m.row,
54 col: m.col,
55 nnz: data.len(),
56 col_ptr,
57 row_ics,
58 data,
59 }
60 }
61
62 pub fn to_dense(&self) -> Matrix {
63 let mut m = zeros(self.row, self.col);
64
65 for j in 0..self.col {
66 for i in self.col_ptr[j]..self.col_ptr[j + 1] {
67 let k = self.row_ics[i];
68 m[(k, j)] = self.data[i];
69 }
70 }
71 m
72 }
73
74 pub fn col_ptr(&self) -> &Vec<usize> {
75 &self.col_ptr
76 }
77
78 pub fn row_ics(&self) -> &Vec<usize> {
79 &self.row_ics
80 }
81
82 pub fn data(&self) -> &Vec<f64> {
83 &self.data
84 }
85
86 pub fn transpose(&self) -> Self {
87 let row = self.row;
88 let col = self.col;
89 let nnz = self.nnz;
90 let col_ptr = self.col_ptr();
91 let row_ics = self.row_ics();
92 let data = self.data();
93 let mut count = vec![0usize; row];
94 let mut result = Self::new(col, row, nnz);
95
96 for i in 0..col {
97 for j in col_ptr[i]..col_ptr[i + 1] {
98 let k = row_ics[j];
99 count[k] += 1;
100 }
101 }
102 for j in 0..row {
103 result.col_ptr[j + 1] = result.col_ptr[j] + count[j];
104 count[j] = 0;
105 }
106 for i in 0..col {
107 for j in col_ptr[i]..col_ptr[i + 1] {
108 let k = row_ics[j];
109 let index = result.col_ptr[k] + count[k];
110 result.row_ics[index] = i;
111 result.data[index] = data[j];
112 count[k] += 1;
113 }
114 }
115 result
116 }
117
118 pub fn t(&self) -> Self {
119 self.transpose()
120 }
121}
122
123impl LinearOp<Vec<f64>, Vec<f64>> for SPMatrix {
124 fn apply(&self, rhs: &Vec<f64>) -> Vec<f64> {
125 let mut y = vec![0f64; self.row];
126 let col_ptr = self.col_ptr();
127 let row_ics = self.row_ics();
128 let data = self.data();
129 for j in 0..self.col {
130 for i in col_ptr[j]..col_ptr[j + 1] {
131 y[row_ics[i]] += data[i] * rhs[j];
132 }
133 }
134 y
135 }
136}
137
138impl LinearAlgebra<Matrix> for SPMatrix {
142 fn back_subs(&self, _b: &[f64]) -> Vec<f64> {
143 unimplemented!()
144 }
145
146 fn forward_subs(&self, _b: &[f64]) -> Vec<f64> {
147 unimplemented!()
148 }
149
150 fn lu(&self) -> PQLU<Matrix> {
151 self.to_dense().lu()
152 }
153
154 fn waz(&self, _d_form: Form) -> Option<WAZD<Matrix>> {
155 unimplemented!()
156 }
157
158 fn qr(&self) -> QR<Matrix> {
159 self.to_dense().qr()
160 }
161
162 fn svd(&self) -> SVD<Matrix> {
163 unimplemented!()
164 }
165
166 #[cfg(feature = "O3")]
167 fn cholesky(&self, _uplo: UPLO) -> Matrix {
168 unimplemented!()
169 }
170
171 fn rref(&self) -> Matrix {
172 self.to_dense().rref()
173 }
174
175 fn det(&self) -> f64 {
176 self.to_dense().det()
177 }
178
179 fn block(&self) -> (Matrix, Matrix, Matrix, Matrix) {
180 self.to_dense().block()
181 }
182
183 fn inv(&self) -> Matrix {
184 self.to_dense().inv()
185 }
186
187 fn pseudo_inv(&self) -> Matrix {
188 self.to_dense().pseudo_inv()
189 }
190
191 fn solve(&self, _b: &[f64], _sk: SolveKind) -> Vec<f64> {
192 unimplemented!()
193 }
194
195 fn solve_mat(&self, _m: &Matrix, _sk: SolveKind) -> Matrix {
196 unimplemented!()
197 }
198
199 fn is_symmetric(&self) -> bool {
200 unimplemented!()
201 }
202}
203
204impl Mul<Vec<f64>> for SPMatrix {
206 type Output = Vec<f64>;
207 fn mul(self, rhs: Vec<f64>) -> Self::Output {
208 self.apply(&rhs)
209 }
210}
211
212impl<'a, 'b> Mul<&'b Vec<f64>> for &'a SPMatrix {
214 type Output = Vec<f64>;
215 fn mul(self, rhs: &'b Vec<f64>) -> Self::Output {
216 self.apply(rhs)
217 }
218}
219
220impl Into<Matrix> for SPMatrix {
221 fn into(self) -> Matrix {
222 self.to_dense()
223 }
224}
225
226impl Into<SPMatrix> for Matrix {
227 fn into(self) -> SPMatrix {
228 SPMatrix::from_dense(&self)
229 }
230}