Skip to content

Commit

Permalink
doc updates
Browse files Browse the repository at this point in the history
  • Loading branch information
goulart-paul committed Feb 2, 2025
1 parent 08430d2 commit d7ca8b4
Show file tree
Hide file tree
Showing 23 changed files with 203 additions and 63 deletions.
4 changes: 2 additions & 2 deletions src/algebra/csc/core.rs
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ pub struct CscMatrix<T = f64> {
pub n: usize,
/// CSC format column pointer.
///
/// Ths field should have length `n+1`. The last entry corresponds
/// This field should have length `n+1`. The last entry corresponds
/// to the the number of nonzeros and should agree with the lengths
/// of the `rowval` and `nzval` fields.
pub colptr: Vec<usize>,
Expand Down Expand Up @@ -419,7 +419,7 @@ where
self.size() == other.size() && self.colptr == other.colptr && self.rowval == other.rowval
}

/// Same as is_equal_sparsity, but returns an error indicating the reason
/// Same as `is_equal_sparsity`, but returns an error indicating the reason
/// for failure if the matrices do not have equivalent sparsity patterns.
pub fn check_equal_sparsity(&self, other: &Self) -> Result<(), SparseFormatError> {
if self.size() != other.size() {
Expand Down
2 changes: 1 addition & 1 deletion src/algebra/dense/types.rs
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ where
{
/// dimensions
pub size: (usize, usize),
/// vector of data in column major formmat
/// vector of data in column major format
pub data: S,
pub(crate) phantom: std::marker::PhantomData<T>,
}
Expand Down
6 changes: 6 additions & 0 deletions src/algebra/error_types.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4,21 +4,27 @@ use thiserror::Error;
#[derive(Error, Debug)]
pub enum MatrixConcatenationError {
#[error("Incompatible dimensions")]
/// Indicates inputs have incompatible dimension
IncompatibleDimension,
}

#[derive(Error, Debug)]
/// Error type returned by sparse matrix assembly operations.
pub enum SparseFormatError {
/// Matrix dimension fields and/or array lengths are incompatible
#[error("Matrix dimension fields and/or array lengths are incompatible")]
IncompatibleDimension,
/// Data is not sorted by row index within each column
#[error("Data is not sorted by row index within each column")]
BadRowOrdering,
#[error("Row value exceeds the matrix row dimension")]
/// Row value exceeds the matrix row dimension
BadRowval,
#[error("Bad column pointer values")]
/// Matrix column pointer values are defective
BadColptr,
#[error("sparsity pattern mismatch")]
/// Operation on matrices that have mismatching sparsity patterns
SparsityMismatch,
}

Expand Down
18 changes: 14 additions & 4 deletions src/algebra/floats.rs
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,10 @@ use crate::algebra::dense::BlasFloatT;
/// Core traits for internal floating point values.
///
/// This trait defines a subset of bounds for `FloatT`, which is preferred
/// throughout for use in the solver. When the "sdp" feature is enabled,
/// throughout for use in the solver. When the `sdp` feature is enabled,
/// `FloatT` is additionally restricted to f32/f64 types supported by BLAS.
/// When the `faer-sparse` feature is enabled, `FloatT` is additionally
/// restricted to types implementing `RealField` from the `faer` crate.
pub trait CoreFloatT:
'static
+ Send
Expand Down Expand Up @@ -51,24 +53,32 @@ impl<T> CoreFloatT for T where

cfg_if::cfg_if! {
if #[cfg(feature="sdp")] {
/// A marker trait implemented on types supported by BLAS (i.e. f32 and f64)
/// when the package is compiled with the "sdp" feature using a BLAS/LAPACK library
#[doc(hidden)]
pub trait MaybeBlasFloatT : BlasFloatT {}
impl<T> MaybeBlasFloatT for T where T: BlasFloatT {}
}
else {
#[doc(hidden)]
pub trait MaybeBlasFloatT {}
impl<T> MaybeBlasFloatT for T {}
}
}

// if "faer" is enabled, we must add an additional bound
// to restrict compilation to Reals implementing SimpleEntity
// to restrict compilation to types implementing RealField

cfg_if::cfg_if! {
if #[cfg(feature="faer-sparse")] {
#[doc(hidden)]
/// A marker trait implemented on types supported by faer-rs
/// when the package is compiled with the "faer-sparse" feature
pub trait MaybeFaerFloatT : faer_traits::RealField {}
impl<T> MaybeFaerFloatT for T where T: faer_traits::RealField {}
}
else {
#[doc(hidden)]
pub trait MaybeFaerFloatT {}
impl<T> MaybeFaerFloatT for T {}
}
Expand All @@ -86,7 +96,7 @@ cfg_if::cfg_if! {
pub trait FloatT: CoreFloatT + MaybeBlasFloatT + MaybeFaerFloatT {}
impl<T> FloatT for T where T: CoreFloatT + MaybeBlasFloatT + MaybeFaerFloatT {}

/// Trait for convering Rust primitives to [`FloatT`](crate::algebra::FloatT)
/// Trait for converting Rust primitives to [`FloatT`](crate::algebra::FloatT)
///
/// This convenience trait is implemented on f32/64 and u32/64. This trait
/// is required internally by the solver for converting constant primitives
Expand All @@ -97,7 +107,7 @@ impl<T> FloatT for T where T: CoreFloatT + MaybeBlasFloatT + MaybeFaerFloatT {}
// NB: `AsFloatT` is a convenience trait for f32/64 and u32/64
// so that we can do things like (2.0).as_T() everywhere on
// constants, rather than the awful T::from_f32(2.0).unwrap()
pub trait AsFloatT<T>: 'static {
pub(crate) trait AsFloatT<T>: 'static {
fn as_T(&self) -> T;
}

Expand Down
6 changes: 3 additions & 3 deletions src/algebra/math_traits.rs
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ pub trait VectorMath<T> {
/// Dot product
fn dot(&self, y: &Self) -> T;

// computes dot(z + αdz,s + αds) without intermediate allocation
/// computes dot(z + αdz,s + αds) without intermediate allocation
fn dot_shifted(z: &[T], s: &[T], dz: &[T], ds: &[T], α: T) -> T;

/// Standard Euclidian or 2-norm distance from `self` to `y`
Expand Down Expand Up @@ -93,7 +93,7 @@ pub trait VectorMath<T> {
/// Inf-norm of an elementwise scaling of `self` by `v`
fn norm_one_scaled(&self, v: &Self) -> T;

// Inf-norm of vector difference
/// Inf-norm of vector difference
fn norm_inf_diff(&self, b: &Self) -> T;

/// Minimum value in vector
Expand Down Expand Up @@ -152,7 +152,7 @@ pub trait MatrixMath<T> {
fn col_norms_no_reset(&self, norms: &mut [T]);

/// Compute columnwise infinity norm operations on
/// a symmstric matrix
/// a symmetric matrix
fn col_norms_sym(&self, norms: &mut [T]);

/// Compute columnwise infinity norm operations on
Expand Down
2 changes: 1 addition & 1 deletion src/algebra/matrix_traits.rs
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ pub trait BlockConcatenate: Sized {
/// Errors if column dimensions are incompatible
fn vcat(A: &Self, B: &Self) -> Result<Self, MatrixConcatenationError>;

/// general block concatentation
/// general block concatenation
fn hvcat(mats: &[&[&Self]]) -> Result<Self, MatrixConcatenationError>;

/// block diagonal concatenation
Expand Down
2 changes: 1 addition & 1 deletion src/algebra/sparsevector/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ pub(crate) struct SparseVector<T = f64> {
pub nzval: Vec<T>,
}

/// Creates a SparseVector from a dense slice.
/// Creates a `SparseVector` from a dense slice.
impl<T> SparseVector<T>
where
T: Num + Copy,
Expand Down
2 changes: 1 addition & 1 deletion src/python/mod.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
//! Clarabel Python interface.
//!
//! This module implements a wrapper for the Rust version of Python using
//! [PYO3](https://pyo3.rs/). To build these wrappers from `cargo`, compile the crate with
//! [PyO3](https://pyo3.rs/). To build these wrappers from `cargo`, compile the crate with
//! `--features python`. This module has no public API.
//!
//! It should not normally be necessary to compile the Python wrapper from
Expand Down
91 changes: 65 additions & 26 deletions src/qdldl/qdldl.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,39 +6,56 @@ use std::iter::zip;
use thiserror::Error;

/// Error codes returnable from [`QDLDLFactorisation`](QDLDLFactorisation) factor operations
#[derive(Error, Debug)]
pub enum QDLDLError {
#[error("Matrix dimension fields are incompatible")]
/// Matrix dimension fields are incompatible
IncompatibleDimension,
#[error("Matrix has a zero column")]
/// Matrix has a zero column
EmptyColumn,
#[error("Matrix is not upper triangular")]
/// Matrix is not upper triangular
NotUpperTriangular,
/// Matrix factorization produced a zero pivot
#[error("Matrix factorization produced a zero pivot")]
ZeroPivot,
#[error("Invalid permutation vector")]
/// Invalid permutation vector supplied
InvalidPermutation,
}

/// Required settings for [`QDLDLFactorisation`](QDLDLFactorisation)
#[derive(Builder, Debug, Clone)]
#[allow(missing_docs)]
/// Required settings for [`QDLDLFactorisation`](QDLDLFactorisation)
pub struct QDLDLSettings<T: FloatT> {
/// "dense scale" parameter for AMD ordering
#[builder(default = "1.0")]
amd_dense_scale: f64,
pub amd_dense_scale: f64,

/// optional ueser-supplied custom permutation vector for the matrix
#[builder(default = "None", setter(strip_option))]
perm: Option<Vec<usize>>,
pub perm: Option<Vec<usize>>,

/// Logical factorisation only, no numerical factorisation
#[builder(default = "false")]
logical: bool,
pub logical: bool,

/// optional user-supplied signs of the diagonal elements of D in LDL^T
#[builder(default = "None", setter(strip_option))]
Dsigns: Option<Vec<i8>>,
pub Dsigns: Option<Vec<i8>>,

/// Enable regularization during factorisation
#[builder(default = "true")]
regularize_enable: bool,
pub regularize_enable: bool,

/// Regularization epsilon parameter
#[builder(default = "(1e-12).as_T()")]
regularize_eps: T,
pub regularize_eps: T,

/// Regularization delta parameter
#[builder(default = "(1e-7).as_T()")]
regularize_delta: T,
pub regularize_delta: T,
}

impl<T> Default for QDLDLSettings<T>
Expand All @@ -51,29 +68,30 @@ where
}

/// Performs $LDL^T$ factorization of a symmetric quasidefinite matrix
#[derive(Debug)]
pub struct QDLDLFactorisation<T = f64> {
// permutation vector
/// permutation vector
pub perm: Vec<usize>,
// inverse permutation
/// inverse permutation
#[allow(dead_code)] //Unused because we call ipermute in solve instead. Keep anyway.
iperm: Vec<usize>,
// lower triangular factor
/// lower triangular factor L in LDL^T
pub L: CscMatrix<T>,
// D and is inverse for A = LDL^T
/// vector of diagonal elements of D in LDL^T
pub D: Vec<T>,
/// vector of reciprocal diagonal elements of D in LDL^T
pub Dinv: Vec<T>,
// workspace data
/// internal workspace data
workspace: QDLDLWorkspace<T>,
// is it logical factorisation only?
is_logical: bool,
/// true of factorisation is symbolic only
is_symbolic: bool,
}

impl<T> QDLDLFactorisation<T>
where
T: FloatT,
{
/// create a new LDL^T factorisation
pub fn new(
Ain: &CscMatrix<T>,
opts: Option<QDLDLSettings<T>>,
Expand All @@ -83,18 +101,21 @@ where
_qdldl_new(Ain, opts)
}

/// returns the number of positive eigenvalues
pub fn positive_inertia(&self) -> usize {
self.workspace.positive_inertia
}
/// returns the number of regularisation shifts
/// that were applied during factorisation
pub fn regularize_count(&self) -> usize {
self.workspace.regularize_count
}

// Solves Ax = b using LDL factors for A.
// Solves in place (x replaces b)
/// Solves Ax = b using LDL factors for A.
/// Solves in place (x replaces b)
pub fn solve(&mut self, b: &mut [T]) {
// bomb if logical factorisation only
assert!(!self.is_logical);
assert!(!self.is_symbolic);

// bomb if b is the wrong size
assert_eq!(b.len(), self.D.len());
Expand All @@ -116,6 +137,8 @@ where
ipermute(b, tmp, &self.perm);
}

/// Update a subset of the values of the matrix to be (re)factored. See [`refactor`](crate::qdldl::QDLDLFactorisation::refactor)
///
pub fn update_values(&mut self, indices: &[usize], values: &[T]) {
let nzval = &mut self.workspace.triuA.nzval; // post perm internal data
let AtoPAPt = &self.workspace.AtoPAPt; //mapping from input matrix entries to triuA
Expand All @@ -125,6 +148,8 @@ where
}
}

/// Update a subset of the values of the matrix to be (re)factored. See [`refactor`](crate::qdldl::QDLDLFactorisation::refactor)
///
pub fn scale_values(&mut self, indices: &[usize], scale: T) {
let nzval = &mut self.workspace.triuA.nzval; // post perm internal data
let AtoPAPt = &self.workspace.AtoPAPt; //mapping from input matrix entries to triuA
Expand All @@ -134,29 +159,43 @@ where
}
}

/// Shifts a subset of the values of the matrix to be (re)factored. The
/// values are offset by `offset`, with `signs` a vector of +/- 1 values
/// indicating the direction of shifts. See [`refactor`](crate::qdldl::QDLDLFactorisation::refactor)
///
pub fn offset_values(&mut self, indices: &[usize], offset: T, signs: &[i8]) {
assert_eq!(indices.len(), signs.len());

let nzval = &mut self.workspace.triuA.nzval; // post perm internal data
let AtoPAPt = &self.workspace.AtoPAPt; //mapping from input matrix entries to triuA

for (&idx, &sign) in zip(indices, signs) {
let sign: T = T::from_i8(sign).unwrap();
nzval[AtoPAPt[idx]] += offset * sign;
match sign.signum() {
1 => {
nzval[AtoPAPt[idx]] += offset;

Check warning on line 175 in src/qdldl/qdldl.rs

View check run for this annotation

Codecov / codecov/patch

src/qdldl/qdldl.rs#L173-L175

Added lines #L173 - L175 were not covered by tests
}
-1 => {
nzval[AtoPAPt[idx]] -= offset;

Check warning on line 178 in src/qdldl/qdldl.rs

View check run for this annotation

Codecov / codecov/patch

src/qdldl/qdldl.rs#L177-L178

Added lines #L177 - L178 were not covered by tests
}
_ => {}

Check warning on line 180 in src/qdldl/qdldl.rs

View check run for this annotation

Codecov / codecov/patch

src/qdldl/qdldl.rs#L180

Added line #L180 was not covered by tests
}
}
}

/// Refactor a matrix after its data has been modified. See [`update_values`](crate::qdldl::QDLDLFactorisation::update_values),
/// [`scale_values`](crate::qdldl::QDLDLFactorisation::scale_values) and [`offset_values`](crate::qdldl::QDLDLFactorisation::offset_values)
///
pub fn refactor(&mut self) -> Result<(), QDLDLError> {
// It never makes sense to call refactor for a logical
// factorization since it will always be the same. Calling
// this function implies that we want a numerical factorization
self.is_logical = false;
self.is_symbolic = false;
_factor(
&mut self.L,
&mut self.D,
&mut self.Dinv,
&mut self.workspace,
self.is_logical,
self.is_symbolic,
)
}
}
Expand Down Expand Up @@ -241,7 +280,7 @@ fn _qdldl_new<T: FloatT>(
D,
Dinv,
workspace,
is_logical: opts.logical,
is_symbolic: opts.logical,
})
}

Expand Down
2 changes: 2 additions & 0 deletions src/solver/core/solver.rs
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,9 @@ pub trait SolverJSONReadWrite<T>: Sized
where
T: FloatT,
{
/// write internal problem data to a JSON file
fn write_to_file(&self, file: &mut std::fs::File) -> Result<(), std::io::Error>;
/// load problem data from a JSON file previously saved using [`write_to_file`](self::SolverJSONReadWrite::write_to_file)
fn read_from_file(
file: &mut std::fs::File,
settings: Option<DefaultSettings<T>>,
Expand Down
Loading

0 comments on commit d7ca8b4

Please sign in to comment.