faer/sparse/
mod.rs

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

faer/sparse/
mod.rs
