peroxide/structure/
sparse.rs

1//! Sparse matrix (CCS format)
2//!
3//! * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007.
4
5use 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};
9//use crate::traits::math::{InnerProduct, LinearOp, Norm, Normed, Vector};
10use 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
138/// Linear algebra for sparse matrix
139///
140/// **Caution** : In every ops in this trait, there is converting process to dense matrix
141impl 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
204/// Matrix multiplication with vector
205impl 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
212/// Reference version of matrix multiplication with vector
213impl<'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}
peroxide/structure/sparse.rs

peroxide/structure/
sparse.rs
