1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
//! Sparse matrices and algorithms for [nalgebra](https://www.nalgebra.org).
//!
//! This crate extends `nalgebra` with sparse matrix formats and operations on sparse matrices.
//!
//! ## Goals
//! The long-term goals for this crate are listed below.
//!
//! - Provide proven sparse matrix formats in an easy-to-use and idiomatic Rust API that
//!   naturally integrates with `nalgebra`.
//! - Provide additional expert-level APIs for fine-grained control over operations.
//! - Integrate well with external sparse matrix libraries.
//! - Provide native Rust high-performance routines, including parallel matrix operations.
//!
//! ## Highlighted current features
//!
//! - [CSR](csr::CsrMatrix), [CSC](csc::CscMatrix) and [COO](coo::CooMatrix) formats, and
//!   [conversions](`convert`) between them.
//! - Common arithmetic operations are implemented. See the [`ops`] module.
//! - Sparsity patterns in CSR and CSC matrices are explicitly represented by the
//!   [SparsityPattern](pattern::SparsityPattern) type, which encodes the invariants of the
//!   associated index data structures.
//! - [proptest strategies](`proptest`) for sparse matrices when the feature
//!   `proptest-support` is enabled.
//! - [matrixcompare support](https://crates.io/crates/matrixcompare) for effortless
//!   (approximate) comparison of matrices in test code (requires the `compare` feature).
//!
//! ## Current state
//!
//! The library is in an early, but usable state. The API has been designed to be extensible,
//! but breaking changes will be necessary to implement several planned features. While it is
//! backed by an extensive test suite, it has yet to be thoroughly battle-tested in real
//! applications. Moreover, the focus so far has been on correctness and API design, with little
//! focus on performance. Future improvements will include incremental performance enhancements.
//!
//! Current limitations:
//!
//! - Limited or no availability of sparse system solvers.
//! - Limited support for complex numbers. Currently only arithmetic operations that do not
//!   rely on particular properties of complex numbers, such as e.g. conjugation, are
//!   supported.
//! - No integration with external libraries.
//!
//! # Usage
//!
//! Add the following to your `Cargo.toml` file:
//!
//! ```toml
//! [dependencies]
//! nalgebra_sparse = "0.1"
//! ```
//!
//! # Supported matrix formats
//!
//! | Format                  | Notes                                        |
//! | ------------------------|--------------------------------------------- |
//! | [COO](`coo::CooMatrix`) | Well-suited for matrix construction. <br /> Ill-suited for algebraic operations. |
//! | [CSR](`csr::CsrMatrix`) | Immutable sparsity pattern, suitable for algebraic operations. <br /> Fast row access. |
//! | [CSC](`csc::CscMatrix`) | Immutable sparsity pattern, suitable for algebraic operations. <br /> Fast column access. |
//!
//! What format is best to use depends on the application. The most common use case for sparse
//! matrices in science is the solution of sparse linear systems. Here we can differentiate between
//! two common cases:
//!
//! - Direct solvers. Typically, direct solvers take their input in CSR or CSC format.
//! - Iterative solvers. Many iterative solvers require only matrix-vector products,
//!   for which the CSR or CSC formats are suitable.
//!
//! The [COO](coo::CooMatrix) format is primarily intended for matrix construction.
//! A common pattern is to use COO for construction, before converting to CSR or CSC for use
//! in a direct solver or for computing matrix-vector products in an iterative solver.
//! Some high-performance applications might also directly manipulate the CSR and/or CSC
//! formats.
//!
//! # Example: COO -> CSR -> matrix-vector product
//!
//! ```
//! use nalgebra_sparse::{coo::CooMatrix, csr::CsrMatrix};
//! use nalgebra::{DMatrix, DVector};
//! use matrixcompare::assert_matrix_eq;
//!
//! // The dense representation of the matrix
//! let dense = DMatrix::from_row_slice(3, 3,
//!     &[1.0, 0.0, 3.0,
//!       2.0, 0.0, 1.3,
//!       0.0, 0.0, 4.1]);
//!
//! // Build the equivalent COO representation. We only add the non-zero values
//! let mut coo = CooMatrix::new(3, 3);
//! // We can add elements in any order. For clarity, we do so in row-major order here.
//! coo.push(0, 0, 1.0);
//! coo.push(0, 2, 3.0);
//! coo.push(1, 0, 2.0);
//! coo.push(1, 2, 1.3);
//! coo.push(2, 2, 4.1);
//!
//! // ... or add entire dense matrices like so:
//! // coo.push_matrix(0, 0, &dense);
//!
//! // The simplest way to construct a CSR matrix is to first construct a COO matrix, and
//! // then convert it to CSR. The `From` trait is implemented for conversions between different
//! // sparse matrix types.
//! // Alternatively, we can construct a matrix directly from the CSR data.
//! // See the docs for CsrMatrix for how to do that.
//! let csr = CsrMatrix::from(&coo);
//!
//! // Let's check that the CSR matrix and the dense matrix represent the same matrix.
//! // We can use macros from the `matrixcompare` crate to easily do this, despite the fact that
//! // we're comparing across two different matrix formats. Note that these macros are only really
//! // appropriate for writing tests, however.
//! assert_matrix_eq!(csr, dense);
//!
//! let x = DVector::from_column_slice(&[1.3, -4.0, 3.5]);
//!
//! // Compute the matrix-vector product y = A * x. We don't need to specify the type here,
//! // but let's just do it to make sure we get what we expect
//! let y: DVector<_> = &csr * &x;
//!
//! // Verify the result with a small element-wise absolute tolerance
//! let y_expected = DVector::from_column_slice(&[11.8, 7.15, 14.35]);
//! assert_matrix_eq!(y, y_expected, comp = abs, tol = 1e-9);
//!
//! // The above expression is simple, and gives easy to read code, but if we're doing this in a
//! // loop, we'll have to keep allocating new vectors. If we determine that this is a bottleneck,
//! // then we can resort to the lower level APIs for more control over the operations
//! {
//!     use nalgebra_sparse::ops::{Op, serial::spmm_csr_dense};
//!     let mut y = y;
//!     // Compute y <- 0.0 * y + 1.0 * csr * dense. We store the result directly in `y`, without
//!     // any intermediate allocations
//!     spmm_csr_dense(0.0, &mut y, 1.0, Op::NoOp(&csr), Op::NoOp(&x));
//!     assert_matrix_eq!(y, y_expected, comp = abs, tol = 1e-9);
//! }
//! ```
#![deny(non_camel_case_types)]
#![deny(unused_parens)]
#![deny(non_upper_case_globals)]
#![deny(unused_qualifications)]
#![deny(unused_results)]
#![deny(missing_docs)]

