diff --git a/.github/workflows/nalgebra-ci-build.yml b/.github/workflows/nalgebra-ci-build.yml index 9e500084..d7027127 100644 --- a/.github/workflows/nalgebra-ci-build.yml +++ b/.github/workflows/nalgebra-ci-build.yml @@ -61,7 +61,7 @@ jobs: steps: - uses: actions/checkout@v2 - name: test - run: cargo test --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,rkyv-safe-deser; + run: cargo test --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,rkyv-safe-deser,rayon; test-nalgebra-glm: runs-on: ubuntu-latest steps: diff --git a/Cargo.toml b/Cargo.toml index 9d0b1a7c..82528ed4 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -34,6 +34,7 @@ libm-force = [ "simba/libm_force" ] macros = [ "nalgebra-macros" ] cuda = [ "cust_core", "simba/cuda" ] + # Conversion convert-mint = [ "mint" ] convert-bytemuck = [ "bytemuck" ] @@ -101,7 +102,7 @@ glam020 = { package = "glam", version = "0.20", optional = true } glam021 = { package = "glam", version = "0.21", optional = true } glam022 = { package = "glam", version = "0.22", optional = true } cust_core = { version = "0.1", optional = true } - +rayon = { version = "1.6", optional = true } [dev-dependencies] serde_json = "1.0" @@ -137,3 +138,5 @@ lto = true [package.metadata.docs.rs] # Enable all the features when building the docs on docs.rs all-features = true +# define the configuration attribute `docsrs` +rustdoc-args = ["--cfg", "docsrs"] diff --git a/src/base/iter.rs b/src/base/iter.rs index 991c24dc..0e4aa8d4 100644 --- a/src/base/iter.rs +++ b/src/base/iter.rs @@ -1,5 +1,11 @@ //! Matrix iterators. +// only enables the `doc_cfg` feature when +// the `docsrs` configuration attribute is defined +#![cfg_attr(docsrs, feature(doc_cfg))] + +use core::fmt::Debug; +use core::ops::Range; use std::iter::FusedIterator; use std::marker::PhantomData; use std::mem; @@ -288,7 +294,6 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut> ExactSizeIte } /* - * * Column iterators. * */ @@ -296,12 +301,33 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut> ExactSizeIte /// An iterator through the columns of a matrix. pub struct ColumnIter<'a, T, R: Dim, C: Dim, S: RawStorage> { mat: &'a Matrix, - curr: usize, + range: Range, } impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> ColumnIter<'a, T, R, C, S> { + /// a new column iterator covering all columns of the matrix pub(crate) fn new(mat: &'a Matrix) -> Self { - ColumnIter { mat, curr: 0 } + ColumnIter { + mat, + range: 0..mat.ncols(), + } + } + + pub(crate) fn split_at(self, index: usize) -> (Self, Self) { + // SAFETY: this makes sur the generated ranges are valid. + let split_pos = (self.range.start + index).min(self.range.end); + + let left_iter = ColumnIter { + mat: self.mat, + range: self.range.start..split_pos, + }; + + let right_iter = ColumnIter { + mat: self.mat, + range: split_pos..self.range.end, + }; + + (left_iter, right_iter) } } @@ -310,9 +336,10 @@ impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> Iterator for ColumnIter #[inline] fn next(&mut self) -> Option { - if self.curr < self.mat.ncols() { - let res = self.mat.column(self.curr); - self.curr += 1; + debug_assert!(self.range.start <= self.range.end); + if self.range.start < self.range.end { + let res = self.mat.column(self.range.start); + self.range.start += 1; Some(res) } else { None @@ -321,15 +348,29 @@ impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> Iterator for ColumnIter #[inline] fn size_hint(&self) -> (usize, Option) { - ( - self.mat.ncols() - self.curr, - Some(self.mat.ncols() - self.curr), - ) + let hint = self.range.len(); + (hint, Some(hint)) } #[inline] fn count(self) -> usize { - self.mat.ncols() - self.curr + self.range.len() + } +} + +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> DoubleEndedIterator + for ColumnIter<'a, T, R, C, S> +{ + fn next_back(&mut self) -> Option { + debug_assert!(self.range.start <= self.range.end); + if !self.range.is_empty() { + self.range.end -= 1; + debug_assert!(self.range.end < self.mat.ncols()); + debug_assert!(self.range.end >= self.range.start); + Some(self.mat.column(self.range.end)) + } else { + None + } } } @@ -338,7 +379,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorage> ExactSizeIterat { #[inline] fn len(&self) -> usize { - self.mat.ncols() - self.curr + self.range.end - self.range.start } } @@ -346,19 +387,39 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorage> ExactSizeIterat #[derive(Debug)] pub struct ColumnIterMut<'a, T, R: Dim, C: Dim, S: RawStorageMut> { mat: *mut Matrix, - curr: usize, + range: Range, phantom: PhantomData<&'a mut Matrix>, } impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> ColumnIterMut<'a, T, R, C, S> { pub(crate) fn new(mat: &'a mut Matrix) -> Self { + let range = 0..mat.ncols(); ColumnIterMut { mat, - curr: 0, - phantom: PhantomData, + range, + phantom: Default::default(), } } + pub(crate) fn split_at(self, index: usize) -> (Self, Self) { + // SAFETY: this makes sur the generated ranges are valid. + let split_pos = (self.range.start + index).min(self.range.end); + + let left_iter = ColumnIterMut { + mat: self.mat, + range: self.range.start..split_pos, + phantom: Default::default(), + }; + + let right_iter = ColumnIterMut { + mat: self.mat, + range: split_pos..self.range.end, + phantom: Default::default(), + }; + + (left_iter, right_iter) + } + fn ncols(&self) -> usize { unsafe { (*self.mat).ncols() } } @@ -370,10 +431,11 @@ impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> Iterator type Item = MatrixViewMut<'a, T, R, U1, S::RStride, S::CStride>; #[inline] - fn next(&mut self) -> Option { - if self.curr < self.ncols() { - let res = unsafe { (*self.mat).column_mut(self.curr) }; - self.curr += 1; + fn next(&'_ mut self) -> Option { + debug_assert!(self.range.start <= self.range.end); + if self.range.start < self.range.end { + let res = unsafe { (*self.mat).column_mut(self.range.start) }; + self.range.start += 1; Some(res) } else { None @@ -382,12 +444,13 @@ impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> Iterator #[inline] fn size_hint(&self) -> (usize, Option) { - (self.ncols() - self.curr, Some(self.ncols() - self.curr)) + let hint = self.range.len(); + (hint, Some(hint)) } #[inline] fn count(self) -> usize { - self.ncols() - self.curr + self.range.len() } } @@ -396,6 +459,22 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut> ExactSizeIte { #[inline] fn len(&self) -> usize { - self.ncols() - self.curr + self.range.len() + } +} + +impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut> DoubleEndedIterator + for ColumnIterMut<'a, T, R, C, S> +{ + fn next_back(&mut self) -> Option { + debug_assert!(self.range.start <= self.range.end); + if !self.range.is_empty() { + self.range.end -= 1; + debug_assert!(self.range.end < self.ncols()); + debug_assert!(self.range.end >= self.range.start); + Some(unsafe { (*self.mat).column_mut(self.range.end) }) + } else { + None + } } } diff --git a/src/base/matrix.rs b/src/base/matrix.rs index 700e2f02..f0fa8d81 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -101,6 +101,7 @@ pub type MatrixCross = /// /// #### Iteration, map, and fold /// - [Iteration on components, rows, and columns `iter`, `column_iter`…](#iteration-on-components-rows-and-columns) +/// - [Parallel iterators using rayon `par_column_iter`, `par_column_iter_mut`…](#parallel-iterators-using-rayon) /// - [Elementwise mapping and folding `map`, `fold`, `zip_map`…](#elementwise-mapping-and-folding) /// - [Folding or columns and rows `compress_rows`, `compress_columns`…](#folding-on-columns-and-rows) /// diff --git a/src/base/mod.rs b/src/base/mod.rs index 4cbcff93..0f09cc33 100644 --- a/src/base/mod.rs +++ b/src/base/mod.rs @@ -42,6 +42,9 @@ mod min_max; /// Mechanisms for working with values that may not be initialized. pub mod uninit; +#[cfg(feature = "rayon")] +pub mod par_iter; + #[cfg(feature = "rkyv-serialize-no-std")] mod rkyv_wrappers; diff --git a/src/base/par_iter.rs b/src/base/par_iter.rs new file mode 100644 index 00000000..af5e1cb7 --- /dev/null +++ b/src/base/par_iter.rs @@ -0,0 +1,285 @@ +//! Parallel iterators for matrices compatible with rayon. + +// only enables the `doc_cfg` feature when +// the `docsrs` configuration attribute is defined +#![cfg_attr(docsrs, feature(doc_cfg))] + +use crate::{ + iter::{ColumnIter, ColumnIterMut}, + Dim, Matrix, MatrixView, MatrixViewMut, RawStorage, RawStorageMut, Scalar, U1, +}; +use rayon::iter::plumbing::Producer; +use rayon::{iter::plumbing::bridge, prelude::*}; + +/// A rayon parallel iterator over the colums of a matrix. It is created +/// using the [`par_column_iter`] method of [`Matrix`]. +/// +/// *Only available if compiled with the feature `rayon`.* +/// [`par_column_iter`]: crate::Matrix::par_column_iter +/// [`Matrix`]: crate::Matrix +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +pub struct ParColumnIter<'a, T, R: Dim, Cols: Dim, S: RawStorage> { + mat: &'a Matrix, +} + +impl<'a, T, R: Dim, Cols: Dim, S: RawStorage> ParColumnIter<'a, T, R, Cols, S> { + /// Create a new parallel iterator for the given matrix. + fn new(matrix: &'a Matrix) -> Self { + Self { mat: matrix } + } +} + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +impl<'a, T, R: Dim, Cols: Dim, S: RawStorage> ParallelIterator + for ParColumnIter<'a, T, R, Cols, S> +where + T: Sync + Send + Scalar, + S: Sync, +{ + type Item = MatrixView<'a, T, R, U1, S::RStride, S::CStride>; + + fn drive_unindexed(self, consumer: Consumer) -> Consumer::Result + where + Consumer: rayon::iter::plumbing::UnindexedConsumer, + { + bridge(self, consumer) + } + + fn opt_len(&self) -> Option { + Some(self.mat.ncols()) + } +} + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +/// *Only available if compiled with the feature `rayon`.* +impl<'a, T, R: Dim, Cols: Dim, S: RawStorage> IndexedParallelIterator + for ParColumnIter<'a, T, R, Cols, S> +where + T: Send + Sync + Scalar, + S: Sync, +{ + fn len(&self) -> usize { + self.mat.ncols() + } + + fn drive>(self, consumer: C) -> C::Result { + bridge(self, consumer) + } + + fn with_producer>( + self, + callback: CB, + ) -> CB::Output { + let producer = ColumnProducer(ColumnIter::new(self.mat)); + callback.callback(producer) + } +} + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +/// A rayon parallel iterator through the mutable columns of a matrix. +/// *Only available if compiled with the feature `rayon`.* +pub struct ParColumnIterMut< + 'a, + T, + R: Dim, + Cols: Dim, + S: RawStorage + RawStorageMut, +> { + mat: &'a mut Matrix, +} + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +/// *only availabe if compiled with the feature `rayon`* +impl<'a, T, R, Cols, S> ParColumnIterMut<'a, T, R, Cols, S> +where + R: Dim, + Cols: Dim, + S: RawStorage + RawStorageMut, +{ + /// create a new parallel iterator for the given matrix. + fn new(mat: &'a mut Matrix) -> Self { + Self { mat } + } +} + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +/// *Only available if compiled with the feature `rayon`* +impl<'a, T, R, Cols, S> ParallelIterator for ParColumnIterMut<'a, T, R, Cols, S> +where + R: Dim, + Cols: Dim, + S: RawStorage + RawStorageMut, + T: Send + Sync + Scalar, + S: Send + Sync, +{ + type Item = MatrixViewMut<'a, T, R, U1, S::RStride, S::CStride>; + fn drive_unindexed(self, consumer: C) -> C::Result + where + C: rayon::iter::plumbing::UnindexedConsumer, + { + bridge(self, consumer) + } + + fn opt_len(&self) -> Option { + Some(self.mat.ncols()) + } +} + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +/// *Only available if compiled with the feature `rayon`* +impl<'a, T, R, Cols, S> IndexedParallelIterator for ParColumnIterMut<'a, T, R, Cols, S> +where + R: Dim, + Cols: Dim, + S: RawStorage + RawStorageMut, + T: Send + Sync + Scalar, + S: Send + Sync, +{ + fn drive>(self, consumer: C) -> C::Result { + bridge(self, consumer) + } + + fn len(&self) -> usize { + self.mat.ncols() + } + + fn with_producer>( + self, + callback: CB, + ) -> CB::Output { + let producer = ColumnProducerMut(ColumnIterMut::new(self.mat)); + callback.callback(producer) + } +} + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +/// # Parallel iterators using `rayon` +/// *Only available if compiled with the feature `rayon`* +impl> Matrix +where + T: Send + Sync + Scalar, + S: Sync, +{ + /// Iterate through the columns of the matrix in parallel using rayon. + /// This iterates over *immutable* references ot the columns of the matrix, + /// if *mutable* access to the columns is required, use [`par_column_iter_mut`] + /// instead. + /// + /// # Example + /// Using parallel column iterators to calculate the sum of the maximum + /// elements in each column: + /// ``` + /// use nalgebra::{dmatrix, DMatrix}; + /// use rayon::prelude::*; + /// + /// let matrix : DMatrix = dmatrix![1.0, 0.0, 5.0; + /// 2.0, 4.0, 1.0; + /// 3.0, 2.0, 2.0; + /// ]; + /// let sum_of_max :f64 = matrix + /// .par_column_iter() + /// .map(|col| col.max()) + /// .sum(); + /// + /// assert_eq!(sum_of_max,3.0 + 4.0 + 5.0); + /// + /// ``` + /// + /// [`par_column_iter_mut`]: crate::Matrix::par_column_iter_mut + pub fn par_column_iter(&self) -> ParColumnIter<'_, T, R, Cols, S> { + ParColumnIter::new(self) + } + + /// Mutably iterate through the columns of this matrix in parallel using rayon. + /// Allows mutable access to the columns in parallel using mutable references. + /// If mutable access to the columns is not required rather use [`par_column_iter`] + /// instead. + /// + /// # Example + /// Normalize each column of a matrix with respect to its own maximum value. + /// + /// ``` + /// use nalgebra::{dmatrix, DMatrix}; + /// use rayon::prelude::*; + /// + /// let mut matrix : DMatrix = dmatrix![ + /// 2.0, 4.0, 6.0; + /// 1.0, 2.0, 3.0; + /// ]; + /// matrix.par_column_iter_mut().for_each(|mut col| col /= col.max()); + /// + /// assert_eq!(matrix, dmatrix![1.0, 1.0, 1.0; 0.5, 0.5, 0.5]); + /// ``` + /// + /// [`par_column_iter`]: crate::Matrix::par_column_iter + pub fn par_column_iter_mut(&mut self) -> ParColumnIterMut<'_, T, R, Cols, S> + where + S: RawStorageMut, + { + ParColumnIterMut::new(self) + } +} + +/// A private helper newtype that wraps the `ColumnIter` and implements +/// the rayon `Producer` trait. It's just here so we don't have to make the +/// rayon trait part of the public interface of the `ColumnIter`. +struct ColumnProducer<'a, T, R: Dim, C: Dim, S: RawStorage>(ColumnIter<'a, T, R, C, S>); + +#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] +/// *only available if compiled with the feature `rayon`* +impl<'a, T, R: Dim, Cols: Dim, S: RawStorage> Producer + for ColumnProducer<'a, T, R, Cols, S> +where + T: Send + Sync + Scalar, + S: Sync, +{ + type Item = MatrixView<'a, T, R, U1, S::RStride, S::CStride>; + type IntoIter = ColumnIter<'a, T, R, Cols, S>; + + #[inline] + fn into_iter(self) -> Self::IntoIter { + self.0 + } + + #[inline] + fn split_at(self, index: usize) -> (Self, Self) { + // The index is relative to the size of this current iterator. + // It will always start at zero so it serves as an offset. + let (left_iter, right_iter) = self.0.split_at(index); + (Self(left_iter), Self(right_iter)) + } +} + +/// See `ColumnProducer`. A private wrapper newtype that keeps the Producer +/// implementation private +struct ColumnProducerMut<'a, T, R: Dim, C: Dim, S: RawStorageMut>( + ColumnIterMut<'a, T, R, C, S>, +); + +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> Producer + for ColumnProducerMut<'a, T, R, C, S> +where + T: Send + Sync + Scalar, + S: Send + Sync, +{ + type Item = MatrixViewMut<'a, T, R, U1, S::RStride, S::CStride>; + type IntoIter = ColumnIterMut<'a, T, R, C, S>; + + fn into_iter(self) -> Self::IntoIter { + self.0 + } + + fn split_at(self, index: usize) -> (Self, Self) { + // The index is relative to the size of this current iterator + // it will always start at zero so it serves as an offset. + let (left_iter, right_iter) = self.0.split_at(index); + (Self(left_iter), Self(right_iter)) + } +} + +/// this implementation is safe because we are enforcing exclusive access +/// to the columns through the active range of the iterator +unsafe impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut> Send + for ColumnIterMut<'a, T, R, C, S> +{ +} diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs index 219845d4..27926a27 100644 --- a/tests/core/matrix.rs +++ b/tests/core/matrix.rs @@ -1136,3 +1136,153 @@ fn omatrix_to_string() { (svec.to_string(), smatr.to_string()) ); } + +#[test] +fn column_iteration() { + // dynamic matrix + let dmat = nalgebra::dmatrix![ + 13,14,15; + 23,24,25; + 33,34,35; + ]; + let mut col_iter = dmat.column_iter(); + assert_eq!(col_iter.next(), Some(dmat.column(0))); + assert_eq!(col_iter.next(), Some(dmat.column(1))); + assert_eq!(col_iter.next(), Some(dmat.column(2))); + assert_eq!(col_iter.next(), None); + + // statically sized matrix + let smat: nalgebra::SMatrix = nalgebra::matrix![1.0, 2.0; 3.0, 4.0]; + let mut col_iter = smat.column_iter(); + assert_eq!(col_iter.next(), Some(smat.column(0))); + assert_eq!(col_iter.next(), Some(smat.column(1))); + assert_eq!(col_iter.next(), None); +} + +#[test] +fn column_iteration_mut() { + let mut dmat = nalgebra::dmatrix![ + 13,14,15; + 23,24,25; + 33,34,35; + ]; + let mut cloned = dmat.clone(); + let mut col_iter = dmat.column_iter_mut(); + assert_eq!(col_iter.next(), Some(cloned.column_mut(0))); + assert_eq!(col_iter.next(), Some(cloned.column_mut(1))); + assert_eq!(col_iter.next(), Some(cloned.column_mut(2))); + assert_eq!(col_iter.next(), None); + + // statically sized matrix + let mut smat: nalgebra::SMatrix = nalgebra::matrix![1.0, 2.0; 3.0, 4.0]; + let mut cloned = smat.clone(); + let mut col_iter = smat.column_iter_mut(); + assert_eq!(col_iter.next(), Some(cloned.column_mut(0))); + assert_eq!(col_iter.next(), Some(cloned.column_mut(1))); + assert_eq!(col_iter.next(), None); +} + +#[test] +fn column_iteration_double_ended() { + let dmat = nalgebra::dmatrix![ + 13,14,15,16,17; + 23,24,25,26,27; + 33,34,35,36,37; + ]; + let mut col_iter = dmat.column_iter(); + assert_eq!(col_iter.next(), Some(dmat.column(0))); + assert_eq!(col_iter.next(), Some(dmat.column(1))); + assert_eq!(col_iter.next_back(), Some(dmat.column(4))); + assert_eq!(col_iter.next_back(), Some(dmat.column(3))); + assert_eq!(col_iter.next(), Some(dmat.column(2))); + assert_eq!(col_iter.next_back(), None); + assert_eq!(col_iter.next(), None); +} + +#[test] +fn column_iterator_double_ended_mut() { + let mut dmat = nalgebra::dmatrix![ + 13,14,15,16,17; + 23,24,25,26,27; + 33,34,35,36,37; + ]; + let mut cloned = dmat.clone(); + let mut col_iter_mut = dmat.column_iter_mut(); + assert_eq!(col_iter_mut.next(), Some(cloned.column_mut(0))); + assert_eq!(col_iter_mut.next(), Some(cloned.column_mut(1))); + assert_eq!(col_iter_mut.next_back(), Some(cloned.column_mut(4))); + assert_eq!(col_iter_mut.next_back(), Some(cloned.column_mut(3))); + assert_eq!(col_iter_mut.next(), Some(cloned.column_mut(2))); + assert_eq!(col_iter_mut.next_back(), None); + assert_eq!(col_iter_mut.next(), None); +} + +#[test] +#[cfg(feature = "rayon")] +fn parallel_column_iteration() { + use nalgebra::dmatrix; + use rayon::prelude::*; + let dmat: DMatrix = dmatrix![ + 13.,14.; + 23.,24.; + 33.,34.; + ]; + let cloned = dmat.clone(); + // test that correct columns are iterated over + dmat.par_column_iter().enumerate().for_each(|(idx, col)| { + assert_eq!(col, cloned.column(idx)); + }); + // test that a more complex expression produces the same + // result as the serial equivalent + let par_result: f64 = dmat.par_column_iter().map(|col| col.norm()).sum(); + let ser_result: f64 = dmat.column_iter().map(|col| col.norm()).sum(); + assert_eq!(par_result, ser_result); + + // repeat this test using mutable iterators + let mut dmat = dmat; + dmat.par_column_iter_mut() + .enumerate() + .for_each(|(idx, col)| { + assert_eq!(col, cloned.column(idx)); + }); + + let par_mut_result: f64 = dmat.par_column_iter_mut().map(|col| col.norm()).sum(); + assert_eq!(par_mut_result, ser_result); +} + +#[test] +#[cfg(feature = "rayon")] +fn column_iteration_mut_double_ended() { + let dmat = nalgebra::dmatrix![ + 13,14,15,16,17; + 23,24,25,26,27; + 33,34,35,36,37; + ]; + let cloned = dmat.clone(); + let mut col_iter = dmat.column_iter(); + assert_eq!(col_iter.next(), Some(cloned.column(0))); + assert_eq!(col_iter.next(), Some(cloned.column(1))); + assert_eq!(col_iter.next_back(), Some(cloned.column(4))); + assert_eq!(col_iter.next_back(), Some(cloned.column(3))); + assert_eq!(col_iter.next(), Some(cloned.column(2))); + assert_eq!(col_iter.next_back(), None); + assert_eq!(col_iter.next(), None); +} + +#[test] +#[cfg(feature = "rayon")] +fn parallel_column_iteration_mut() { + use rayon::prelude::*; + let mut first = DMatrix::::zeros(400, 300); + let mut second = DMatrix::::zeros(400, 300); + first + .column_iter_mut() + .enumerate() + .for_each(|(idx, mut col)| col[idx] = 1.); + second + .par_column_iter_mut() + .enumerate() + .for_each(|(idx, mut col)| col[idx] = 1.); + assert_eq!(first, second); + assert_eq!(second, DMatrix::identity(400, 300)); +}