From ef3477f411278c3e97845f0002a4c7138b71a331 Mon Sep 17 00:00:00 2001 From: Andreas Longva Date: Tue, 19 Jan 2021 14:15:19 +0100 Subject: [PATCH] Remove Zero bound for transpose and impl SparsityPattern::transpose --- nalgebra-sparse/src/convert/impl_std_ops.rs | 4 +- nalgebra-sparse/src/convert/serial.rs | 99 ++++------------- nalgebra-sparse/src/cs.rs | 111 ++++++++++++++++++-- nalgebra-sparse/src/csc.rs | 4 +- nalgebra-sparse/src/csr.rs | 4 +- nalgebra-sparse/src/pattern.rs | 19 ++++ 6 files changed, 151 insertions(+), 90 deletions(-) diff --git a/nalgebra-sparse/src/convert/impl_std_ops.rs b/nalgebra-sparse/src/convert/impl_std_ops.rs index 0a3fbdf6..66f70408 100644 --- a/nalgebra-sparse/src/convert/impl_std_ops.rs +++ b/nalgebra-sparse/src/convert/impl_std_ops.rs @@ -107,7 +107,7 @@ impl<'a, T> From<&'a CscMatrix> for DMatrix impl<'a, T> From<&'a CscMatrix> for CsrMatrix where - T: Scalar + Zero + T: Scalar { fn from(matrix: &'a CscMatrix) -> Self { convert_csc_csr(matrix) @@ -116,7 +116,7 @@ impl<'a, T> From<&'a CscMatrix> for CsrMatrix impl<'a, T> From<&'a CsrMatrix> for CscMatrix where - T: Scalar + Zero + T: Scalar { fn from(matrix: &'a CsrMatrix) -> Self { convert_csr_csc(matrix) diff --git a/nalgebra-sparse/src/convert/serial.rs b/nalgebra-sparse/src/convert/serial.rs index f82020f4..4ddefcca 100644 --- a/nalgebra-sparse/src/convert/serial.rs +++ b/nalgebra-sparse/src/convert/serial.rs @@ -1,12 +1,15 @@ //! TODO -use crate::coo::CooMatrix; -use crate::csr::CsrMatrix; -use nalgebra::{DMatrix, Scalar, Matrix, Dim, ClosedAdd}; -use nalgebra::storage::Storage; +use std::ops::Add; + use num_traits::Zero; -use std::ops::{Add}; +use nalgebra::{ClosedAdd, Dim, DMatrix, Matrix, Scalar}; +use nalgebra::storage::Storage; + +use crate::coo::CooMatrix; +use crate::cs; use crate::csc::CscMatrix; +use crate::csr::CsrMatrix; /// TODO pub fn convert_dense_coo(dense: &Matrix) -> CooMatrix @@ -192,13 +195,13 @@ pub fn convert_dense_csc(dense: &Matrix) -> CscMatrix /// TODO pub fn convert_csr_csc(csr: &CsrMatrix) -> CscMatrix where - T: Scalar + Zero + T: Scalar { - let (offsets, indices, values) = transpose_cs(csr.nrows(), - csr.ncols(), - csr.row_offsets(), - csr.col_indices(), - csr.values()); + let (offsets, indices, values) = cs::transpose_cs(csr.nrows(), + csr.ncols(), + csr.row_offsets(), + csr.col_indices(), + csr.values()); // TODO: Avoid data validity check? CscMatrix::try_from_csc_data(csr.nrows(), csr.ncols(), offsets, indices, values) @@ -208,13 +211,13 @@ where /// TODO pub fn convert_csc_csr(csc: &CscMatrix) -> CsrMatrix where - T: Scalar + Zero + T: Scalar { - let (offsets, indices, values) = transpose_cs(csc.ncols(), - csc.nrows(), - csc.col_offsets(), - csc.row_indices(), - csc.values()); + let (offsets, indices, values) = cs::transpose_cs(csc.ncols(), + csc.nrows(), + csc.col_offsets(), + csc.row_indices(), + csc.values()); // TODO: Avoid data validity check? CsrMatrix::try_from_csr_data(csc.nrows(), csc.ncols(), offsets, indices, values) @@ -326,7 +329,7 @@ fn coo_to_unsorted_cs( major_offsets[*major_idx] += 1; } - convert_counts_to_offsets(major_offsets); + cs::convert_counts_to_offsets(major_offsets); { // TODO: Instead of allocating a whole new vector storing the current counts, @@ -377,66 +380,6 @@ fn sort_lane( apply_permutation(values_result, values, permutation); } -/// Transposes the compressed format. -/// -/// This means that major and minor roles are switched. This is used for converting between CSR -/// and CSC formats. -fn transpose_cs(major_dim: usize, - minor_dim: usize, - source_major_offsets: &[usize], - source_minor_indices: &[usize], - values: &[T]) - -> (Vec, Vec, Vec) -where - T: Scalar + Zero -{ - assert_eq!(source_major_offsets.len(), major_dim + 1); - assert_eq!(source_minor_indices.len(), values.len()); - let nnz = values.len(); - - // Count the number of occurences of each minor index - let mut minor_counts = vec![0; minor_dim]; - for minor_idx in source_minor_indices { - minor_counts[*minor_idx] += 1; - } - convert_counts_to_offsets(&mut minor_counts); - let mut target_offsets = minor_counts; - target_offsets.push(nnz); - let mut target_indices = vec![usize::MAX; nnz]; - let mut target_values = vec![T::zero(); nnz]; - - // Keep track of how many entries we have placed in each target major lane - let mut current_target_major_counts = vec![0; minor_dim]; - - for source_major_idx in 0 .. major_dim { - let source_lane_begin = source_major_offsets[source_major_idx]; - let source_lane_end = source_major_offsets[source_major_idx + 1]; - let source_lane_indices = &source_minor_indices[source_lane_begin .. source_lane_end]; - let source_lane_values = &values[source_lane_begin .. source_lane_end]; - - for (&source_minor_idx, val) in source_lane_indices.iter().zip(source_lane_values) { - // Compute the offset in the target data for this particular source entry - let target_lane_count = &mut current_target_major_counts[source_minor_idx]; - let entry_offset = target_offsets[source_minor_idx] + *target_lane_count; - target_indices[entry_offset] = source_major_idx; - target_values[entry_offset] = val.inlined_clone(); - *target_lane_count += 1; - } - } - - (target_offsets, target_indices, target_values) -} - -fn convert_counts_to_offsets(counts: &mut [usize]) { - // Convert the counts to an offset - let mut offset = 0; - for i_offset in counts.iter_mut() { - let count = *i_offset; - *i_offset = offset; - offset += count; - } -} - // TODO: Move this into `utils` or something? fn apply_permutation(out_slice: &mut [T], in_slice: &[T], permutation: &[usize]) { assert_eq!(out_slice.len(), in_slice.len()); diff --git a/nalgebra-sparse/src/cs.rs b/nalgebra-sparse/src/cs.rs index 634ce413..30d73007 100644 --- a/nalgebra-sparse/src/cs.rs +++ b/nalgebra-sparse/src/cs.rs @@ -1,12 +1,14 @@ -use crate::pattern::SparsityPattern; -use crate::{SparseEntry, SparseEntryMut}; - -use std::sync::Arc; -use std::ops::Range; use std::mem::replace; +use std::ops::Range; +use std::sync::Arc; + use num_traits::One; + use nalgebra::Scalar; +use crate::{SparseEntry, SparseEntryMut}; +use crate::pattern::SparsityPattern; + /// An abstract compressed matrix. /// /// For the time being, this is only used internally to share implementation between @@ -396,4 +398,101 @@ impl<'a, T> CsLaneMut<'a, T> { self.values, global_minor_index) } -} \ No newline at end of file +} + +/// Helper struct for working with uninitialized data in vectors. +/// TODO: This doesn't belong here. +struct UninitVec { + vec: Vec +} + +impl UninitVec { + pub fn from_len(len: usize) -> Self { + Self { + vec: Vec::with_capacity(len) + } + } + + /// Sets the element associated with the given index to the provided value. + /// + /// Must be called exactly once per index, otherwise results in undefined behavior. + pub unsafe fn set(&mut self, index: usize, value: T) { + self.vec.as_mut_ptr().add(index).write(value) + } + + /// Marks the vector data as initialized by returning a full vector. + /// + /// It is undefined behavior to call this function unless *all* elements have been written to + /// exactly once. + pub unsafe fn assume_init(mut self) -> Vec { + self.vec.set_len(self.vec.capacity()); + self.vec + } +} + +/// Transposes the compressed format. +/// +/// This means that major and minor roles are switched. This is used for converting between CSR +/// and CSC formats. +pub fn transpose_cs( + major_dim: usize, + minor_dim: usize, + source_major_offsets: &[usize], + source_minor_indices: &[usize], + values: &[T]) + -> (Vec, Vec, Vec) +where + T: Scalar +{ + assert_eq!(source_major_offsets.len(), major_dim + 1); + assert_eq!(source_minor_indices.len(), values.len()); + let nnz = values.len(); + + // Count the number of occurences of each minor index + let mut minor_counts = vec![0; minor_dim]; + for minor_idx in source_minor_indices { + minor_counts[*minor_idx] += 1; + } + convert_counts_to_offsets(&mut minor_counts); + let mut target_offsets = minor_counts; + target_offsets.push(nnz); + let mut target_indices = vec![usize::MAX; nnz]; + + // We have to use uninitialized storage, because we don't have any kind of "default" value + // available for `T`. Unfortunately this necessitates some small amount of unsafe code + let mut target_values = UninitVec::from_len(nnz); + + // Keep track of how many entries we have placed in each target major lane + let mut current_target_major_counts = vec![0; minor_dim]; + + for source_major_idx in 0 .. major_dim { + let source_lane_begin = source_major_offsets[source_major_idx]; + let source_lane_end = source_major_offsets[source_major_idx + 1]; + let source_lane_indices = &source_minor_indices[source_lane_begin .. source_lane_end]; + let source_lane_values = &values[source_lane_begin .. source_lane_end]; + + for (&source_minor_idx, val) in source_lane_indices.iter().zip(source_lane_values) { + // Compute the offset in the target data for this particular source entry + let target_lane_count = &mut current_target_major_counts[source_minor_idx]; + let entry_offset = target_offsets[source_minor_idx] + *target_lane_count; + target_indices[entry_offset] = source_major_idx; + unsafe { target_values.set(entry_offset, val.inlined_clone()); } + *target_lane_count += 1; + } + } + + // At this point, we should have written to each element in target_values exactly once, + // so initialization should be sound + let target_values = unsafe { target_values.assume_init() }; + (target_offsets, target_indices, target_values) +} + +pub fn convert_counts_to_offsets(counts: &mut [usize]) { + // Convert the counts to an offset + let mut offset = 0; + for i_offset in counts.iter_mut() { + let count = *i_offset; + *i_offset = offset; + offset += count; + } +} diff --git a/nalgebra-sparse/src/csc.rs b/nalgebra-sparse/src/csc.rs index 1d8b8970..c3d79508 100644 --- a/nalgebra-sparse/src/csc.rs +++ b/nalgebra-sparse/src/csc.rs @@ -7,7 +7,7 @@ use crate::cs::{CsMatrix, CsLane, CsLaneMut, CsLaneIter, CsLaneIterMut}; use std::sync::Arc; use std::slice::{IterMut, Iter}; -use num_traits::{Zero, One}; +use num_traits::{One}; use nalgebra::Scalar; /// A CSC representation of a sparse matrix. @@ -368,7 +368,7 @@ impl CscMatrix { impl CscMatrix where - T: Scalar + Zero + T: Scalar { /// Compute the transpose of the matrix. pub fn transpose(&self) -> CscMatrix { diff --git a/nalgebra-sparse/src/csr.rs b/nalgebra-sparse/src/csr.rs index 1f621c86..81f8191c 100644 --- a/nalgebra-sparse/src/csr.rs +++ b/nalgebra-sparse/src/csr.rs @@ -5,7 +5,7 @@ use crate::csc::CscMatrix; use crate::cs::{CsMatrix, CsLaneIterMut, CsLaneIter, CsLane, CsLaneMut}; use nalgebra::Scalar; -use num_traits::{Zero, One}; +use num_traits::{One}; use std::sync::Arc; use std::slice::{IterMut, Iter}; @@ -368,7 +368,7 @@ impl CsrMatrix { impl CsrMatrix where - T: Scalar + Zero + T: Scalar { /// Compute the transpose of the matrix. pub fn transpose(&self) -> CsrMatrix { diff --git a/nalgebra-sparse/src/pattern.rs b/nalgebra-sparse/src/pattern.rs index f95f23c8..6a0cf851 100644 --- a/nalgebra-sparse/src/pattern.rs +++ b/nalgebra-sparse/src/pattern.rs @@ -2,6 +2,7 @@ use crate::SparseFormatError; use std::fmt; use std::error::Error; +use crate::cs::transpose_cs; /// A representation of the sparsity pattern of a CSR or CSC matrix. /// @@ -204,6 +205,24 @@ impl SparsityPattern { pub fn disassemble(self) -> (Vec, Vec) { (self.major_offsets, self.minor_indices) } + + /// TODO + pub fn transpose(&self) -> Self { + // By using unit () values, we can use the same routines as for CSR/CSC matrices + let values = vec![(); self.nnz()]; + let (new_offsets, new_indices, _) = transpose_cs( + self.major_dim(), + self.minor_dim(), + self.major_offsets(), + self.minor_indices(), + &values); + // TODO: Skip checks + Self::try_from_offsets_and_indices(self.minor_dim(), + self.major_dim(), + new_offsets, + new_indices) + .expect("Internal error: Transpose should never fail.") + } } /// Error type for `SparsityPattern` format errors.