pub extern crate nalgebra as na;
pub mod convert;
pub mod coo;
pub mod csc;
pub mod csr;
pub mod factorization;
pub mod ops;
pub mod pattern;

pub(crate) mod cs;

#[cfg(feature = "proptest-support")]
pub mod proptest;

#[cfg(feature = "compare")]
mod matrixcompare;

use num_traits::Zero;
use std::error::Error;
use std::fmt;

pub use self::coo::CooMatrix;
pub use self::csc::CscMatrix;
pub use self::csr::CsrMatrix;

/// Errors produced by functions that expect well-formed sparse format data.
#[derive(Debug)]
pub struct SparseFormatError {
    kind: SparseFormatErrorKind,
    // Currently we only use an underlying error for generating the `Display` impl
    error: Box<dyn Error>,
}

impl SparseFormatError {
    /// The type of error.
    #[must_use]
    pub fn kind(&self) -> &SparseFormatErrorKind {
        &self.kind
    }

    pub(crate) fn from_kind_and_error(kind: SparseFormatErrorKind, error: Box<dyn Error>) -> Self {
        Self { kind, error }
    }

    /// Helper functionality for more conveniently creating errors.
    pub(crate) fn from_kind_and_msg(kind: SparseFormatErrorKind, msg: &'static str) -> Self {
        Self::from_kind_and_error(kind, Box::<dyn Error>::from(msg))
    }
}

/// The type of format error described by a [SparseFormatError](struct.SparseFormatError.html).
#[non_exhaustive]
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum SparseFormatErrorKind {
    /// Indicates that the index data associated with the format contains at least one index
    /// out of bounds.
    IndexOutOfBounds,

    /// Indicates that the provided data contains at least one duplicate entry, and the
    /// current format does not support duplicate entries.
    DuplicateEntry,

    /// Indicates that the provided data for the format does not conform to the high-level
    /// structure of the format.
    ///
    /// For example, the arrays defining the format data might have incompatible sizes.
    InvalidStructure,
}

impl fmt::Display for SparseFormatError {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "{}", self.error)
    }
}

impl Error for SparseFormatError {}

/// An entry in a sparse matrix.
///
/// Sparse matrices do not store all their entries explicitly. Therefore, entry (i, j) in the matrix
/// can either be a reference to an explicitly stored element, or it is implicitly zero.
#[derive(Debug, PartialEq, Eq)]
pub enum SparseEntry<'a, T> {
    /// The entry is a reference to an explicitly stored element.
    ///
    /// Note that the naming here is a misnomer: The element can still be zero, even though it
    /// is explicitly stored (a so-called "explicit zero").
    NonZero(&'a T),
    /// The entry is implicitly zero, i.e. it is not explicitly stored.
    Zero,
}

impl<'a, T: Clone + Zero> SparseEntry<'a, T> {
    /// Returns the value represented by this entry.
    ///
    /// Either clones the underlying reference or returns zero if the entry is not explicitly
    /// stored.
    pub fn into_value(self) -> T {
        match self {
            SparseEntry::NonZero(value) => value.clone(),
            SparseEntry::Zero => T::zero(),
        }
    }
}

/// A mutable entry in a sparse matrix.
///
/// See also `SparseEntry`.
#[derive(Debug, PartialEq, Eq)]
pub enum SparseEntryMut<'a, T> {
    /// The entry is a mutable reference to an explicitly stored element.
    ///
    /// Note that the naming here is a misnomer: The element can still be zero, even though it
    /// is explicitly stored (a so-called "explicit zero").
    NonZero(&'a mut T),
    /// The entry is implicitly zero i.e. it is not explicitly stored.
    Zero,
}

impl<'a, T: Clone + Zero> SparseEntryMut<'a, T> {
    /// Returns the value represented by this entry.
    ///
    /// Either clones the underlying reference or returns zero if the entry is not explicitly
    /// stored.
    pub fn into_value(self) -> T {
        match self {
            SparseEntryMut::NonZero(value) => value.clone(),
            SparseEntryMut::Zero => T::zero(),
        }
    }
}