diff --git a/nalgebra-sparse/Cargo.toml b/nalgebra-sparse/Cargo.toml index 3f9f8c8a..27d90a90 100644 --- a/nalgebra-sparse/Cargo.toml +++ b/nalgebra-sparse/Cargo.toml @@ -21,3 +21,7 @@ matrixcompare-core = { version = "0.1.0", optional = true } itertools = "0.9" matrixcompare = { version = "0.2.0", features = [ "proptest-support" ] } nalgebra = { version="0.23", path = "../", features = ["compare"] } + +[package.metadata.docs.rs] +# Enable certain features when building docs for docs.rs +features = [ "proptest-support", "compare" ] \ No newline at end of file diff --git a/nalgebra-sparse/src/convert/mod.rs b/nalgebra-sparse/src/convert/mod.rs index cf627dfa..77388b22 100644 --- a/nalgebra-sparse/src/convert/mod.rs +++ b/nalgebra-sparse/src/convert/mod.rs @@ -1,4 +1,39 @@ -//! TODO +//! Routines for converting between sparse matrix formats. +//! +//! Most users should instead use the provided `From` implementations to convert between matrix +//! formats. Note that `From` implementations may not be available between all combinations of +//! sparse matrices. +//! +//! The following example illustrates how to convert between matrix formats with the `From` +//! implementations. +//! +//! ```rust +//! use nalgebra_sparse::{csr::CsrMatrix, csc::CscMatrix, coo::CooMatrix}; +//! use nalgebra::DMatrix; +//! +//! // Conversion from dense +//! let dense = DMatrix::::identity(9, 8); +//! let csr = CsrMatrix::from(&dense); +//! let csc = CscMatrix::from(&dense); +//! let coo = CooMatrix::from(&dense); +//! +//! // CSR <-> CSC +//! let _ = CsrMatrix::from(&csc); +//! let _ = CscMatrix::from(&csr); +//! +//! // CSR <-> COO +//! let _ = CooMatrix::from(&csr); +//! let _ = CsrMatrix::from(&coo); +//! +//! // CSC <-> COO +//! let _ = CooMatrix::from(&csc); +//! let _ = CscMatrix::from(&coo); +//! ``` +//! +//! The routines available here are able to provide more specialized APIs, giving +//! more control over the conversion process. The routines are organized by backends. +//! Currently, only the [`serial`] backend is available. +//! In the future, backends that offer parallel routines may become available. pub mod serial; diff --git a/nalgebra-sparse/src/convert/serial.rs b/nalgebra-sparse/src/convert/serial.rs index 4ddefcca..67a57966 100644 --- a/nalgebra-sparse/src/convert/serial.rs +++ b/nalgebra-sparse/src/convert/serial.rs @@ -1,4 +1,8 @@ -//! TODO +//! Serial routines for converting between matrix formats. +//! +//! All routines in this module are single-threaded. At present these routines offer no +//! advantage over using the [`From`] trait, but future changes to the API might offer more +//! control to the user. use std::ops::Add; use num_traits::Zero; @@ -11,7 +15,7 @@ use crate::cs; use crate::csc::CscMatrix; use crate::csr::CsrMatrix; -/// TODO +/// Converts a dense matrix to [`CooMatrix`]. pub fn convert_dense_coo(dense: &Matrix) -> CooMatrix where T: Scalar + Zero, @@ -33,9 +37,7 @@ where coo } -/// TODO -/// -/// TODO: What should the actual trait bounds be? +/// Converts a [`CooMatrix`] to a dense matrix. pub fn convert_coo_dense(coo: &CooMatrix) -> DMatrix where T: Scalar + Zero + ClosedAdd, @@ -47,7 +49,7 @@ where output } -/// TODO +/// Converts a [`CooMatrix`] to a [`CsrMatrix`]. pub fn convert_coo_csr(coo: &CooMatrix) -> CsrMatrix where T: Scalar + Zero @@ -63,7 +65,7 @@ where .expect("Internal error: Invalid CSR data during COO->CSR conversion") } -/// TODO +/// Converts a [`CsrMatrix`] to a [`CooMatrix`]. pub fn convert_csr_coo(csr: &CsrMatrix) -> CooMatrix { let mut result = CooMatrix::new(csr.nrows(), csr.ncols()); @@ -73,7 +75,7 @@ pub fn convert_csr_coo(csr: &CsrMatrix) -> CooMatrix result } -/// TODO +/// Converts a [`CsrMatrix`] to a dense matrix. pub fn convert_csr_dense(csr:& CsrMatrix) -> DMatrix where T: Scalar + ClosedAdd + Zero @@ -87,7 +89,7 @@ where output } -/// TODO +/// Converts a dense matrix to a [`CsrMatrix`]. pub fn convert_dense_csr(dense: &Matrix) -> CsrMatrix where T: Scalar + Zero, @@ -120,7 +122,7 @@ where .expect("Internal error: Invalid CsrMatrix format during dense-> CSR conversion") } -/// TODO +/// Converts a [`CooMatrix`] to a [`CscMatrix`]. pub fn convert_coo_csc(coo: &CooMatrix) -> CscMatrix where T: Scalar + Zero @@ -136,7 +138,7 @@ where .expect("Internal error: Invalid CSC data during COO->CSC conversion") } -/// TODO +/// Converts a [`CscMatrix`] to a [`CooMatrix`]. pub fn convert_csc_coo(csc: &CscMatrix) -> CooMatrix where T: Scalar @@ -148,7 +150,7 @@ where coo } -/// TODO +/// Converts a [`CscMatrix`] to a dense matrix. pub fn convert_csc_dense(csc: &CscMatrix) -> DMatrix where T: Scalar + ClosedAdd + Zero @@ -162,7 +164,7 @@ where output } -/// TODO +/// Converts a dense matrix to a [`CscMatrix`]. pub fn convert_dense_csc(dense: &Matrix) -> CscMatrix where T: Scalar + Zero, @@ -192,7 +194,7 @@ pub fn convert_dense_csc(dense: &Matrix) -> CscMatrix .expect("Internal error: Invalid CscMatrix format during dense-> CSC conversion") } -/// TODO +/// Converts a [`CsrMatrix`] to a [`CscMatrix`]. pub fn convert_csr_csc(csr: &CsrMatrix) -> CscMatrix where T: Scalar @@ -208,7 +210,7 @@ where .expect("Internal error: Invalid CSC data during CSR->CSC conversion") } -/// TODO +/// Converts a [`CscMatrix`] to a [`CsrMatrix`]. pub fn convert_csc_csr(csc: &CscMatrix) -> CsrMatrix where T: Scalar diff --git a/nalgebra-sparse/src/factorization/cholesky.rs b/nalgebra-sparse/src/factorization/cholesky.rs index 8ffbbe3d..1c4f95a8 100644 --- a/nalgebra-sparse/src/factorization/cholesky.rs +++ b/nalgebra-sparse/src/factorization/cholesky.rs @@ -1,6 +1,3 @@ -// TODO: Remove this allowance -#![allow(missing_docs)] - use crate::pattern::SparsityPattern; use crate::csc::CscMatrix; use core::{mem, iter}; @@ -9,6 +6,10 @@ use std::fmt::{Display, Formatter}; use crate::ops::serial::spsolve_csc_lower_triangular; use crate::ops::Op; +/// A symbolic sparse Cholesky factorization of a CSC matrix. +/// +/// The symbolic factorization computes the sparsity pattern of `L`, the Cholesky factor. +#[derive(Debug, Clone, PartialEq, Eq)] pub struct CscSymbolicCholesky { // Pattern of the original matrix that was decomposed m_pattern: SparsityPattern, @@ -18,6 +19,14 @@ pub struct CscSymbolicCholesky { } impl CscSymbolicCholesky { + /// Compute the symbolic factorization for a sparsity pattern belonging to a CSC matrix. + /// + /// The sparsity pattern must be symmetric. However, this is not enforced, and it is the + /// responsibility of the user to ensure that this property holds. + /// + /// # Panics + /// + /// Panics if the sparsity pattern is not square. pub fn factor(pattern: SparsityPattern) -> Self { assert_eq!(pattern.major_dim(), pattern.minor_dim(), "Major and minor dimensions must be the same (square matrix)."); @@ -29,11 +38,27 @@ impl CscSymbolicCholesky { } } + /// The pattern of the Cholesky factor `L`. pub fn l_pattern(&self) -> &SparsityPattern { &self.l_pattern } } +/// A sparse Cholesky factorization `A = L L^T` of a [`CscMatrix`]. +/// +/// The factor `L` is a sparse, lower-triangular matrix. See the article on [Wikipedia] for +/// more information. +/// +/// The implementation is a port of the `CsCholesky` implementation in `nalgebra`. It is similar +/// to Tim Davis' [`CSparse`]. The current implementation performs no fill-in reduction, and can +/// therefore be expected to produce much too dense Cholesky factors for many matrices. +/// It is therefore not currently recommended to use this implementation for serious projects. +/// +/// [`CSparse`]: https://epubs.siam.org/doi/book/10.1137/1.9780898718881 +/// [Wikipedia]: https://en.wikipedia.org/wiki/Cholesky_decomposition +// TODO: We should probably implement PartialEq/Eq, but in that case we'd probably need a +// custom implementation, due to the need to exclude the workspace arrays +#[derive(Debug, Clone)] pub struct CscCholesky { // Pattern of the original matrix m_pattern: SparsityPattern, @@ -45,6 +70,7 @@ pub struct CscCholesky { #[derive(Debug, PartialEq, Eq, Clone)] #[non_exhaustive] +/// Possible errors produced by the Cholesky factorization. pub enum CholeskyError { /// The matrix is not positive definite. NotPositiveDefinite, @@ -59,7 +85,24 @@ impl Display for CholeskyError { impl std::error::Error for CholeskyError {} impl CscCholesky { - pub fn factor_numerical(symbolic: CscSymbolicCholesky, values: &[T]) -> Result { + /// Computes the numerical Cholesky factorization associated with the given + /// symbolic factorization and the provided values. + /// + /// The values correspond to the non-zero values of the CSC matrix for which the + /// symbolic factorization was computed. + /// + /// # Errors + /// + /// Returns an error if the numerical factorization fails. This can occur if the matrix is not + /// symmetric positive definite. + /// + /// # Panics + /// + /// Panics if the number of values differ from the number of non-zeros of the sparsity pattern + /// of the matrix that was symbolically factored. + pub fn factor_numerical(symbolic: CscSymbolicCholesky, values: &[T]) + -> Result + { assert_eq!(symbolic.l_pattern.nnz(), symbolic.u_pattern.nnz(), "u is just the transpose of l, so should have the same nnz"); @@ -84,19 +127,50 @@ impl CscCholesky { Ok(factorization) } + /// Computes the Cholesky factorization of the provided matrix. + /// + /// The matrix must be symmetric positive definite. Symmetry is not checked, and it is up + /// to the user to enforce this property. + /// + /// # Errors + /// + /// Returns an error if the numerical factorization fails. This can occur if the matrix is not + /// symmetric positive definite. + /// + /// # Panics + /// + /// Panics if the matrix is not square. + /// + /// TODO: Take matrix by value or not? pub fn factor(matrix: &CscMatrix) -> Result { let symbolic = CscSymbolicCholesky::factor(matrix.pattern().clone()); Self::factor_numerical(symbolic, matrix.values()) } + /// Re-computes the factorization for a new set of non-zero values. + /// + /// This is useful when the values of a matrix changes, but the sparsity pattern remains + /// constant. + /// + /// # Errors + /// + /// Returns an error if the numerical factorization fails. This can occur if the matrix is not + /// symmetric positive definite. + /// + /// # Panics + /// + /// Panics if the number of values does not match the number of non-zeros in the sparsity + /// pattern. pub fn refactor(&mut self, values: &[T]) -> Result<(), CholeskyError> { self.decompose_left_looking(values) } + /// Returns a reference to the Cholesky factor `L`. pub fn l(&self) -> &CscMatrix { &self.l_factor } + /// Returns the Cholesky factor `L`. pub fn take_l(self) -> CscMatrix { self.l_factor } @@ -178,6 +252,11 @@ impl CscCholesky { Ok(()) } + /// Solves the system `A X = B`, where `X` and `B` are dense matrices. + /// + /// # Panics + /// + /// Panics if `B` is not square. pub fn solve<'a>(&'a self, b: impl Into>) -> DMatrix { let b = b.into(); let mut output = b.clone_owned(); @@ -185,6 +264,13 @@ impl CscCholesky { output } + /// Solves the system `AX = B`, where `X` and `B` are dense matrices. + /// + /// The result is stored in-place in `b`. + /// + /// # Panics + /// + /// Panics if `b` is not square. pub fn solve_mut<'a>(&'a self, b: impl Into>) { let expect_msg = "If the Cholesky factorization succeeded,\ @@ -201,9 +287,6 @@ impl CscCholesky { } } - - - fn reach( pattern: &SparsityPattern, j: usize, diff --git a/nalgebra-sparse/src/factorization/mod.rs b/nalgebra-sparse/src/factorization/mod.rs index 6596f530..b4adcba2 100644 --- a/nalgebra-sparse/src/factorization/mod.rs +++ b/nalgebra-sparse/src/factorization/mod.rs @@ -1,4 +1,6 @@ //! Matrix factorization for sparse matrices. +//! +//! Currently, the only factorization provided here is the [`CscCholesky`] factorization. mod cholesky; pub use cholesky::*; \ No newline at end of file diff --git a/nalgebra-sparse/src/lib.rs b/nalgebra-sparse/src/lib.rs index fe4255a8..73c0d88b 100644 --- a/nalgebra-sparse/src/lib.rs +++ b/nalgebra-sparse/src/lib.rs @@ -158,17 +158,26 @@ impl fmt::Display for SparseFormatError { impl Error for SparseFormatError {} -/// TODO +/// 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> { - /// TODO + /// 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), - /// TODO + /// The entry is implicitly zero, i.e. it is not explicitly stored. Zero } impl<'a, T: Clone + Zero> SparseEntry<'a, T> { - /// TODO + /// Returns the value represented by this entry. + /// + /// Either clones the underlying reference or returns zero if the entry is not explicitly + /// stored. pub fn to_value(self) -> T { match self { SparseEntry::NonZero(value) => value.clone(), @@ -177,17 +186,25 @@ impl<'a, T: Clone + Zero> SparseEntry<'a, T> { } } -/// TODO +/// A mutable entry in a sparse matrix. +/// +/// See also `SparseEntry`. #[derive(Debug, PartialEq, Eq)] pub enum SparseEntryMut<'a, T> { - /// TODO + /// 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), - /// TODO + /// The entry is implicitly zero i.e. it is not explicitly stored. Zero } impl<'a, T: Clone + Zero> SparseEntryMut<'a, T> { - /// TODO + /// Returns the value represented by this entry. + /// + /// Either clones the underlying reference or returns zero if the entry is not explicitly + /// stored. pub fn to_value(self) -> T { match self { SparseEntryMut::NonZero(value) => value.clone(), diff --git a/nalgebra-sparse/src/ops/mod.rs b/nalgebra-sparse/src/ops/mod.rs index dffcfa22..2857533f 100644 --- a/nalgebra-sparse/src/ops/mod.rs +++ b/nalgebra-sparse/src/ops/mod.rs @@ -1,24 +1,37 @@ -//! TODO +//! Sparse matrix arithmetic operations. +//! +//! TODO: Explain that users should prefer to use std ops unless they need to get more performance +//! +//! The available operations are organized by backend. Currently, only the [`serial`] backend +//! is available. In the future, backends that expose parallel operations may become available. +//! +//! Many routines are able to implicitly transpose matrices involved in the operation. +//! For example, the routine [`spadd_csr_prealloc`](serial::spadd_csr_prealloc) performs the +//! operation `C <- beta * C + alpha * op(A)`. Here `op(A)` indicates that the matrix `A` can +//! either be used as-is or transposed. The notation `op(A)` is represented in code by the +//! [`Op`] enum. mod impl_std_ops; pub mod serial; -/// TODO +/// Determines whether a matrix should be transposed in a given operation. +/// +/// See the [module-level documentation](crate::ops) for the purpose of this enum. #[derive(Debug, Copy, Clone, PartialEq, Eq)] pub enum Op { - /// TODO + /// Indicates that the matrix should be used as-is. NoOp(T), - /// TODO + /// Indicates that the matrix should be transposed. Transpose(T), } impl Op { - /// TODO + /// Returns a reference to the inner value that the operation applies to. pub fn inner_ref(&self) -> &T { self.as_ref().into_inner() } - /// TODO + /// Returns an `Op` applied to a reference of the inner value. pub fn as_ref(&self) -> Op<&T> { match self { Op::NoOp(obj) => Op::NoOp(&obj), @@ -26,15 +39,14 @@ impl Op { } } - /// TODO + /// Converts the underlying data type. pub fn convert(self) -> Op where T: Into { self.map_same_op(T::into) } - /// TODO - /// TODO: Rewrite the other functions by leveraging this one + /// Transforms the inner value with the provided function, but preserves the operation. pub fn map_same_op U>(self, f: F) -> Op { match self { Op::NoOp(obj) => Op::NoOp(f(obj)), @@ -42,7 +54,7 @@ impl Op { } } - /// TODO + /// Consumes the `Op` and returns the inner value. pub fn into_inner(self) -> T { match self { Op::NoOp(obj) | Op::Transpose(obj) => obj, diff --git a/nalgebra-sparse/src/ops/serial/cs.rs b/nalgebra-sparse/src/ops/serial/cs.rs index 9a4f7e34..dd73c304 100644 --- a/nalgebra-sparse/src/ops/serial/cs.rs +++ b/nalgebra-sparse/src/ops/serial/cs.rs @@ -1,13 +1,13 @@ use crate::cs::CsMatrix; use crate::ops::Op; -use crate::ops::serial::{OperationErrorType, OperationError}; +use crate::ops::serial::{OperationErrorKind, OperationError}; use nalgebra::{Scalar, ClosedAdd, ClosedMul, DMatrixSliceMut, DMatrixSlice}; use num_traits::{Zero, One}; use crate::SparseEntryMut; fn spmm_cs_unexpected_entry() -> OperationError { - OperationError::from_type_and_message( - OperationErrorType::InvalidPattern, + OperationError::from_kind_and_message( + OperationErrorKind::InvalidPattern, String::from("Found unexpected entry that is not present in `c`.")) } @@ -58,8 +58,8 @@ pub fn spmm_cs_prealloc( } fn spadd_cs_unexpected_entry() -> OperationError { - OperationError::from_type_and_message( - OperationErrorType::InvalidPattern, + OperationError::from_kind_and_message( + OperationErrorKind::InvalidPattern, String::from("Found entry in `op(a)` that is not present in `c`.")) } diff --git a/nalgebra-sparse/src/ops/serial/mod.rs b/nalgebra-sparse/src/ops/serial/mod.rs index ea6e6d2b..38ee266a 100644 --- a/nalgebra-sparse/src/ops/serial/mod.rs +++ b/nalgebra-sparse/src/ops/serial/mod.rs @@ -1,4 +1,12 @@ -//! TODO +//! Serial sparse matrix arithmetic routines. +//! +//! All routines are single-threaded. +//! +//! Some operations have the `prealloc` suffix. This means that they expect that the sparsity +//! pattern of the output matrix has already been pre-allocated: that is, the pattern of the result +//! of the operation fits entirely in the output pattern. In the future, there will also be +//! some operations which will be able to dynamically adapt the output pattern to fit the +//! result, but these have yet to be implemented. #[macro_use] macro_rules! assert_compatible_spmm_dims { @@ -58,24 +66,32 @@ pub use csc::*; pub use csr::*; pub use pattern::*; -/// TODO +/// A description of the error that occurred during an arithmetic operation. #[derive(Clone, Debug)] pub struct OperationError { - error_type: OperationErrorType, + error_kind: OperationErrorKind, message: String } -/// TODO +/// The different kinds of operation errors that may occur. #[non_exhaustive] #[derive(Clone, Debug)] -pub enum OperationErrorType { - /// TODO +pub enum OperationErrorKind { + /// Indicates that one or more sparsity patterns involved in the operation violate the + /// expectations of the routine. + /// + /// For example, this could indicate that the sparsity pattern of the output is not able to + /// contain the result of the operation. InvalidPattern, } impl OperationError { - /// TODO - pub fn from_type_and_message(error_type: OperationErrorType, message: String) -> Self { - Self { error_type, message } + fn from_kind_and_message(error_type: OperationErrorKind, message: String) -> Self { + Self { error_kind: error_type, message } + } + + /// The operation error kind. + pub fn kind(&self) -> &OperationErrorKind { + &self.error_kind } } \ No newline at end of file diff --git a/nalgebra-sparse/src/proptest.rs b/nalgebra-sparse/src/proptest.rs index 143696ad..43b80896 100644 --- a/nalgebra-sparse/src/proptest.rs +++ b/nalgebra-sparse/src/proptest.rs @@ -1,6 +1,10 @@ -//! TODO +//! Functionality for integrating `nalgebra-sparse` with `proptest`. //! -//! TODO: Clarify that this module needs proptest-support feature +//! **This module is only available if the `proptest-support` feature is enabled**. +//! +//! The strategies provided here are generally expected to be able to generate the entire range +//! of possible outputs given the constraints on dimensions and values. However, there are no +//! particular guarantees on the distribution of possible values. // Contains some patched code from proptest that we can remove in the (hopefully near) future. // See docs in file for more details. @@ -139,7 +143,12 @@ fn sparse_triplet_strategy(value_strategy: T, .prop_shuffle() } -/// TODO +/// A strategy for producing COO matrices without duplicate entries. +/// +/// The values of the matrix are picked from the provided `value_strategy`, while the size of the +/// generated matrices is determined by the ranges `rows` and `cols`. The number of explicitly +/// stored entries is bounded from above by `max_nonzeros`. Note that the matrix might still +/// contain explicitly stored zeroes if the value strategy is capable of generating zero values. pub fn coo_no_duplicates( value_strategy: T, rows: impl Into, @@ -177,10 +186,17 @@ where }) } -/// TODO +/// A strategy for producing COO matrices with duplicate entries. /// -/// TODO: Write note on how this strategy only maintains the constraints on values -/// for each triplet, but does not consider the sum of triplets +/// The values of the matrix are picked from the provided `value_strategy`, while the size of the +/// generated matrices is determined by the ranges `rows` and `cols`. Note that the values +/// only apply to individual entries, and since this strategy can generate duplicate entries, +/// the matrix will generally have values outside the range determined by `value_strategy` when +/// converted to other formats, since the duplicate entries are summed together in this case. +/// +/// The number of explicitly stored entries is bounded from above by `max_nonzeros`. The maximum +/// number of duplicate entries is determined by `max_duplicates`. Note that the matrix might still +/// contain explicitly stored zeroes if the value strategy is capable of generating zero values. pub fn coo_with_duplicates( value_strategy: T, rows: impl Into, @@ -255,7 +271,7 @@ where .expect("Internal error: Generated sparsity pattern is invalid") } -/// TODO +/// A strategy for generating sparsity patterns. pub fn sparsity_pattern( major_lanes: impl Into, minor_lanes: impl Into, @@ -286,7 +302,7 @@ pub fn sparsity_pattern( }) } -/// TODO +/// A strategy for generating CSR matrices. pub fn csr(value_strategy: T, rows: impl Into, cols: impl Into, @@ -310,7 +326,7 @@ where }) } -/// TODO +/// A strategy for generating CSC matrices. pub fn csc(value_strategy: T, rows: impl Into, cols: impl Into,