From c08ceb0a68658911ea17ec5eaa5f2878660902f7 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Wed, 26 Aug 2020 22:35:21 +0200 Subject: [PATCH 01/11] Bind SuiteSparse's CAMD fill-in reducing library --- Cargo.toml | 1 + .../suitesparse_camd_sys/.gitignore | 1 + .../suitesparse_camd_sys/Cargo.toml | 12 ++ .../suitesparse_camd_sys/src/lib.rs | 191 ++++++++++++++++++ 4 files changed, 205 insertions(+) create mode 100644 suitesparse_bindings/suitesparse_camd_sys/.gitignore create mode 100644 suitesparse_bindings/suitesparse_camd_sys/Cargo.toml create mode 100644 suitesparse_bindings/suitesparse_camd_sys/src/lib.rs diff --git a/Cargo.toml b/Cargo.toml index 5f4f443f..cbeb924a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -51,6 +51,7 @@ members = [ "sprs-ldl", "suitesparse_bindings/suitesparse_ldl_sys", "suitesparse_bindings/sprs_suitesparse_ldl", + "suitesparse_bindings/suitesparse_camd_sys", "sprs-rand", "sprs-benches", ] diff --git a/suitesparse_bindings/suitesparse_camd_sys/.gitignore b/suitesparse_bindings/suitesparse_camd_sys/.gitignore new file mode 100644 index 00000000..2f7896d1 --- /dev/null +++ b/suitesparse_bindings/suitesparse_camd_sys/.gitignore @@ -0,0 +1 @@ +target/ diff --git a/suitesparse_bindings/suitesparse_camd_sys/Cargo.toml b/suitesparse_bindings/suitesparse_camd_sys/Cargo.toml new file mode 100644 index 00000000..b16a14fe --- /dev/null +++ b/suitesparse_bindings/suitesparse_camd_sys/Cargo.toml @@ -0,0 +1,12 @@ +[package] +name = "suitesparse_camd_sys" +version = "0.1.0" +authors = ["Vincent Barrielle "] +edition = "2018" +description = "Raw bindings to SuiteSparse's CAMD algorithm" +license = "MIT OR Apache-2.0" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +libc = "0.2.74" diff --git a/suitesparse_bindings/suitesparse_camd_sys/src/lib.rs b/suitesparse_bindings/suitesparse_camd_sys/src/lib.rs new file mode 100644 index 00000000..62821c42 --- /dev/null +++ b/suitesparse_bindings/suitesparse_camd_sys/src/lib.rs @@ -0,0 +1,191 @@ +#[cfg(target_os = "windows")] +pub type SuiteSparseLong = libc::c_longlong; +#[cfg(not(target_os = "windows"))] +pub type SuiteSparseLong = libc::c_long; + +#[link(name = "camd")] +extern "C" { + /// Find a permutation matrix P, represented by the permutation indices + /// `p`, which reduces the fill-in of the symmetric sparse matrix A + /// (represented by its index pointer `ap` and indices pointer `ai`) + /// in Cholesky factorization (ie, the number of nonzeros of the Cholesky + /// factorization of P A P^T is less than for the Cholesky factorization + /// of A). + /// + /// If A is not symmetric, the ordering will be computed for A + A^T + /// + /// # Safety and constraints + /// - `n` is the number of rows and columns of the matrix A and should be + /// positive + /// - `ap` must be an array of length `n + 1`. + /// - `ai` must be an array of length `ap[n]`. The performance is better + /// if `ai[ap[i]..ap[i+1]]` is always sorted and without duplicates. + /// - `p` must be an array of length `n`. Its contents can be uninitialized. + /// - `control` must be an array of size `CAMD_CONTROL`. + /// - `info` must be an array of size `CAMD_INFO`. + /// - `constraint` must be either the null pointer, or an array of size `n`. + pub fn camd_order( + n: libc::c_int, + ap: *const libc::c_int, + ai: *const libc::c_int, + p: *mut libc::c_int, + control: *mut libc::c_double, + info: *mut libc::c_double, + constraint: *mut libc::c_int, + ) -> libc::c_int; + + /// Long version of `camd_order`, see its documentation. + pub fn camd_l_order( + n: SuiteSparseLong, + ap: *const SuiteSparseLong, + ai: *const SuiteSparseLong, + p: *mut SuiteSparseLong, + control: *mut libc::c_double, + info: *mut libc::c_double, + constraint: *mut SuiteSparseLong, + ) -> SuiteSparseLong; + + /// Checks if the matrix A represented by its index pointer array + /// `ap` and its indices array `ai` is suitable to pass to `camd_order`. + /// + /// Will return `CAMD_OK` if the matrix is suitable, and `CAMD_OK_BUT_JUMBLED` + /// if the matrix has unsorted or duplicate row indices in one or more columns. + /// The matrix must be square, ie `n_rows == n_cols`. + /// + /// Otherwise `CAMD_INVALID` will be returned. + pub fn camd_valid( + n_rows: libc::c_int, + n_cols: libc::c_int, + ap: *const libc::c_int, + ai: *const libc::c_int, + ) -> libc::c_int; + + /// Long version of `camd_valid`, see its documentation. + pub fn camd_l_valid( + n_rows: SuiteSparseLong, + n_cols: SuiteSparseLong, + ap: *const SuiteSparseLong, + ai: *const SuiteSparseLong, + ) -> SuiteSparseLong; + + /// Check if the array `constraint`, of size `n`, is valid as input + /// to `camd_order`. Returns `1` if valid, `0` otherwise. + pub fn camd_cvalid(n: libc::c_int, constraint: *const libc::c_int); + + /// Long version of `camd_cvalid`, see its documentation. + pub fn camd_l_cvalid( + n: SuiteSparseLong, + constraint: *const SuiteSparseLong, + ); + + /// Fill the `control` array of size `CAMD_CONTROL` with default values + pub fn camd_defaults(control: *mut libc::c_double); + + /// Fill the `control` array of size `CAMD_CONTROL` with default values + pub fn camd_l_defaults(control: *mut libc::c_double); + + /// Pretty print the `control` array of size `CAMD_CONTROL` + pub fn camd_control(control: *const libc::c_double); + + /// Pretty print the `control` array of size `CAMD_CONTROL` + pub fn camd_l_control(control: *const libc::c_double); + + /// Pretty print the `info` array of size `CAMD_INFO` + pub fn camd_info(info: *const libc::c_double); + + /// Pretty print the `info` array of size `CAMD_INFO` + pub fn camd_l_info(info: *const libc::c_double); +} + +pub const CAMD_CONTROL: usize = 5; +pub const CAMD_INFO: usize = 20; + +pub const CAMD_DENSE: usize = 0; +pub const CAMD_AGGRESSIVE: usize = 1; + +pub const CAMD_STATUS: usize = 0; +pub const CAMD_N: usize = 1; +pub const CAMD_NZ: usize = 2; +pub const CAMD_SYMMETRY: usize = 3; +pub const CAMD_NZDIAG: usize = 4; +pub const CAMD_NZ_A_PLUS_AT: usize = 5; +pub const CAMD_NDENSE: usize = 6; +pub const CAMD_MEMORY: usize = 7; +pub const CAMD_NCMPA: usize = 8; +pub const CAMD_LNZ: usize = 9; +pub const CAMD_NDIV: usize = 10; +pub const CAMD_NMULTSUBS_LDL: usize = 11; +pub const CAMD_NMULTSUBS_LU: usize = 12; +pub const CAMD_DMAX: usize = 13; + +pub const CAMD_OK: isize = 0; +pub const CAMD_OUT_OF_MEMORY: isize = -1; +pub const CAMD_INVALID: isize = -2; +pub const CAMD_OK_BUT_JUMBLED: isize = 1; + +#[cfg(test)] +mod tests { + #[test] + fn camd_valid() { + // | 0 1 3 | + // | 1 1 0 | + // | 3 0 0 | + let n: libc::c_int = 3; + let ap = &[0, 2, 4, 5]; + let ai = &[1, 2, 0, 1, 0]; + let valid = + unsafe { super::camd_valid(n, n, ap.as_ptr(), ai.as_ptr()) }; + assert_eq!(valid, super::CAMD_OK as libc::c_int); + let ai = &[2, 1, 0, 1, 0]; + let valid = + unsafe { super::camd_valid(n, n, ap.as_ptr(), ai.as_ptr()) }; + assert_eq!(valid, super::CAMD_OK_BUT_JUMBLED as libc::c_int); + let ai = &[1, 2, 0, 1, 1]; + let valid = + unsafe { super::camd_valid(n, n, ap.as_ptr(), ai.as_ptr()) }; + assert_eq!(valid, super::CAMD_OK as libc::c_int); + let valid = + unsafe { super::camd_valid(n, n + 1, ap.as_ptr(), ai.as_ptr()) }; + assert_eq!(valid, super::CAMD_INVALID as libc::c_int); + let valid = unsafe { + super::camd_valid(n + 1, n + 1, ap.as_ptr(), ai.as_ptr()) + }; + assert_eq!(valid, super::CAMD_INVALID as libc::c_int); + + // long version test, only test once as we now have the behavior tested + // (we only want tot test the binding is correct now). + use super::SuiteSparseLong as Long; + let n: Long = 3; + let ap: &[Long] = &[0, 2, 4, 5]; + let ai: &[Long] = &[1, 2, 0, 1, 0]; + let valid = + unsafe { super::camd_l_valid(n, n, ap.as_ptr(), ai.as_ptr()) }; + assert_eq!(valid, super::CAMD_OK as Long); + } + + #[test] + fn camd_order() { + // | 0 1 3 | + // | 1 1 0 | + // | 3 0 0 | + let n: libc::c_int = 3; + let ap = &[0, 2, 4, 5]; + let ai = &[1, 2, 0, 1, 0]; + let mut perm = [0; 3]; + let mut control = [0.; super::CAMD_CONTROL]; + let mut info = [0.; super::CAMD_INFO]; + let constraint: *const libc::c_int = std::ptr::null(); + let res = unsafe { + super::camd_order( + n, + ap.as_ptr(), + ai.as_ptr(), + perm.as_mut_ptr(), + control.as_mut_ptr(), + info.as_mut_ptr(), + constraint as *mut libc::c_int, + ) + }; + assert_eq!(res, super::CAMD_OK as libc::c_int); + } +} From f1f7ab6ea53badc5a82ea5789f688beb9191ba8c Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Wed, 2 Sep 2020 23:02:15 +0200 Subject: [PATCH 02/11] Add structure only sparse matrix type This type will make it easier to write symbolic algorithms without worrying about the scalar type. --- src/lib.rs | 7 ++++--- src/sparse.rs | 14 ++++++++++---- src/sparse/csmat.rs | 22 ++++++++++++++++++++++ 3 files changed, 36 insertions(+), 7 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index e62526cf..b3cfbd9e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -93,9 +93,10 @@ pub use crate::sparse::{ csmat::CsIter, csmat::OuterIterator, csmat::OuterIteratorMut, csmat::OuterIteratorPerm, kronecker::kronecker_product, CsMat, CsMatBase, CsMatI, CsMatVecView, CsMatView, CsMatViewI, CsMatViewMut, CsMatViewMutI, - CsVec, CsVecBase, CsVecI, CsVecView, CsVecViewI, CsVecViewMut, - CsVecViewMutI, SparseMat, TriMat, TriMatBase, TriMatI, TriMatIter, - TriMatView, TriMatViewI, TriMatViewMut, TriMatViewMutI, + CsStructure, CsStructureI, CsStructureView, CsStructureViewI, CsVec, + CsVecBase, CsVecI, CsVecView, CsVecViewI, CsVecViewMut, CsVecViewMutI, + SparseMat, TriMat, TriMatBase, TriMatI, TriMatIter, TriMatView, + TriMatViewI, TriMatViewMut, TriMatViewMutI, }; pub use crate::sparse::symmetric::is_symmetric; diff --git a/src/sparse.rs b/src/sparse.rs index e175b9a2..7b4ddf7f 100644 --- a/src/sparse.rs +++ b/src/sparse.rs @@ -111,6 +111,11 @@ pub type CsMatViewMut<'a, N> = CsMatViewMutI<'a, N, usize>; // FIXME: a fixed size array would be better, but no Deref impl pub type CsMatVecView<'a, N> = CsMatVecView_<'a, N, usize>; +pub type CsStructureViewI<'a, I, Iptr = I> = CsMatViewI<'a, (), I, Iptr>; +pub type CsStructureView<'a> = CsStructureViewI<'a, usize>; +pub type CsStructureI = CsMatI<(), I, Iptr>; +pub type CsStructure = CsStructureI; + /// A sparse vector, storing the indices of its non-zero data. /// /// A `CsVec` represents a sparse vector by storing a sorted `indices()` array @@ -235,10 +240,11 @@ pub struct TriMatIter { mod prelude { pub use super::{ CsMat, CsMatBase, CsMatI, CsMatVecView, CsMatVecView_, CsMatView, - CsMatViewI, CsMatViewMut, CsMatViewMutI, CsVec, CsVecBase, CsVecI, - CsVecView, CsVecViewI, CsVecViewMut, CsVecViewMutI, SparseMat, TriMat, - TriMatBase, TriMatI, TriMatIter, TriMatView, TriMatViewI, - TriMatViewMut, TriMatViewMutI, + CsMatViewI, CsMatViewMut, CsMatViewMutI, CsStructure, CsStructureI, + CsStructureView, CsStructureViewI, CsVec, CsVecBase, CsVecI, CsVecView, + CsVecViewI, CsVecViewMut, CsVecViewMutI, SparseMat, TriMat, TriMatBase, + TriMatI, TriMatIter, TriMatView, TriMatViewI, TriMatViewMut, + TriMatViewMutI, }; } diff --git a/src/sparse/csmat.rs b/src/sparse/csmat.rs index dad2993b..b268d220 100644 --- a/src/sparse/csmat.rs +++ b/src/sparse/csmat.rs @@ -1127,6 +1127,28 @@ where } } + pub fn structure_view(&self) -> CsStructureViewI { + // Safety: std::slice::from_raw_parts requires its passed + // pointer to be valid for the whole length of the slice. We have a + // zero-sized type, so the length is zero, and since we cast + // a non-null pointer, the pointer is valid as all pointers to zero-sized + // types are valid if they are not null. + let zst_data = unsafe { + std::slice::from_raw_parts( + self.data.as_ptr() as *const (), + self.data.len(), + ) + }; + CsStructureViewI { + storage: self.storage, + nrows: self.nrows, + ncols: self.ncols, + indptr: &self.indptr[..], + indices: &self.indices[..], + data: zst_data, + } + } + pub fn to_dense(&self) -> Array where N: Clone + Zero, From 8e0b5c279e787a407ec7b879fcd1917d90588e10 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Wed, 2 Sep 2020 23:04:48 +0200 Subject: [PATCH 03/11] Expose integer type used when binding SuiteSparse's CAMD --- .../suitesparse_camd_sys/src/lib.rs | 50 +++++++++++-------- 1 file changed, 28 insertions(+), 22 deletions(-) diff --git a/suitesparse_bindings/suitesparse_camd_sys/src/lib.rs b/suitesparse_bindings/suitesparse_camd_sys/src/lib.rs index 62821c42..31662b0f 100644 --- a/suitesparse_bindings/suitesparse_camd_sys/src/lib.rs +++ b/suitesparse_bindings/suitesparse_camd_sys/src/lib.rs @@ -3,6 +3,8 @@ pub type SuiteSparseLong = libc::c_longlong; #[cfg(not(target_os = "windows"))] pub type SuiteSparseLong = libc::c_long; +pub type SuiteSparseInt = libc::c_int; + #[link(name = "camd")] extern "C" { /// Find a permutation matrix P, represented by the permutation indices @@ -25,14 +27,14 @@ extern "C" { /// - `info` must be an array of size `CAMD_INFO`. /// - `constraint` must be either the null pointer, or an array of size `n`. pub fn camd_order( - n: libc::c_int, - ap: *const libc::c_int, - ai: *const libc::c_int, - p: *mut libc::c_int, + n: SuiteSparseInt, + ap: *const SuiteSparseInt, + ai: *const SuiteSparseInt, + p: *mut SuiteSparseInt, control: *mut libc::c_double, info: *mut libc::c_double, - constraint: *mut libc::c_int, - ) -> libc::c_int; + constraint: *mut SuiteSparseInt, + ) -> SuiteSparseInt; /// Long version of `camd_order`, see its documentation. pub fn camd_l_order( @@ -54,11 +56,11 @@ extern "C" { /// /// Otherwise `CAMD_INVALID` will be returned. pub fn camd_valid( - n_rows: libc::c_int, - n_cols: libc::c_int, - ap: *const libc::c_int, - ai: *const libc::c_int, - ) -> libc::c_int; + n_rows: SuiteSparseInt, + n_cols: SuiteSparseInt, + ap: *const SuiteSparseInt, + ai: *const SuiteSparseInt, + ) -> SuiteSparseInt; /// Long version of `camd_valid`, see its documentation. pub fn camd_l_valid( @@ -70,7 +72,7 @@ extern "C" { /// Check if the array `constraint`, of size `n`, is valid as input /// to `camd_order`. Returns `1` if valid, `0` otherwise. - pub fn camd_cvalid(n: libc::c_int, constraint: *const libc::c_int); + pub fn camd_cvalid(n: SuiteSparseInt, constraint: *const SuiteSparseInt); /// Long version of `camd_cvalid`, see its documentation. pub fn camd_l_cvalid( @@ -125,32 +127,33 @@ pub const CAMD_OK_BUT_JUMBLED: isize = 1; #[cfg(test)] mod tests { + use super::SuiteSparseInt; #[test] fn camd_valid() { // | 0 1 3 | // | 1 1 0 | // | 3 0 0 | - let n: libc::c_int = 3; + let n: SuiteSparseInt = 3; let ap = &[0, 2, 4, 5]; let ai = &[1, 2, 0, 1, 0]; let valid = unsafe { super::camd_valid(n, n, ap.as_ptr(), ai.as_ptr()) }; - assert_eq!(valid, super::CAMD_OK as libc::c_int); + assert_eq!(valid, super::CAMD_OK as SuiteSparseInt); let ai = &[2, 1, 0, 1, 0]; let valid = unsafe { super::camd_valid(n, n, ap.as_ptr(), ai.as_ptr()) }; - assert_eq!(valid, super::CAMD_OK_BUT_JUMBLED as libc::c_int); + assert_eq!(valid, super::CAMD_OK_BUT_JUMBLED as SuiteSparseInt); let ai = &[1, 2, 0, 1, 1]; let valid = unsafe { super::camd_valid(n, n, ap.as_ptr(), ai.as_ptr()) }; - assert_eq!(valid, super::CAMD_OK as libc::c_int); + assert_eq!(valid, super::CAMD_OK as SuiteSparseInt); let valid = unsafe { super::camd_valid(n, n + 1, ap.as_ptr(), ai.as_ptr()) }; - assert_eq!(valid, super::CAMD_INVALID as libc::c_int); + assert_eq!(valid, super::CAMD_INVALID as SuiteSparseInt); let valid = unsafe { super::camd_valid(n + 1, n + 1, ap.as_ptr(), ai.as_ptr()) }; - assert_eq!(valid, super::CAMD_INVALID as libc::c_int); + assert_eq!(valid, super::CAMD_INVALID as SuiteSparseInt); // long version test, only test once as we now have the behavior tested // (we only want tot test the binding is correct now). @@ -168,13 +171,16 @@ mod tests { // | 0 1 3 | // | 1 1 0 | // | 3 0 0 | - let n: libc::c_int = 3; + let n: SuiteSparseInt = 3; let ap = &[0, 2, 4, 5]; let ai = &[1, 2, 0, 1, 0]; let mut perm = [0; 3]; let mut control = [0.; super::CAMD_CONTROL]; let mut info = [0.; super::CAMD_INFO]; - let constraint: *const libc::c_int = std::ptr::null(); + let constraint: *const SuiteSparseInt = std::ptr::null(); + unsafe { + super::camd_defaults(control.as_mut_ptr()); + } let res = unsafe { super::camd_order( n, @@ -183,9 +189,9 @@ mod tests { perm.as_mut_ptr(), control.as_mut_ptr(), info.as_mut_ptr(), - constraint as *mut libc::c_int, + constraint as *mut SuiteSparseInt, ) }; - assert_eq!(res, super::CAMD_OK as libc::c_int); + assert_eq!(res, super::CAMD_OK as SuiteSparseInt); } } From f17893d1372cfbe5c338d9b6ccb8c9fad2cf6d75 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Wed, 2 Sep 2020 23:07:05 +0200 Subject: [PATCH 04/11] Add a high-level binding for SuiteSparse's CAMD --- Cargo.toml | 1 + .../sprs_suitesparse_camd/.gitignore | 1 + .../sprs_suitesparse_camd/Cargo.toml | 17 ++++ .../sprs_suitesparse_camd/src/lib.rs | 88 +++++++++++++++++++ 4 files changed, 107 insertions(+) create mode 100644 suitesparse_bindings/sprs_suitesparse_camd/.gitignore create mode 100644 suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml create mode 100644 suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs diff --git a/Cargo.toml b/Cargo.toml index cbeb924a..04435f81 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -52,6 +52,7 @@ members = [ "suitesparse_bindings/suitesparse_ldl_sys", "suitesparse_bindings/sprs_suitesparse_ldl", "suitesparse_bindings/suitesparse_camd_sys", + "suitesparse_bindings/sprs_suitesparse_camd", "sprs-rand", "sprs-benches", ] diff --git a/suitesparse_bindings/sprs_suitesparse_camd/.gitignore b/suitesparse_bindings/sprs_suitesparse_camd/.gitignore new file mode 100644 index 00000000..2f7896d1 --- /dev/null +++ b/suitesparse_bindings/sprs_suitesparse_camd/.gitignore @@ -0,0 +1 @@ +target/ diff --git a/suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml b/suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml new file mode 100644 index 00000000..80acb8e1 --- /dev/null +++ b/suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml @@ -0,0 +1,17 @@ +[package] +name = "sprs_suitesparse_camd" +description = "sprs bindings to the suitesparse camd fill-in reducting ordering" +version = "0.1.0" +authors = ["Vincent Barrielle "] +edition = "2018" +keywords = ["sparse", "matrix", "fill-in", "permutation", "suitesparse"] + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies.sprs] +version = "0.8.0" +path = "../.." + +[dependencies.suitesparse_camd_sys] +path = "../suitesparse_camd_sys/" +version = "0.1.0" diff --git a/suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs b/suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs new file mode 100644 index 00000000..6bada7e2 --- /dev/null +++ b/suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs @@ -0,0 +1,88 @@ +use sprs::errors::SprsError; +use sprs::{CsStructureI, CsStructureViewI, PermOwnedI, SpIndex}; +use suitesparse_camd_sys::*; + +/// Find a permutation matrix P which reduces the fill-in of the square +/// sparse matrix `mat` in Cholesky factorization (ie, the number of nonzeros +/// of the Cholesky factorization of P A P^T is less than for the Cholesky +/// factorization of A). +/// +/// If A is not symmetric, the ordering will be computed for A + A^T +/// +/// # Errors +/// +/// This function will error if the passed matrix is not square. +pub fn try_camd( + mat: CsStructureViewI, +) -> Result, SprsError> +where + I: SpIndex, + Iptr: SpIndex, +{ + let n = mat.rows(); + if n != mat.cols() { + return Err(SprsError::IllegalArguments( + "Input to camd must be square", + )); + } + let mut perm: Vec = vec![0; n]; + let mut control = [0.; CAMD_CONTROL]; + let mut info = [0.; CAMD_INFO]; + let constraint: *const SuiteSparseInt = std::ptr::null(); + // TODO use the long version when required by the size of the matrix + let mat: CsStructureI = + mat.to_other_types(); + let res = unsafe { + camd_order( + n as SuiteSparseInt, + mat.indptr().as_ptr(), + mat.indices().as_ptr(), + perm.as_mut_ptr(), + control.as_mut_ptr(), + info.as_mut_ptr(), + constraint as *mut SuiteSparseInt, + ) + }; + // CsMat invariants guarantee sorted and non duplicate indices so this + // should not happen. + if res as isize != CAMD_OK { + return Err(SprsError::NonSortedIndices); + } + let perm = perm.iter().map(|&i| I::from_usize(i as usize)).collect(); + Ok(PermOwnedI::new(perm)) +} + +/// Find a permutation matrix P which reduces the fill-in of the square +/// sparse matrix `mat` in Cholesky factorization (ie, the number of nonzeros +/// of the Cholesky factorization of P A P^T is less than for the Cholesky +/// factorization of A). +/// +/// If A is not symmetric, the ordering will be computed for A + A^T +/// +/// # Panics +/// +/// This function will panic if the passed matrix is not square. +pub fn camd(mat: CsStructureViewI) -> PermOwnedI +where + I: SpIndex, + Iptr: SpIndex, +{ + try_camd(mat).unwrap() +} + +#[cfg(test)] +mod tests { + use sprs::CsMatI; + + #[test] + fn try_camd() { + let mat = CsMatI::new_csc( + (4, 4), + vec![0, 2, 4, 6, 8], + vec![0, 3, 1, 2, 1, 2, 0, 3], + vec![1., 2., 21., 6., 6., 2., 2., 8.], + ); + let res = super::try_camd(mat.structure_view()); + assert!(res.is_ok()); + } +} From 89ffe05aefd00f2289ce5392d5fc6b0daf65c674 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Mon, 7 Sep 2020 13:36:33 +0000 Subject: [PATCH 05/11] Call camd_l_long when matrix size is longer than i32::MAX --- .../sprs_suitesparse_camd/src/lib.rs | 54 +++++++++++++------ 1 file changed, 37 insertions(+), 17 deletions(-) diff --git a/suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs b/suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs index 6bada7e2..5f5c4d84 100644 --- a/suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs +++ b/suitesparse_bindings/sprs_suitesparse_camd/src/lib.rs @@ -25,30 +25,50 @@ where "Input to camd must be square", )); } - let mut perm: Vec = vec![0; n]; let mut control = [0.; CAMD_CONTROL]; let mut info = [0.; CAMD_INFO]; - let constraint: *const SuiteSparseInt = std::ptr::null(); - // TODO use the long version when required by the size of the matrix - let mat: CsStructureI = - mat.to_other_types(); - let res = unsafe { - camd_order( - n as SuiteSparseInt, - mat.indptr().as_ptr(), - mat.indices().as_ptr(), - perm.as_mut_ptr(), - control.as_mut_ptr(), - info.as_mut_ptr(), - constraint as *mut SuiteSparseInt, - ) + let (camd_res, perm) = if n <= SuiteSparseInt::MAX as usize { + let constraint: *const SuiteSparseInt = std::ptr::null(); + let mat: CsStructureI = + mat.to_other_types(); + let mut perm: Vec = vec![0; n]; + let camd_res = unsafe { + camd_order( + n as SuiteSparseInt, + mat.indptr().as_ptr(), + mat.indices().as_ptr(), + perm.as_mut_ptr(), + control.as_mut_ptr(), + info.as_mut_ptr(), + constraint as *mut SuiteSparseInt, + ) as isize + }; + let perm = perm.iter().map(|&i| I::from_usize(i as usize)).collect(); + (camd_res, perm) + } else { + let constraint: *const SuiteSparseLong = std::ptr::null(); + let mat: CsStructureI = + mat.to_other_types(); + let mut perm: Vec = vec![0; n]; + let camd_res = unsafe { + camd_l_order( + n as SuiteSparseLong, + mat.indptr().as_ptr(), + mat.indices().as_ptr(), + perm.as_mut_ptr(), + control.as_mut_ptr(), + info.as_mut_ptr(), + constraint as *mut SuiteSparseLong, + ) + } as isize; + let perm = perm.iter().map(|&i| I::from_usize(i as usize)).collect(); + (camd_res, perm) }; // CsMat invariants guarantee sorted and non duplicate indices so this // should not happen. - if res as isize != CAMD_OK { + if camd_res != CAMD_OK { return Err(SprsError::NonSortedIndices); } - let perm = perm.iter().map(|&i| I::from_usize(i as usize)).collect(); Ok(PermOwnedI::new(perm)) } From 3dbab0b69f02c2f00b7d7fd2f8668fd20f37db95 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Thu, 10 Sep 2020 21:06:45 +0000 Subject: [PATCH 06/11] Make fill in reduction enum non exhaustive This will avoid version churn --- sprs-ldl/src/lib.rs | 5 +++++ src/lib.rs | 1 + 2 files changed, 6 insertions(+) diff --git a/sprs-ldl/src/lib.rs b/sprs-ldl/src/lib.rs index fa6f709a..7db7bd87 100644 --- a/sprs-ldl/src/lib.rs +++ b/sprs-ldl/src/lib.rs @@ -145,6 +145,11 @@ impl Ldl { FillInReduction::ReverseCuthillMcKee => { sprs::linalg::reverse_cuthill_mckee(mat).perm } + _ => { + unreachable!( + "Unhandled method, report a bug at https://github.com/vbarrielle/sprs/issues/199" + ) + } } } diff --git a/src/lib.rs b/src/lib.rs index b3cfbd9e..0490de45 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -153,6 +153,7 @@ pub use PermutationCheck::*; /// The different kinds of fill-in-reduction algorithms supported by sprs #[derive(Copy, Clone, Eq, PartialEq, Debug)] +#[non_exhaustive] pub enum FillInReduction { NoReduction, ReverseCuthillMcKee, From d74e9d68db7e0b47ed697fd14ce52bcbf3054afb Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Thu, 10 Sep 2020 21:08:07 +0000 Subject: [PATCH 07/11] Bump versions to reflect breaking changes Though this commit will not be the one where I publish, this will come later, but I'm sure I won't forget about this breaking change now. --- Cargo.toml | 2 +- changelog.rst | 4 ++++ sprs-benches/Cargo.toml | 4 ++-- sprs-ldl/Cargo.toml | 6 +++--- sprs-rand/Cargo.toml | 4 ++-- suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml | 2 +- suitesparse_bindings/sprs_suitesparse_ldl/Cargo.toml | 4 ++-- 7 files changed, 15 insertions(+), 11 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 04435f81..7134a1ae 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -2,7 +2,7 @@ name = "sprs" description = "A sparse matrix library" -version = "0.8.1" +version = "0.9.0" authors = ["Vincent Barrielle "] edition = "2018" diff --git a/changelog.rst b/changelog.rst index a9bfe01d..8f520fb2 100644 --- a/changelog.rst +++ b/changelog.rst @@ -2,6 +2,10 @@ Changelog ========= +- 0.9.0 + - Make FillInReduction enum non exhaustive to prevent excessive breakage + when new algorithms are implemented. **breaking change** + - Make rayon optional **breaking change** - 0.8.1 - Expose the ``num_kinds`` module to allow generic usage of matrix market serialization functions in client crates diff --git a/sprs-benches/Cargo.toml b/sprs-benches/Cargo.toml index 24d59c58..8dd8bf73 100644 --- a/sprs-benches/Cargo.toml +++ b/sprs-benches/Cargo.toml @@ -5,11 +5,11 @@ authors = ["Vincent Barrielle "] edition = "2018" [dependencies.sprs] -version = "0.8.0" +version = "0.9.0" path = ".." [dependencies.sprs-rand] -version = "0.1.0" +version = "0.2.0" path = "../sprs-rand" [dependencies.plotters] diff --git a/sprs-ldl/Cargo.toml b/sprs-ldl/Cargo.toml index a8a1da04..748bf341 100644 --- a/sprs-ldl/Cargo.toml +++ b/sprs-ldl/Cargo.toml @@ -2,7 +2,7 @@ name = "sprs-ldl" description = "Sparse cholesky factorization" -version = "0.6.1" +version = "0.7.0" authors = ["Vincent Barrielle"] edition = "2018" @@ -18,10 +18,10 @@ num-traits = "0.1.32" [dependencies.sprs] -version = "0.8.0" +version = "0.9.0" path = ".." [dependencies.sprs_suitesparse_ldl] -version = "0.4.0" +version = "0.5.0" path = "../suitesparse_bindings/sprs_suitesparse_ldl" optional = true diff --git a/sprs-rand/Cargo.toml b/sprs-rand/Cargo.toml index 21c5e815..bdc8916e 100644 --- a/sprs-rand/Cargo.toml +++ b/sprs-rand/Cargo.toml @@ -1,7 +1,7 @@ [package] name = "sprs-rand" description = "Random sparse matrix generation" -version = "0.1.0" +version = "0.2.0" authors = ["Vincent Barrielle "] license = "MIT OR Apache-2.0" edition = "2018" @@ -13,5 +13,5 @@ rand_distr = "0.2.2" rand_pcg = "0.2.1" [dependencies.sprs] -version = "0.8.0" +version = "0.9.0" path = ".." diff --git a/suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml b/suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml index 80acb8e1..9d59ddc2 100644 --- a/suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml +++ b/suitesparse_bindings/sprs_suitesparse_camd/Cargo.toml @@ -9,7 +9,7 @@ keywords = ["sparse", "matrix", "fill-in", "permutation", "suitesparse"] # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [dependencies.sprs] -version = "0.8.0" +version = "0.9.0" path = "../.." [dependencies.suitesparse_camd_sys] diff --git a/suitesparse_bindings/sprs_suitesparse_ldl/Cargo.toml b/suitesparse_bindings/sprs_suitesparse_ldl/Cargo.toml index 9e7aefb4..348281df 100644 --- a/suitesparse_bindings/sprs_suitesparse_ldl/Cargo.toml +++ b/suitesparse_bindings/sprs_suitesparse_ldl/Cargo.toml @@ -1,7 +1,7 @@ [package] name = "sprs_suitesparse_ldl" description = "sprs bindings to the suitesparse ldl solver" -version = "0.4.0" +version = "0.5.0" authors = ["Vincent Barrielle "] license = "MIT OR Apache-2.0" repository = "https://github.com/vbarrielle/sprs" @@ -17,5 +17,5 @@ path = "../suitesparse_ldl_sys/" version = "0.2.0" [dependencies.sprs] -version = "0.8.0" +version = "0.9.0" path = "../.." From 17053dfe3d01b6e2589d595aa60f60959fa0e933 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Thu, 10 Sep 2020 21:09:23 +0000 Subject: [PATCH 08/11] Relax bounds in Ldl::perm Since we can now only use the structure, we can actually remove constraints on N. --- sprs-ldl/src/lib.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sprs-ldl/src/lib.rs b/sprs-ldl/src/lib.rs index 7db7bd87..ec48af6f 100644 --- a/sprs-ldl/src/lib.rs +++ b/sprs-ldl/src/lib.rs @@ -138,12 +138,11 @@ impl Ldl { pub fn perm(&self, mat: CsMatViewI) -> PermOwnedI where I: SpIndex, - N: Copy + PartialEq, { match self.fill_red_method { FillInReduction::NoReduction => PermOwnedI::identity(mat.rows()), FillInReduction::ReverseCuthillMcKee => { - sprs::linalg::reverse_cuthill_mckee(mat).perm + sprs::linalg::reverse_cuthill_mckee(mat.structure_view()).perm } _ => { unreachable!( From 42e84f23476c580f3f4d4206f46089644dda5615 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Thu, 10 Sep 2020 21:10:28 +0000 Subject: [PATCH 09/11] Integrate SuiteSparse's CAMD in Ldl builder struct --- sprs-ldl/Cargo.toml | 8 ++++++++ sprs-ldl/src/lib.rs | 50 +++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 1 + 3 files changed, 59 insertions(+) diff --git a/sprs-ldl/Cargo.toml b/sprs-ldl/Cargo.toml index 748bf341..2f3667c9 100644 --- a/sprs-ldl/Cargo.toml +++ b/sprs-ldl/Cargo.toml @@ -25,3 +25,11 @@ path = ".." version = "0.5.0" path = "../suitesparse_bindings/sprs_suitesparse_ldl" optional = true + +[dependencies.sprs_suitesparse_camd] +version = "0.1.0" +path = "../suitesparse_bindings/sprs_suitesparse_camd" +optional = true + +[dev-dependencies] +ndarray = ">=0.11.0,<0.14" diff --git a/sprs-ldl/src/lib.rs b/sprs-ldl/src/lib.rs index ec48af6f..5db51096 100644 --- a/sprs-ldl/src/lib.rs +++ b/sprs-ldl/src/lib.rs @@ -144,6 +144,12 @@ impl Ldl { FillInReduction::ReverseCuthillMcKee => { sprs::linalg::reverse_cuthill_mckee(mat.structure_view()).perm } + FillInReduction::CAMDSuiteSparse => { + #[cfg(not(feature = "sprs_suitesparse_camd"))] + panic!("Unavailable without the `sprs_suitesparse_camd` feature"); + #[cfg(feature = "sprs_suitesparse_camd")] + sprs_suitesparse_camd::camd(mat.structure_view()) + } _ => { unreachable!( "Unhandled method, report a bug at https://github.com/vbarrielle/sprs/issues/199" @@ -865,4 +871,48 @@ mod test { let x = ldlt.solve(&b); assert_eq!(x, x0); } + + #[cfg(feature = "sprs_suitesparse_camd")] + #[test] + fn camd_ldl_solve() { + // 0 - A - 2 - 3 + // | \ | \ | / | + // 7 - 5 - 6 - 4 + // | / | / | \ | + // 8 - 9 - 1 - E + #[rustfmt::skip] + let triangles = ndarray::arr2( + &[[0, 7, 5], + [0, 5, 10], + [10, 5, 6], + [10, 6, 2], + [2, 6, 3], + [3, 6, 4], + [7, 8, 5], + [5, 8, 9], + [5, 9, 6], + [6, 9, 1], + [6, 1, 11], + [6, 11, 4]], + ); + let lap_mat = + sprs::special_mats::tri_mesh_graph_laplacian(12, triangles.view()); + let ldlt_camd = super::Ldl::new() + .check_symmetry(super::SymmetryCheck::DontCheckSymmetry) + .fill_in_reduction(super::FillInReduction::CAMDSuiteSparse) + .numeric(lap_mat.view()) + .unwrap(); + let ldlt_cuthill = super::Ldl::new() + .check_symmetry(super::SymmetryCheck::DontCheckSymmetry) + .fill_in_reduction(super::FillInReduction::ReverseCuthillMcKee) + .numeric(lap_mat.view()) + .unwrap(); + let ldlt_raw = super::Ldl::new() + .check_symmetry(super::SymmetryCheck::DontCheckSymmetry) + .fill_in_reduction(super::FillInReduction::NoReduction) + .numeric(lap_mat.view()) + .unwrap(); + assert!(ldlt_camd.nnz() < ldlt_raw.nnz()); + assert!(ldlt_camd.nnz() < ldlt_cuthill.nnz()); + } } diff --git a/src/lib.rs b/src/lib.rs index 0490de45..e567b13b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -157,6 +157,7 @@ pub use PermutationCheck::*; pub enum FillInReduction { NoReduction, ReverseCuthillMcKee, + CAMDSuiteSparse, } #[cfg(feature = "approx")] From 94ab762dcffb07b1dc63ca4520841429de0395fb Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Fri, 11 Sep 2020 08:41:18 +0000 Subject: [PATCH 10/11] Remove nightly-only checks from travis-ci These are no longer necessary as pyo3 builds on stable now. --- .travis.yml | 1 - .travis_nightly | 13 ------------- 2 files changed, 14 deletions(-) delete mode 100755 .travis_nightly diff --git a/.travis.yml b/.travis.yml index 75538da0..7d9549b4 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,7 +19,6 @@ script: - cargo test --release --all --verbose - cargo run --example heat - ./.travis_rustfmt - - ./.travis_nightly notifications: email: diff --git a/.travis_nightly b/.travis_nightly deleted file mode 100755 index 268de7ce..00000000 --- a/.travis_nightly +++ /dev/null @@ -1,13 +0,0 @@ -#!/bin/bash - -set -ev - -if [ -z ${WITH_NIGHTLY+x} ]; then - echo "Not on nightly channel: skipping nightly only parts." -else - cd sprs-benches - cargo build --features nightly - cd .. - cd sprs-ldl - cargo test --features sprs_suitesparse_ldl -fi From c17ea190efb7696a544c966b5cb335e2205ba981 Mon Sep 17 00:00:00 2001 From: Vincent Barrielle Date: Fri, 11 Sep 2020 09:01:18 +0000 Subject: [PATCH 11/11] Exclude suitesparse CAMD from ci builds other than ubuntu At some point I'll need to figure out a proper way to use suitesparse on platforms other than linux. --- .github/workflows/ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 899dbad2..8c843dd6 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -105,7 +105,7 @@ jobs: override: true - name: Build (exclude suitesparse) - run: cargo build --verbose --workspace --exclude suitesparse_ldl_sys --exclude sprs_suitesparse_ldl + run: cargo build --verbose --workspace --exclude suitesparse_ldl_sys --exclude sprs_suitesparse_ldl --exclude sprs_suitesparse_camd --exclude suitesparse_camd_sys if: matrix.os != 'ubuntu-18.04' - name: Build (all) @@ -113,7 +113,7 @@ jobs: if: matrix.os == 'ubuntu-18.04' - name: Test (exclude suitesparse) - run: cargo test --verbose --workspace --exclude suitesparse_ldl_sys --exclude sprs_suitesparse_ldl + run: cargo test --verbose --workspace --exclude suitesparse_ldl_sys --exclude sprs_suitesparse_ldl --exclude sprs_suitesparse_camd --exclude suitesparse_camd_sys if: matrix.os != 'ubuntu-18.04' - name: Test (all) @@ -121,7 +121,7 @@ jobs: if: matrix.os == 'ubuntu-18.04' - name: Test (suitesparse integration) - run: cargo test -p sprs-ldl --features sprs_suitesparse_ldl -Zpackage-features + run: cargo test -p sprs-ldl --features "sprs_suitesparse_ldl sprs_suitesparse_camd" -Zpackage-features if: matrix.os == 'ubuntu-18.04' && matrix.build == 'nightly' benches: