Merge pull request #1165 from geo-ant/feature/parallel-column-iterators

Parallel Column Iterators with Rayon

.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:

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"]

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<T, R, C>> ExactSizeIte
 }
 
 /*
- *
  * Column iterators.
  *
  */
@@ -296,12 +301,33 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut<T, R, C>> ExactSizeIte
 /// An iterator through the columns of a matrix.
 pub struct ColumnIter<'a, T, R: Dim, C: Dim, S: RawStorage<T, R, C>> {
     mat: &'a Matrix<T, R, C, S>,
-    curr: usize,
+    range: Range<usize>,
 }
 
 impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage<T, R, C>> ColumnIter<'a, T, R, C, S> {
+    /// a new column iterator covering all columns of the matrix
     pub(crate) fn new(mat: &'a Matrix<T, R, C, S>) -> 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<T, R, C>> Iterator for ColumnIter
 
     #[inline]
     fn next(&mut self) -> Option<Self::Item> {
-        if self.curr < self.mat.ncols() {
+        debug_assert!(self.range.start <= self.range.end);
-            let res = self.mat.column(self.curr);
+        if self.range.start < self.range.end {
-            self.curr += 1;
+            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<T, R, C>> Iterator for ColumnIter
 
     #[inline]
     fn size_hint(&self) -> (usize, Option<usize>) {
-        (
+        let hint = self.range.len();
-            self.mat.ncols() - self.curr,
+        (hint, Some(hint))
-            Some(self.mat.ncols() - self.curr),
-        )
     }
 
     #[inline]
     fn count(self) -> usize {
-        self.mat.ncols() - self.curr
+        self.range.len()
+    }
+}
+
+impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage<T, R, C>> DoubleEndedIterator
+    for ColumnIter<'a, T, R, C, S>
+{
+    fn next_back(&mut self) -> Option<Self::Item> {
+        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<T, R, C>> 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<T, R, C>> ExactSizeIterat
 #[derive(Debug)]
 pub struct ColumnIterMut<'a, T, R: Dim, C: Dim, S: RawStorageMut<T, R, C>> {
     mat: *mut Matrix<T, R, C, S>,
-    curr: usize,
+    range: Range<usize>,
     phantom: PhantomData<&'a mut Matrix<T, R, C, S>>,
 }
 
 impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut<T, R, C>> ColumnIterMut<'a, T, R, C, S> {
     pub(crate) fn new(mat: &'a mut Matrix<T, R, C, S>) -> Self {
+        let range = 0..mat.ncols();
         ColumnIterMut {
             mat,
-            curr: 0,
+            range,
-            phantom: PhantomData,
+            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<T, R, C>> Iterator
     type Item = MatrixViewMut<'a, T, R, U1, S::RStride, S::CStride>;
 
     #[inline]
-    fn next(&mut self) -> Option<Self::Item> {
+    fn next(&'_ mut self) -> Option<Self::Item> {
-        if self.curr < self.ncols() {
+        debug_assert!(self.range.start <= self.range.end);
-            let res = unsafe { (*self.mat).column_mut(self.curr) };
+        if self.range.start < self.range.end {
-            self.curr += 1;
+            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<T, R, C>> Iterator
 
     #[inline]
     fn size_hint(&self) -> (usize, Option<usize>) {
-        (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<T, R, C>> ExactSizeIte
 {
     #[inline]
     fn len(&self) -> usize {
-        self.ncols() - self.curr
+        self.range.len()
+    }
+}
+
+impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut<T, R, C>> DoubleEndedIterator
+    for ColumnIterMut<'a, T, R, C, S>
+{
+    fn next_back(&mut self) -> Option<Self::Item> {
+        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
+        }
     }
 }

src/base/matrix.rs

@@ -101,6 +101,7 @@ pub type MatrixCross<T, R1, C1, R2, C2> =
 ///
 /// #### Iteration, map, and fold
 /// - [Iteration on components, rows, and columns <span style="float:right;">`iter`, `column_iter`…</span>](#iteration-on-components-rows-and-columns)
+/// - [Parallel iterators using rayon <span style="float:right;">`par_column_iter`, `par_column_iter_mut`…</span>](#parallel-iterators-using-rayon)
 /// - [Elementwise mapping and folding <span style="float:right;">`map`, `fold`, `zip_map`…</span>](#elementwise-mapping-and-folding)
 /// - [Folding or columns and rows <span style="float:right;">`compress_rows`, `compress_columns`…</span>](#folding-on-columns-and-rows)
 ///

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;
 

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<T, R, Cols>> {
+    mat: &'a Matrix<T, R, Cols, S>,
+}
+
+impl<'a, T, R: Dim, Cols: Dim, S: RawStorage<T, R, Cols>> ParColumnIter<'a, T, R, Cols, S> {
+    /// Create a new parallel iterator for the given matrix.
+    fn new(matrix: &'a Matrix<T, R, Cols, S>) -> Self {
+        Self { mat: matrix }
+    }
+}
+
+#[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))]
+impl<'a, T, R: Dim, Cols: Dim, S: RawStorage<T, R, Cols>> 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<Consumer>(self, consumer: Consumer) -> Consumer::Result
+    where
+        Consumer: rayon::iter::plumbing::UnindexedConsumer<Self::Item>,
+    {
+        bridge(self, consumer)
+    }
+
+    fn opt_len(&self) -> Option<usize> {
+        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<T, R, Cols>> IndexedParallelIterator
+    for ParColumnIter<'a, T, R, Cols, S>
+where
+    T: Send + Sync + Scalar,
+    S: Sync,
+{
+    fn len(&self) -> usize {
+        self.mat.ncols()
+    }
+
+    fn drive<C: rayon::iter::plumbing::Consumer<Self::Item>>(self, consumer: C) -> C::Result {
+        bridge(self, consumer)
+    }
+
+    fn with_producer<CB: rayon::iter::plumbing::ProducerCallback<Self::Item>>(
+        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<T, R, Cols> + RawStorageMut<T, R, Cols>,
+> {
+    mat: &'a mut Matrix<T, R, Cols, S>,
+}
+
+#[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<T, R, Cols> + RawStorageMut<T, R, Cols>,
+{
+    /// create a new parallel iterator for the given matrix.
+    fn new(mat: &'a mut Matrix<T, R, Cols, S>) -> 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<T, R, Cols> + RawStorageMut<T, R, Cols>,
+    T: Send + Sync + Scalar,
+    S: Send + Sync,
+{
+    type Item = MatrixViewMut<'a, T, R, U1, S::RStride, S::CStride>;
+    fn drive_unindexed<C>(self, consumer: C) -> C::Result
+    where
+        C: rayon::iter::plumbing::UnindexedConsumer<Self::Item>,
+    {
+        bridge(self, consumer)
+    }
+
+    fn opt_len(&self) -> Option<usize> {
+        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<T, R, Cols> + RawStorageMut<T, R, Cols>,
+    T: Send + Sync + Scalar,
+    S: Send + Sync,
+{
+    fn drive<C: rayon::iter::plumbing::Consumer<Self::Item>>(self, consumer: C) -> C::Result {
+        bridge(self, consumer)
+    }
+
+    fn len(&self) -> usize {
+        self.mat.ncols()
+    }
+
+    fn with_producer<CB: rayon::iter::plumbing::ProducerCallback<Self::Item>>(
+        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<T, R: Dim, Cols: Dim, S: RawStorage<T, R, Cols>> Matrix<T, R, Cols, S>
+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<f64> = 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<f64> = 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<T, R, Cols>,
+    {
+        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<T, R, C>>(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<T, R, Cols>> 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<T, R, C>>(
+    ColumnIterMut<'a, T, R, C, S>,
+);
+
+impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut<T, R, C>> 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<T, R, C>> Send
+    for ColumnIterMut<'a, T, R, C, S>
+{
+}

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<f64, 2, 2> = 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<f64, 2, 2> = 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<f64> = 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::<f32>::zeros(400, 300);
+    let mut second = DMatrix::<f32>::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));
+}