faer/sparse/
mod.rs

1//! sparse matrix data structures
2//!
3//! most sparse matrix algorithms accept matrices in sparse column-oriented
4//! format. this format represents each column of the matrix by storing the row
5//! indices of its non-zero elements, as well as their values
6//!
7//! the indices and the values are each stored in a contiguous slice (or group
8//! of slices for arbitrary values). in order to specify where each column
9//! starts and ends, a slice of size `ncols + 1` stores the start of each
10//! column, with the last element being equal to the total number of non-zeros
11//! (or the capacity in uncompressed mode)
12//!
13//! # example
14//!
15//! consider the 4-by-5 matrix:
16//! ```notcode
17//! [[10.0, 0.0, 12.0, -1.0, 13.0]
18//!  [ 0.0, 0.0, 25.0, -2.0,  0.0]
19//!  [ 1.0, 0.0,  0.0,  0.0,  0.0]
20//!  [ 4.0, 0.0,  0.0,  0.0,  5.0]]
21//! ```
22//!
23//! the matrix is stored as follows (with `|` marking column boundaries and `,`
24//! separating intra-column values):
25//! ```notcode
26//! column pointers: | 0                | 3 | 3           | 5           | 7           | 9
27//!
28//! row indices    : |    0 ,   2 ,   3 |   |    0 ,    1 |    0 ,    1 |    0 ,    3 |
29//! values         : | 10.0 , 1.0 , 4.0 |   | 12.0 , 25.0 | -1.0 , -2.0 | 13.0 ,  5.0 |
30//! ```
31//!
32//! this is the standard compressed sparse column (csc) layout, see
33//! [wikipedia.org/wiki/Sparse_matrix](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_(CSC_or_CCS))
34mod csc;
35mod csr;
36pub(crate) const NONE: usize = usize::MAX;
37/// sparse linear algebra module.
38/// contains low level routines and the implementation of their corresponding
39/// high level wrappers
40pub mod linalg;
41/// sparse matrix binary and ternary operation implementations
42pub mod ops;
43use crate::internal_prelude_sp::Index;
44pub use csc::{
45	SparseColMat, SparseColMatMut, SparseColMatRef, SymbolicSparseColMat,
46	SymbolicSparseColMatRef, numeric as csc_numeric, symbolic as csc_symbolic,
47};
48pub use csr::{
49	SparseRowMat, SparseRowMatMut, SparseRowMatRef, SymbolicSparseRowMat,
50	SymbolicSparseRowMatRef, numeric as csr_numeric, symbolic as csr_symbolic,
51};
52use reborrow::*;
53extern crate alloc;
54/// pair of indices with `C`-compatible layout
55#[derive(Copy, Clone, Debug, PartialEq, Eq)]
56#[repr(C)]
57pub struct Pair<Row, Col> {
58	/// row index
59	pub row: Row,
60	/// column index
61	pub col: Col,
62}
63/// triplet of indices and value with `C`-compatible layout
64#[derive(Copy, Clone, Debug)]
65#[repr(C)]
66pub struct Triplet<Row, Col, T> {
67	/// row index
68	pub row: Row,
69	/// column index
70	pub col: Col,
71	/// value
72	pub val: T,
73}
74impl<Row, Col> Pair<Row, Col> {
75	/// creates a new pair of indices
76	#[inline]
77	pub const fn new(row: Row, col: Col) -> Self {
78		Pair { row, col }
79	}
80}
81impl<Row, Col, T> Triplet<Row, Col, T> {
82	/// creates a new pair of indices and value
83	#[inline]
84	pub const fn new(row: Row, col: Col, val: T) -> Self {
85		Triplet { row, col, val }
86	}
87}
88/// errors that can occur in sparse algorithms
89#[derive(Copy, Clone, Debug, Eq, PartialEq, Ord, PartialOrd)]
90#[non_exhaustive]
91pub enum FaerError {
92	/// an index exceeding the maximum value (`I::Signed::MAX` for a given index
93	/// type `I`)
94	IndexOverflow,
95	/// memory allocation failed
96	OutOfMemory,
97}
98impl From<dyn_stack::mem::AllocError> for FaerError {
99	#[inline]
100	fn from(value: dyn_stack::mem::AllocError) -> Self {
101		_ = value;
102		FaerError::OutOfMemory
103	}
104}
105impl From<alloc::collections::TryReserveError> for FaerError {
106	#[inline]
107	fn from(value: alloc::collections::TryReserveError) -> Self {
108		_ = value;
109		FaerError::OutOfMemory
110	}
111}
112impl core::fmt::Display for FaerError {
113	#[inline]
114	fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
115		core::fmt::Debug::fmt(self, f)
116	}
117}
118impl core::error::Error for FaerError {}
119/// errors that can occur during the creation of sparse matrices from user input
120#[derive(Copy, Clone, Debug, Eq, PartialEq, Ord, PartialOrd)]
121pub enum CreationError {
122	/// generic error (allocation or index overflow)
123	Generic(FaerError),
124	/// matrix index out-of-bounds error
125	OutOfBounds {
126		/// row of the out-of-bounds index
127		row: usize,
128		/// column of the out-of-bounds index
129		col: usize,
130	},
131}
132impl From<FaerError> for CreationError {
133	#[inline]
134	fn from(value: FaerError) -> Self {
135		Self::Generic(value)
136	}
137}
138impl core::fmt::Display for CreationError {
139	#[inline]
140	fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
141		core::fmt::Debug::fmt(self, f)
142	}
143}
144impl core::error::Error for CreationError {}
145#[inline(always)]
146pub(crate) fn windows2<I>(
147	slice: &[I],
148) -> impl DoubleEndedIterator<Item = &[I; 2]> {
149	slice.windows(2).map(
150		#[inline(always)]
151		|window| unsafe { &*(window.as_ptr() as *const [I; 2]) },
152	)
153}
154#[inline]
155#[track_caller]
156pub(crate) fn try_zeroed<I: bytemuck::Pod>(
157	n: usize,
158) -> Result<alloc::vec::Vec<I>, FaerError> {
159	let mut v = alloc::vec::Vec::new();
160	v.try_reserve_exact(n).map_err(|_| FaerError::OutOfMemory)?;
161	unsafe {
162		core::ptr::write_bytes::<I>(v.as_mut_ptr(), 0u8, n);
163		v.set_len(n);
164	}
165	Ok(v)
166}
167#[inline]
168#[track_caller]
169pub(crate) fn try_collect<I: IntoIterator>(
170	iter: I,
171) -> Result<alloc::vec::Vec<I::Item>, FaerError> {
172	let iter = iter.into_iter();
173	let mut v = alloc::vec::Vec::new();
174	v.try_reserve_exact(iter.size_hint().0)
175		.map_err(|_| FaerError::OutOfMemory)?;
176	v.extend(iter);
177	Ok(v)
178}
179/// the order values should be read in, when constructing/filling from indices
180/// and values
181///
182/// allows separately creating the symbolic structure and filling the numerical
183/// values
184#[derive(Debug, Clone)]
185pub struct Argsort<I> {
186	idx: alloc::vec::Vec<I>,
187	all_nnz: usize,
188	nnz: usize,
189}
190/// algorithmic primitives for sparse matrices
191pub mod utils;
192impl<I: Index, T> core::ops::Index<(usize, usize)>
193	for SparseColMatRef<'_, I, T>
194{
195	type Output = T;
196
197	#[track_caller]
198	fn index(&self, (row, col): (usize, usize)) -> &Self::Output {
199		self.get(row, col).unwrap()
200	}
201}
202impl<I: Index, T> core::ops::Index<(usize, usize)>
203	for SparseRowMatRef<'_, I, T>
204{
205	type Output = T;
206
207	#[track_caller]
208	fn index(&self, (row, col): (usize, usize)) -> &Self::Output {
209		self.get(row, col).unwrap()
210	}
211}
212impl<I: Index, T> core::ops::Index<(usize, usize)>
213	for SparseColMatMut<'_, I, T>
214{
215	type Output = T;
216
217	#[track_caller]
218	fn index(&self, (row, col): (usize, usize)) -> &Self::Output {
219		self.rb().get(row, col).unwrap()
220	}
221}
222impl<I: Index, T> core::ops::Index<(usize, usize)>
223	for SparseRowMatMut<'_, I, T>
224{
225	type Output = T;
226
227	#[track_caller]
228	fn index(&self, (row, col): (usize, usize)) -> &Self::Output {
229		self.rb().get(row, col).unwrap()
230	}
231}
232impl<I: Index, T> core::ops::IndexMut<(usize, usize)>
233	for SparseColMatMut<'_, I, T>
234{
235	#[track_caller]
236	fn index_mut(&mut self, (row, col): (usize, usize)) -> &mut Self::Output {
237		self.rb_mut().get_mut(row, col).unwrap()
238	}
239}
240impl<I: Index, T> core::ops::IndexMut<(usize, usize)>
241	for SparseRowMatMut<'_, I, T>
242{
243	#[track_caller]
244	fn index_mut(&mut self, (row, col): (usize, usize)) -> &mut Self::Output {
245		self.rb_mut().get_mut(row, col).unwrap()
246	}
247}
248impl<I: Index, T> core::ops::Index<(usize, usize)> for SparseColMat<I, T> {
249	type Output = T;
250
251	#[track_caller]
252	fn index(&self, (row, col): (usize, usize)) -> &Self::Output {
253		self.rb().get(row, col).unwrap()
254	}
255}
256impl<I: Index, T> core::ops::Index<(usize, usize)> for SparseRowMat<I, T> {
257	type Output = T;
258
259	#[track_caller]
260	fn index(&self, (row, col): (usize, usize)) -> &Self::Output {
261		self.rb().get(row, col).unwrap()
262	}
263}
264impl<I: Index, T> core::ops::IndexMut<(usize, usize)> for SparseColMat<I, T> {
265	#[track_caller]
266	fn index_mut(&mut self, (row, col): (usize, usize)) -> &mut Self::Output {
267		self.rb_mut().get_mut(row, col).unwrap()
268	}
269}
270impl<I: Index, T> core::ops::IndexMut<(usize, usize)> for SparseRowMat<I, T> {
271	#[track_caller]
272	fn index_mut(&mut self, (row, col): (usize, usize)) -> &mut Self::Output {
273		self.rb_mut().get_mut(row, col).unwrap()
274	}
275}
276#[cfg(test)]
277mod tests {
278	use super::*;
279	use crate::assert;
280	#[test]
281	fn test_from_indices() {
282		let nrows = 5;
283		let ncols = 4;
284		let indices = &[
285			Pair::new(0, 0),
286			Pair::new(1, 2),
287			Pair::new(0, 0),
288			Pair::new(1, 1),
289			Pair::new(0, 1),
290			Pair::new(3, 3),
291			Pair::new(3, 3usize),
292		];
293		let values = &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0f64];
294		let triplets = &[
295			Triplet::new(0, 0, 1.0),
296			Triplet::new(1, 2, 2.0),
297			Triplet::new(0, 0, 3.0),
298			Triplet::new(1, 1, 4.0),
299			Triplet::new(0, 1, 5.0),
300			Triplet::new(3, 3, 6.0),
301			Triplet::new(3, 3usize, 7.0_f64),
302		];
303		{
304			let mat = SymbolicSparseColMat::try_new_from_indices(
305				nrows, ncols, indices,
306			);
307			assert!(mat.is_ok());
308			let (mat, order) = mat.unwrap();
309			assert!(mat.nrows() == nrows);
310			assert!(mat.ncols() == ncols);
311			assert!(mat.col_ptr() == &[0, 1, 3, 4, 5]);
312			assert!(mat.col_nnz() == None);
313			assert!(mat.row_idx() == &[0, 0, 1, 1, 3]);
314			let mat =
315				SparseColMat::<_, f64>::new_from_argsort(mat, &order, values)
316					.unwrap();
317			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
318		}
319		{
320			let mat =
321				SparseColMat::try_new_from_triplets(nrows, ncols, triplets);
322			assert!(mat.is_ok());
323			let mat = mat.unwrap();
324			assert!(mat.nrows() == nrows);
325			assert!(mat.ncols() == ncols);
326			assert!(mat.col_ptr() == &[0, 1, 3, 4, 5]);
327			assert!(mat.col_nnz() == None);
328			assert!(mat.row_idx() == &[0, 0, 1, 1, 3]);
329			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
330		}
331		{
332			let mat = SymbolicSparseRowMat::try_new_from_indices(
333				nrows, ncols, indices,
334			);
335			assert!(mat.is_ok());
336			let (mat, order) = mat.unwrap();
337			assert!(mat.nrows() == nrows);
338			assert!(mat.ncols() == ncols);
339			assert!(mat.row_ptr() == &[0, 2, 4, 4, 5, 5]);
340			assert!(mat.row_nnz() == None);
341			assert!(mat.col_idx() == &[0, 1, 1, 2, 3]);
342			let mat =
343				SparseRowMat::<_, f64>::new_from_argsort(mat, &order, values)
344					.unwrap();
345			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
346		}
347		{
348			let mat =
349				SparseRowMat::try_new_from_triplets(nrows, ncols, triplets);
350			assert!(mat.is_ok());
351			let mat = mat.unwrap();
352			assert!(mat.nrows() == nrows);
353			assert!(mat.ncols() == ncols);
354			assert!(mat.row_ptr() == &[0, 2, 4, 4, 5, 5]);
355			assert!(mat.row_nnz() == None);
356			assert!(mat.col_idx() == &[0, 1, 1, 2, 3]);
357			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
358		}
359	}
360	#[test]
361	fn test_from_nonnegative_indices() {
362		let nrows = 5;
363		let ncols = 4;
364		let indices = &[
365			Pair::new(0, 0),
366			Pair::new(1, 2),
367			Pair::new(0, 0),
368			Pair::new(1, 1),
369			Pair::new(0, 1),
370			Pair::new(-1, 2),
371			Pair::new(-2, 1),
372			Pair::new(-3, -4),
373			Pair::new(3, 3),
374			Pair::new(3, 3isize),
375		];
376		let values = &[
377			1.0,
378			2.0,
379			3.0,
380			4.0,
381			5.0,
382			f64::NAN,
383			f64::NAN,
384			f64::NAN,
385			6.0,
386			7.0f64,
387		];
388		let triplets = &[
389			Triplet::new(0, 0, 1.0),
390			Triplet::new(1, 2, 2.0),
391			Triplet::new(0, 0, 3.0),
392			Triplet::new(1, 1, 4.0),
393			Triplet::new(0, 1, 5.0),
394			Triplet::new(-1, 2, f64::NAN),
395			Triplet::new(-2, 1, f64::NAN),
396			Triplet::new(-3, -4, f64::NAN),
397			Triplet::new(3, 3, 6.0),
398			Triplet::new(3, 3isize, 7.0_f64),
399		];
400		{
401			let mat =
402				SymbolicSparseColMat::<usize>::try_new_from_nonnegative_indices(
403					nrows, ncols, indices,
404				);
405			assert!(mat.is_ok());
406			let (mat, order) = mat.unwrap();
407			assert!(mat.nrows() == nrows);
408			assert!(mat.ncols() == ncols);
409			assert!(mat.col_ptr() == &[0, 1, 3, 4, 5]);
410			assert!(mat.col_nnz() == None);
411			assert!(mat.row_idx() == &[0, 0, 1, 1, 3]);
412			let mat =
413				SparseColMat::<_, f64>::new_from_argsort(mat, &order, values)
414					.unwrap();
415			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
416		}
417		{
418			let mat =
419				SparseColMat::<usize, _>::try_new_from_nonnegative_triplets(
420					nrows, ncols, triplets,
421				);
422			assert!(mat.is_ok());
423			let mat = mat.unwrap();
424			assert!(mat.nrows() == nrows);
425			assert!(mat.ncols() == ncols);
426			assert!(mat.col_ptr() == &[0, 1, 3, 4, 5]);
427			assert!(mat.col_nnz() == None);
428			assert!(mat.row_idx() == &[0, 0, 1, 1, 3]);
429			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
430		}
431		{
432			let mat =
433				SymbolicSparseRowMat::<usize>::try_new_from_nonnegative_indices(
434					nrows, ncols, indices,
435				);
436			assert!(mat.is_ok());
437			let (mat, order) = mat.unwrap();
438			assert!(mat.nrows() == nrows);
439			assert!(mat.ncols() == ncols);
440			assert!(mat.row_ptr() == &[0, 2, 4, 4, 5, 5]);
441			assert!(mat.row_nnz() == None);
442			assert!(mat.col_idx() == &[0, 1, 1, 2, 3]);
443			let mat =
444				SparseRowMat::<_, f64>::new_from_argsort(mat, &order, values)
445					.unwrap();
446			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
447		}
448		{
449			let mat =
450				SparseRowMat::<usize, _>::try_new_from_nonnegative_triplets(
451					nrows, ncols, triplets,
452				);
453			assert!(mat.is_ok());
454			let mat = mat.unwrap();
455			assert!(mat.nrows() == nrows);
456			assert!(mat.ncols() == ncols);
457			assert!(mat.row_ptr() == &[0, 2, 4, 4, 5, 5]);
458			assert!(mat.row_nnz() == None);
459			assert!(mat.col_idx() == &[0, 1, 1, 2, 3]);
460			assert!(mat.val() == &[1.0 + 3.0, 5.0, 4.0, 2.0, 6.0 + 7.0]);
461		}
462	}
463	#[test]
464	fn test_from_indices_oob_row() {
465		let nrows = 5;
466		let ncols = 4;
467		let indices = &[
468			Pair::new(0, 0),
469			Pair::new(1, 2),
470			Pair::new(0, 0),
471			Pair::new(1, 1),
472			Pair::new(0, 1),
473			Pair::new(3, 3),
474			Pair::new(3, 3),
475			Pair::new(5, 3usize),
476		];
477		let err =
478			SymbolicSparseColMat::try_new_from_indices(nrows, ncols, indices);
479		assert!(err.is_err());
480		let err = err.unwrap_err();
481		assert!(err == CreationError::OutOfBounds { row: 5, col: 3 });
482	}
483	#[test]
484	fn test_from_indices_oob_col() {
485		let nrows = 5;
486		let ncols = 4;
487		let indices = &[
488			Pair::new(0, 0),
489			Pair::new(1, 2),
490			Pair::new(0, 0),
491			Pair::new(1, 1),
492			Pair::new(0, 1),
493			Pair::new(3, 3),
494			Pair::new(3, 3),
495			Pair::new(2, 4usize),
496		];
497		let err =
498			SymbolicSparseColMat::try_new_from_indices(nrows, ncols, indices);
499		assert!(err.is_err());
500		let err = err.unwrap_err();
501		assert!(err == CreationError::OutOfBounds { row: 2, col: 4 });
502	}
503}
faer/sparse/mod.rs

faer/sparse/
mod.rs
