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

faer/sparse/
mod.rs
