Merge branch 'dimforge:dev' into dev

2022-09-14 18:12:42 +02:00 · 2022-09-14 18:12:42 +02:00 · 5a87884b50
parent 743eaa4747 7adecdef6c
commit 5a87884b50
28 changed files with 931 additions and 139 deletions
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@ -4,6 +4,21 @@ documented here.

 This project adheres to [Semantic Versioning](https://semver.org/).

+## [0.31.1] (31 July 2022)
+
+### Modified
+- Improve performances of multiplication of two sparse matrices.
+
+### Added
+- Add `Matrix::from_row_iterator` to build a matrix from an iterator yielding components in row-major order.
+- Add support for conversion from/to types of `glam` 0.21.
+- `nalgebra-sparse`: add support for the matrix-market export of sparse matrices.
+- `nalgebra-lapack`: add a `GE` for solving the generalized eigenvalues problem.
+
+### Fixed
+- Fix `Rotation3::from_matrix` and `UnitQuaternion::from_matrix` when the input matrix is already a valid
+  rotation matrix.
+
 ## [0.31.0] (30 Apr. 2022)

 ### Breaking changes
--- a/Cargo.toml
+++ b/Cargo.toml
@ -1,6 +1,6 @@
 [package]
 name    = "nalgebra"
-version = "0.31.0"
+version = "0.31.1"
 authors = [ "Sébastien Crozet <developer@crozet.re>" ]

 description = "General-purpose linear algebra library with transformations and statically-sized or dynamically-sized matrices."
@ -44,6 +44,7 @@ convert-glam017 = [ "glam017" ]
 convert-glam018 = [ "glam018" ]
 convert-glam019 = [ "glam019" ]
 convert-glam020 = [ "glam020" ]
+convert-glam021 = [ "glam021" ]

 # Serialization
 ## To use serde in a #[no-std] environment, enable the
@ -95,6 +96,7 @@ glam017        = { package = "glam", version = "0.17", optional = true }
 glam018        = { package = "glam", version = "0.18", optional = true }
 glam019        = { package = "glam", version = "0.19", optional = true }
 glam020        = { package = "glam", version = "0.20", optional = true }
+glam021        = { package = "glam", version = "0.21", optional = true }
 cust_core      = { version = "0.1", optional = true }


@ -130,6 +132,6 @@ required-features = ["rand"]
 lto = true

 [package.metadata.docs.rs]
-# Enable certain features when building docs for docs.rs
-features = [ "proptest-support", "compare", "macros", "rand" ]
+# Enable all the features when building the docs on docs.rs
+all-features = true

--- a/nalgebra-sparse/Cargo.toml
+++ b/nalgebra-sparse/Cargo.toml
@ -1,6 +1,6 @@
 [package]
 name = "nalgebra-sparse"
-version = "0.7.0"
+version = "0.7.1"
 authors = [ "Andreas Longva", "Sébastien Crozet <developer@crozet.re>" ]
 edition = "2018"
 description = "Sparse matrix computation based on nalgebra."
@ -36,8 +36,9 @@ serde = { version = "1.0", default-features = false, features = [ "derive" ], op
 itertools = "0.10"
 matrixcompare = { version = "0.3.0", features = [ "proptest-support" ] }
 nalgebra = { version="0.31", path = "../", features = ["compare"] }
+tempfile = "3.3"
 serde_json = "1.0"

 [package.metadata.docs.rs]
 # Enable certain features when building docs for docs.rs
-features = [ "proptest-support", "compare" ]
+features = [ "proptest-support", "compare", "io"]
--- a/nalgebra-sparse/src/coo.rs
+++ b/nalgebra-sparse/src/coo.rs
@ -211,6 +211,13 @@ impl<T> CooMatrix<T> {
        self.values.push(v);
    }

+    /// Clear all triplets from the matrix.
+    pub fn clear_triplets(&mut self) {
+        self.col_indices.clear();
+        self.row_indices.clear();
+        self.values.clear();
+    }
+
    /// The number of rows in the matrix.
    #[inline]
    #[must_use]
--- a/nalgebra-sparse/src/io/matrix_market.rs
+++ b/nalgebra-sparse/src/io/matrix_market.rs
@ -1,9 +1,9 @@
 //! Implementation of matrix market io code.
 //!
 //! See the [website](https://math.nist.gov/MatrixMarket/formats.html) or the [paper](https://www.researchgate.net/publication/2630533_The_Matrix_Market_Exchange_Formats_Initial_Design) for more details about matrix market.
-use crate::coo::CooMatrix;
 use crate::SparseFormatError;
 use crate::SparseFormatErrorKind;
+use crate::{CooMatrix, CscMatrix, CsrMatrix};
 use nalgebra::Complex;
 use pest::iterators::Pairs;
 use pest::Parser;
@ -12,7 +12,8 @@ use std::convert::Infallible;
 use std::convert::TryFrom;
 use std::fmt;
 use std::fmt::Formatter;
-use std::fs;
+use std::fs::{self, File};
+use std::io::{BufWriter, Write};
 use std::num::ParseIntError;
 use std::num::TryFromIntError;
 use std::path::Path;
@ -267,7 +268,7 @@ impl fmt::Display for MatrixMarketError {
                write!(f, "InvalidHeader,")?;
            }
            MatrixMarketErrorKind::EntryMismatch => {
-                write!(f, "EntryNumUnmatched,")?;
+                write!(f, "EntryMismatch,")?;
            }
            MatrixMarketErrorKind::TypeMismatch => {
                write!(f, "TypeMismatch,")?;
@ -288,7 +289,7 @@ impl fmt::Display for MatrixMarketError {
                write!(f, "NotLowerTriangle,")?;
            }
            MatrixMarketErrorKind::NonSquare => {
-                write!(f, "NotSquareMatrix,")?;
+                write!(f, "NonSquare,")?;
            }
        }
        write!(f, " message: {}", self.message)
@ -506,6 +507,21 @@ mod internal {
        fn negative(self) -> Result<Self, MatrixMarketError>;
        /// When matrix is a Hermitian matrix, it will convert itself to its conjugate.
        fn conjugate(self) -> Result<Self, MatrixMarketError>;
+        /// Returns the name of SupportedMatrixMarketScalar, used when write the matrix
+        fn typename() -> &'static str;
+        /// Write the data self to w
+        fn write_matrix_market<W: std::fmt::Write>(&self, w: W) -> Result<(), std::fmt::Error>;
+    }
+
+    pub trait SupportedMatrixMarketExport<T: SupportedMatrixMarketScalar> {
+        /// iterate over triplets
+        fn triplet_iter(&self) -> Box<dyn Iterator<Item = (usize, usize, &T)> + '_>;
+        /// number of rows
+        fn nrows(&self) -> usize;
+        /// number of columns
+        fn ncols(&self) -> usize;
+        /// number of non-zeros
+        fn nnz(&self) -> usize;
    }
 }

@ -557,6 +573,17 @@ macro_rules! mm_int_impl {
            fn negative(self) -> Result<Self, MatrixMarketError> {
                Ok(-self)
            }
+            #[inline]
+            fn typename() -> &'static str {
+                "integer"
+            }
+            #[inline]
+            fn write_matrix_market<W: std::fmt::Write>(
+                &self,
+                mut w: W,
+            ) -> Result<(), std::fmt::Error> {
+                write!(w, "{}", self)
+            }
        }
    };
 }
@ -602,6 +629,17 @@ macro_rules! mm_real_impl {
            fn negative(self) -> Result<Self, MatrixMarketError> {
                Ok(-self)
            }
+            #[inline]
+            fn typename() -> &'static str {
+                "real"
+            }
+            #[inline]
+            fn write_matrix_market<W: std::fmt::Write>(
+                &self,
+                mut w: W,
+            ) -> Result<(), std::fmt::Error> {
+                write!(w, "{}", self)
+            }
        }
    };
 }
@ -648,6 +686,17 @@ macro_rules! mm_complex_impl {
            fn negative(self) -> Result<Self, MatrixMarketError> {
                Ok(-self)
            }
+            #[inline]
+            fn typename() -> &'static str {
+                "complex"
+            }
+            #[inline]
+            fn write_matrix_market<W: std::fmt::Write>(
+                &self,
+                mut w: W,
+            ) -> Result<(), std::fmt::Error> {
+                write!(w, "{} {}", self.re, self.im)
+            }
        }
    };
 }
@ -697,6 +746,17 @@ macro_rules! mm_pattern_impl {
                    format!("Pattern type has no negative"),
                ))
            }
+            #[inline]
+            fn typename() -> &'static str {
+                "pattern"
+            }
+            #[inline]
+            fn write_matrix_market<W: std::fmt::Write>(
+                &self,
+                mut _w: W,
+            ) -> Result<(), std::fmt::Error> {
+                Ok(())
+            }
        }
    };
 }
@ -715,6 +775,46 @@ mm_complex_impl!(f64);

 mm_pattern_impl!(());

+/// A marker trait for sparse matrix types that can be exported to the matrix market format.
+///
+/// This is a sealed trait; it cannot be implemented by external crates. This is done in order to prevent leaking
+/// some of the implementation details we currently rely on. We may relax this restriction in the future.
+pub trait MatrixMarketExport<T: MatrixMarketScalar>:
+    internal::SupportedMatrixMarketExport<T>
+{
+}
+
+macro_rules! mm_matrix_impl {
+    ($T_MATRIX:ty) => {
+        impl<T: MatrixMarketScalar> MatrixMarketExport<T> for $T_MATRIX {}
+
+        impl<T: internal::SupportedMatrixMarketScalar> internal::SupportedMatrixMarketExport<T>
+            for $T_MATRIX
+        {
+            #[inline]
+            fn triplet_iter(&self) -> Box<dyn Iterator<Item = (usize, usize, &T)> + '_> {
+                Box::new(self.triplet_iter())
+            }
+            #[inline]
+            fn nrows(&self) -> usize {
+                self.nrows()
+            }
+            #[inline]
+            fn ncols(&self) -> usize {
+                self.ncols()
+            }
+            #[inline]
+            fn nnz(&self) -> usize {
+                self.nnz()
+            }
+        }
+    };
+}
+
+mm_matrix_impl!(CooMatrix<T>);
+mm_matrix_impl!(CsrMatrix<T>);
+mm_matrix_impl!(CscMatrix<T>);
+
 #[derive(Parser)]
 #[grammar = "io/matrix_market.pest"]
 struct MatrixMarketParser;
@ -1329,3 +1429,123 @@ fn next_dense_coordinate(
        }
    }
 }
+
+/// Save a sparse matrix as a Matrix Market format string.
+///
+/// The exporter only writes the matrix into `coordinate` and `general` format.
+///
+///
+/// Examples
+/// --------
+/// ```
+/// # use nalgebra_sparse::CooMatrix;
+/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+/// use nalgebra_sparse::io::{save_to_matrix_market_str};
+/// let expected_str = r#"%%matrixmarket matrix coordinate integer general
+/// % matrixmarket file generated by nalgebra-sparse.
+/// 5 4 2
+/// 1 1 10
+/// 2 3 5
+/// "#;
+/// let row_indices = vec![0,1];
+/// let col_indices = vec![0,2];
+/// let values = vec![10,5];
+/// let matrix = CooMatrix::try_from_triplets(5,4,row_indices,col_indices,values)?;
+/// let generated_matrixmarket_str = save_to_matrix_market_str(&matrix);
+/// assert_eq!(expected_str,generated_matrixmarket_str);
+/// # Ok(()) }
+/// ```
+pub fn save_to_matrix_market_str<T, S>(sparse_matrix: &S) -> String
+where
+    T: MatrixMarketScalar,
+    S: MatrixMarketExport<T>,
+{
+    let mut bytes = Vec::<u8>::new();
+    // This will call impl<A: Allocator> Write for Vec<u8, A>
+    // The vector will grow as needed.
+    // So, unwrap here won't cause any issue.
+    save_to_matrix_market(&mut bytes, sparse_matrix).unwrap();
+
+    String::from_utf8(bytes)
+        .expect("Unexpected non UTF-8 data was generated when export to matrix market string")
+}
+
+/// Save a sparse matrix to a Matrix Market format file.
+///
+/// The exporter only saves the matrix with the `coordinate` and `general` matrix market formats.
+///
+/// Errors
+/// --------
+///
+/// See [MatrixMarketErrorKind] for a list of possible error conditions.
+///
+/// Examples
+/// --------
+/// ```no_run
+/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+/// use nalgebra_sparse::io::{save_to_matrix_market_file,load_coo_from_matrix_market_str};
+/// let str = r#"
+/// %%matrixmarket matrix coordinate integer general
+/// 5 4 2
+/// 1 1 10
+/// 2 3 5
+/// "#;
+/// let matrix = load_coo_from_matrix_market_str::<i32>(&str)?;
+/// save_to_matrix_market_file(&matrix,"path/to/matrix.mtx")?;
+/// # Ok(()) }
+/// ```
+pub fn save_to_matrix_market_file<T, S, P>(sparse_matrix: &S, path: P) -> Result<(), std::io::Error>
+where
+    T: MatrixMarketScalar,
+    S: MatrixMarketExport<T>,
+    P: AsRef<Path>,
+{
+    let file = File::create(path)?;
+    let mut file = BufWriter::new(file);
+    save_to_matrix_market(&mut file, sparse_matrix)?;
+    // Quote from BufWriter doc.
+    // > It is critical to call flush before BufWriter<W> is dropped. Though dropping will attempt to flush the contents of the buffer, any errors that happen in the process of dropping will be ignored. Calling flush ensures that the buffer is empty and thus dropping will not even attempt file operations.
+    file.flush()
+        .expect("Unexpected error when flushing the buffer data to File");
+    Ok(())
+}
+
+/// Save a sparse matrix to an [std::io::Write] instance.
+///
+/// This is the most general save functionality. See [save_to_matrix_market_file] and
+/// [save_to_matrix_market_str] for higher-level functionality.
+pub fn save_to_matrix_market<T, S, W>(mut w: W, sparse_matrix: &S) -> Result<(), std::io::Error>
+where
+    T: MatrixMarketScalar,
+    S: MatrixMarketExport<T>,
+    W: Write,
+{
+    // write header
+    writeln!(
+        w,
+        "%%matrixmarket matrix coordinate {} general",
+        T::typename()
+    )?;
+
+    //write comment
+    writeln!(w, "% matrixmarket file generated by nalgebra-sparse.")?;
+
+    // write shape information
+    writeln!(
+        w,
+        "{} {} {}",
+        sparse_matrix.nrows(),
+        sparse_matrix.ncols(),
+        sparse_matrix.nnz()
+    )?;
+
+    //write triplets
+    let mut buffer = String::new();
+    for (r, c, d) in sparse_matrix.triplet_iter() {
+        buffer.clear();
+        d.write_matrix_market(&mut buffer)
+            .expect("Unexpected format error was generated when write to String");
+        writeln!(w, "{} {} {}", r + 1, c + 1, buffer)?;
+    }
+    Ok(())
+}
--- a/nalgebra-sparse/src/io/mod.rs
+++ b/nalgebra-sparse/src/io/mod.rs
@ -6,7 +6,7 @@
 //!
 //! | Format                                          |  Import    |   Export   |
 //! | ------------------------------------------------|------------|------------|
-//! | [Matrix market](#matrix-market-format)          |  Yes       |    No      |
+//! | [Matrix market](#matrix-market-format)          |  Yes       |    Yes     |
 //!
 //! [Matrix market]: https://math.nist.gov/MatrixMarket/formats.html
 //!
@ -19,10 +19,12 @@
 //! which also uses the Matrix Market file format.
 //!
 //! We currently offer functionality for importing a Matrix market file to an instance of a
-//! [CooMatrix](crate::CooMatrix) through the function [load_coo_from_matrix_market_file]. It is also possible to load
-//! a matrix stored in the matrix market format with the function [load_coo_from_matrix_market_str].
-//!
-//! Export is currently not implemented, but [planned](https://github.com/dimforge/nalgebra/issues/1037).
+//! [CooMatrix](crate::CooMatrix) through the function [load_coo_from_matrix_market_file],
+//! as well as functionality for writing various sparse matrices to the matrix market format
+//! through [save_to_matrix_market_file]. It is also possible to load
+//! a matrix stored as a string in the matrix market format with the function
+//! [load_coo_from_matrix_market_str], or similarly write to a string with
+//! [save_to_matrix_market_str].
 //!
 //! Our implementation is based on the [format description](https://math.nist.gov/MatrixMarket/formats.html)
 //! on the Matrix Market website and the
@ -32,7 +34,8 @@
 //! > "*The Matrix Market Exchange Formats: Initial Design.*" (1996).

 pub use self::matrix_market::{
-    load_coo_from_matrix_market_file, load_coo_from_matrix_market_str, MatrixMarketError,
-    MatrixMarketErrorKind, MatrixMarketScalar,
+    load_coo_from_matrix_market_file, load_coo_from_matrix_market_str, save_to_matrix_market,
+    save_to_matrix_market_file, save_to_matrix_market_str, MatrixMarketError,
+    MatrixMarketErrorKind, MatrixMarketExport, MatrixMarketScalar,
 };
 mod matrix_market;
--- a/nalgebra-sparse/src/lib.rs
+++ b/nalgebra-sparse/src/lib.rs
@ -143,8 +143,6 @@
 )]

 pub extern crate nalgebra as na;
-#[cfg(feature = "io")]
-extern crate pest;
 #[macro_use]
 #[cfg(feature = "io")]
 extern crate pest_derive;
--- a/nalgebra-sparse/src/ops/impl_std_ops.rs
+++ b/nalgebra-sparse/src/ops/impl_std_ops.rs
@ -3,7 +3,7 @@ use crate::csr::CsrMatrix;

 use crate::ops::serial::{
    spadd_csc_prealloc, spadd_csr_prealloc, spadd_pattern, spmm_csc_dense, spmm_csc_pattern,
-    spmm_csc_prealloc, spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc,
+    spmm_csc_prealloc_unchecked, spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc_unchecked,
 };
 use crate::ops::Op;
 use nalgebra::allocator::Allocator;
@ -112,9 +112,9 @@ macro_rules! impl_spmm {
    }
 }

-impl_spmm!(CsrMatrix, spmm_csr_pattern, spmm_csr_prealloc);
+impl_spmm!(CsrMatrix, spmm_csr_pattern, spmm_csr_prealloc_unchecked);
 // Need to switch order of operations for CSC pattern
-impl_spmm!(CscMatrix, spmm_csc_pattern, spmm_csc_prealloc);
+impl_spmm!(CscMatrix, spmm_csc_pattern, spmm_csc_prealloc_unchecked);

 /// Implements Scalar * Matrix operations for *concrete* scalar types. The reason this is necessary
 /// is that we are not able to implement Mul<Matrix<T>> for all T generically due to orphan rules.
--- a/nalgebra-sparse/src/ops/serial/cs.rs
+++ b/nalgebra-sparse/src/ops/serial/cs.rs
@ -20,6 +20,51 @@ fn spmm_cs_unexpected_entry() -> OperationError {
 /// reversed (since transpose(AB) = transpose(B) * transpose(A) and CSC(A) = transpose(CSR(A)).
 ///
 /// We assume here that the matrices have already been verified to be dimensionally compatible.
+pub fn spmm_cs_prealloc_unchecked<T>(
+    beta: T,
+    c: &mut CsMatrix<T>,
+    alpha: T,
+    a: &CsMatrix<T>,
+    b: &CsMatrix<T>,
+) -> Result<(), OperationError>
+where
+    T: Scalar + ClosedAdd + ClosedMul + Zero + One,
+{
+    assert_eq!(c.pattern().major_dim(), a.pattern().major_dim());
+    assert_eq!(c.pattern().minor_dim(), b.pattern().minor_dim());
+    let some_val = Zero::zero();
+    let mut scratchpad_values: Vec<T> = vec![some_val; b.pattern().minor_dim()];
+    for i in 0..c.pattern().major_dim() {
+        let a_lane_i = a.get_lane(i).unwrap();
+
+        let mut c_lane_i = c.get_lane_mut(i).unwrap();
+
+        for (&k, a_ik) in a_lane_i.minor_indices().iter().zip(a_lane_i.values()) {
+            let b_lane_k = b.get_lane(k).unwrap();
+            let alpha_aik = alpha.clone() * a_ik.clone();
+            for (j, b_kj) in b_lane_k.minor_indices().iter().zip(b_lane_k.values()) {
+                // use a dense scatter vector to accumulate non-zeros quickly
+                unsafe {
+                    *scratchpad_values.get_unchecked_mut(*j) += alpha_aik.clone() * b_kj.clone();
+                }
+            }
+        }
+
+        //Get indices from C pattern and gather from the dense scratchpad_values
+        let (indices, values) = c_lane_i.indices_and_values_mut();
+        values
+            .iter_mut()
+            .zip(indices)
+            .for_each(|(output_ref, index)| unsafe {
+                *output_ref = beta.clone() * output_ref.clone()
+                    + scratchpad_values.get_unchecked(*index).clone();
+                *scratchpad_values.get_unchecked_mut(*index) = Zero::zero();
+            });
+    }
+
+    Ok(())
+}
+
 pub fn spmm_cs_prealloc<T>(
    beta: T,
    c: &mut CsMatrix<T>,
--- a/nalgebra-sparse/src/ops/serial/csc.rs
+++ b/nalgebra-sparse/src/ops/serial/csc.rs
@ -1,5 +1,7 @@
 use crate::csc::CscMatrix;
-use crate::ops::serial::cs::{spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc};
+use crate::ops::serial::cs::{
+    spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc, spmm_cs_prealloc_unchecked,
+};
 use crate::ops::serial::{OperationError, OperationErrorKind};
 use crate::ops::Op;
 use nalgebra::{ClosedAdd, ClosedMul, DMatrixSlice, DMatrixSliceMut, RealField, Scalar};
@ -83,14 +85,65 @@ where
 {
    assert_compatible_spmm_dims!(c, a, b);

-    use Op::{NoOp, Transpose};
+    use Op::NoOp;

    match (&a, &b) {
        (NoOp(ref a), NoOp(ref b)) => {
            // Note: We have to reverse the order for CSC matrices
            spmm_cs_prealloc(beta, &mut c.cs, alpha, &b.cs, &a.cs)
        }
-        _ => {
+        _ => spmm_csc_transposed(beta, c, alpha, a, b, spmm_csc_prealloc),
+    }
+}
+
+/// Faster sparse-sparse matrix multiplication, `C <- beta * C + alpha * op(A) * op(B)`.
+/// This will not return an error even if the patterns don't match.
+/// Should be used for situations where pattern creation immediately preceeds multiplication.
+///
+/// Panics if the dimensions of the matrices involved are not compatible with the expression.
+pub fn spmm_csc_prealloc_unchecked<T>(
+    beta: T,
+    c: &mut CscMatrix<T>,
+    alpha: T,
+    a: Op<&CscMatrix<T>>,
+    b: Op<&CscMatrix<T>>,
+) -> Result<(), OperationError>
+where
+    T: Scalar + ClosedAdd + ClosedMul + Zero + One,
+{
+    assert_compatible_spmm_dims!(c, a, b);
+
+    use Op::NoOp;
+
+    match (&a, &b) {
+        (NoOp(ref a), NoOp(ref b)) => {
+            // Note: We have to reverse the order for CSC matrices
+            spmm_cs_prealloc_unchecked(beta, &mut c.cs, alpha, &b.cs, &a.cs)
+        }
+        _ => spmm_csc_transposed(beta, c, alpha, a, b, spmm_csc_prealloc_unchecked),
+    }
+}
+
+fn spmm_csc_transposed<T, F>(
+    beta: T,
+    c: &mut CscMatrix<T>,
+    alpha: T,
+    a: Op<&CscMatrix<T>>,
+    b: Op<&CscMatrix<T>>,
+    spmm_kernel: F,
+) -> Result<(), OperationError>
+where
+    T: Scalar + ClosedAdd + ClosedMul + Zero + One,
+    F: Fn(
+        T,
+        &mut CscMatrix<T>,
+        T,
+        Op<&CscMatrix<T>>,
+        Op<&CscMatrix<T>>,
+    ) -> Result<(), OperationError>,
+{
+    use Op::{NoOp, Transpose};
+
    // Currently we handle transposition by explicitly precomputing transposed matrices
    // and calling the operation again without transposition
    let a_ref: &CscMatrix<T> = a.inner_ref();
@ -101,15 +154,10 @@ where
            (NoOp(_), NoOp(_)) => unreachable!(),
            (Transpose(ref a), NoOp(_)) => (Owned(a.transpose()), Borrowed(b_ref)),
            (NoOp(_), Transpose(ref b)) => (Borrowed(a_ref), Owned(b.transpose())),
-                    (Transpose(ref a), Transpose(ref b)) => {
-                        (Owned(a.transpose()), Owned(b.transpose()))
-                    }
+            (Transpose(ref a), Transpose(ref b)) => (Owned(a.transpose()), Owned(b.transpose())),
        }
    };
-
-            spmm_csc_prealloc(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
-        }
-    }
+    spmm_kernel(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
 }

 /// Solve the lower triangular system `op(L) X = B`.
--- a/nalgebra-sparse/src/ops/serial/csr.rs
+++ b/nalgebra-sparse/src/ops/serial/csr.rs
@ -1,5 +1,7 @@
 use crate::csr::CsrMatrix;
-use crate::ops::serial::cs::{spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc};
+use crate::ops::serial::cs::{
+    spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc, spmm_cs_prealloc_unchecked,
+};
 use crate::ops::serial::OperationError;
 use crate::ops::Op;
 use nalgebra::{ClosedAdd, ClosedMul, DMatrixSlice, DMatrixSliceMut, Scalar};
@ -77,15 +79,63 @@ where
 {
    assert_compatible_spmm_dims!(c, a, b);

-    use Op::{NoOp, Transpose};
+    use Op::NoOp;

    match (&a, &b) {
        (NoOp(ref a), NoOp(ref b)) => spmm_cs_prealloc(beta, &mut c.cs, alpha, &a.cs, &b.cs),
-        _ => {
+        _ => spmm_csr_transposed(beta, c, alpha, a, b, spmm_csr_prealloc),
+    }
+}
+
+/// Faster sparse-sparse matrix multiplication, `C <- beta * C + alpha * op(A) * op(B)`.
+/// This will not return an error even if the patterns don't match.
+/// Should be used for situations where pattern creation immediately preceeds multiplication.
+///
+/// Panics if the dimensions of the matrices involved are not compatible with the expression.
+pub fn spmm_csr_prealloc_unchecked<T>(
+    beta: T,
+    c: &mut CsrMatrix<T>,
+    alpha: T,
+    a: Op<&CsrMatrix<T>>,
+    b: Op<&CsrMatrix<T>>,
+) -> Result<(), OperationError>
+where
+    T: Scalar + ClosedAdd + ClosedMul + Zero + One,
+{
+    assert_compatible_spmm_dims!(c, a, b);
+
+    use Op::NoOp;
+
+    match (&a, &b) {
+        (NoOp(ref a), NoOp(ref b)) => {
+            spmm_cs_prealloc_unchecked(beta, &mut c.cs, alpha, &a.cs, &b.cs)
+        }
+        _ => spmm_csr_transposed(beta, c, alpha, a, b, spmm_csr_prealloc_unchecked),
+    }
+}
+
+fn spmm_csr_transposed<T, F>(
+    beta: T,
+    c: &mut CsrMatrix<T>,
+    alpha: T,
+    a: Op<&CsrMatrix<T>>,
+    b: Op<&CsrMatrix<T>>,
+    spmm_kernel: F,
+) -> Result<(), OperationError>
+where
+    T: Scalar + ClosedAdd + ClosedMul + Zero + One,
+    F: Fn(
+        T,
+        &mut CsrMatrix<T>,
+        T,
+        Op<&CsrMatrix<T>>,
+        Op<&CsrMatrix<T>>,
+    ) -> Result<(), OperationError>,
+{
+    use Op::{NoOp, Transpose};
+
    // Currently we handle transposition by explicitly precomputing transposed matrices
    // and calling the operation again without transposition
-            // TODO: At least use workspaces to allow control of allocations. Maybe
-            // consider implementing certain patterns (like A^T * B) explicitly
    let a_ref: &CsrMatrix<T> = a.inner_ref();
    let b_ref: &CsrMatrix<T> = b.inner_ref();
    let (a, b) = {
@ -94,13 +144,8 @@ where
            (NoOp(_), NoOp(_)) => unreachable!(),
            (Transpose(ref a), NoOp(_)) => (Owned(a.transpose()), Borrowed(b_ref)),
            (NoOp(_), Transpose(ref b)) => (Borrowed(a_ref), Owned(b.transpose())),
-                    (Transpose(ref a), Transpose(ref b)) => {
-                        (Owned(a.transpose()), Owned(b.transpose()))
-                    }
+            (Transpose(ref a), Transpose(ref b)) => (Owned(a.transpose()), Owned(b.transpose())),
        }
    };
-
-            spmm_csr_prealloc(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
-        }
-    }
+    spmm_kernel(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
 }
--- a/nalgebra-sparse/tests/unit_tests/coo.rs
+++ b/nalgebra-sparse/tests/unit_tests/coo.rs
@ -226,6 +226,29 @@ fn coo_push_valid_entries() {
    );
 }

+#[test]
+fn coo_clear_triplets_valid_entries() {
+    let mut coo = CooMatrix::new(3, 3);
+
+    coo.push(0, 0, 1);
+    coo.push(0, 0, 2);
+    coo.push(2, 2, 3);
+    assert_eq!(
+        coo.triplet_iter().collect::<Vec<_>>(),
+        vec![(0, 0, &1), (0, 0, &2), (2, 2, &3)]
+    );
+    coo.clear_triplets();
+    assert_eq!(coo.triplet_iter().collect::<Vec<_>>(), vec![]);
+    // making sure everyhting works after clearing
+    coo.push(0, 0, 1);
+    coo.push(0, 0, 2);
+    coo.push(2, 2, 3);
+    assert_eq!(
+        coo.triplet_iter().collect::<Vec<_>>(),
+        vec![(0, 0, &1), (0, 0, &2), (2, 2, &3)]
+    );
+}
+
 #[test]
 fn coo_push_out_of_bounds_entries() {
    {
--- a/nalgebra-sparse/tests/unit_tests/matrix_market.rs
+++ b/nalgebra-sparse/tests/unit_tests/matrix_market.rs
@ -1,12 +1,21 @@
 use matrixcompare::assert_matrix_eq;
-use nalgebra::dmatrix;
+use nalgebra::matrix;
 use nalgebra::Complex;
-use nalgebra_sparse::io::load_coo_from_matrix_market_str;
+use nalgebra_sparse::io::{
+    load_coo_from_matrix_market_file, load_coo_from_matrix_market_str, save_to_matrix_market_file,
+    save_to_matrix_market_str,
+};
+use nalgebra_sparse::proptest::coo_no_duplicates;
 use nalgebra_sparse::CooMatrix;
+use proptest::prelude::*;
+use tempfile::tempdir;
+
+type C64 = Complex<f64>;
+type C32 = Complex<f32>;

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_sparse_real_general_empty() {
+fn test_matrixmarket_load_sparse_real_general_empty() {
    // Test several valid zero-shapes of a sparse matrix
    let shapes = vec![ (0, 0), (1, 0), (0, 1) ];
    let strings: Vec<String> = shapes
@ -24,7 +33,7 @@ fn test_matrixmarket_sparse_real_general_empty() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_dense_real_general_empty() {
+fn test_matrixmarket_load_dense_real_general_empty() {
    // Test several valid zero-shapes of a dense matrix
    let shapes = vec![ (0, 0), (1, 0), (0, 1) ];
    let strings: Vec<String> = shapes
@ -42,7 +51,7 @@ fn test_matrixmarket_dense_real_general_empty() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_sparse_real_general() {
+fn test_matrixmarket_load_sparse_real_general() {
    let file_str = r#"
 %%MatrixMarket matrix CoOrdinate real general
 % This is also an example of free-format features.
@ -89,7 +98,7 @@ fn test_matrixmarket_sparse_real_general() {
    5     5   1.200e+01
 "#;
    let sparse_mat = load_coo_from_matrix_market_str::<f32>(file_str).unwrap();
-    let expected = dmatrix![
+    let expected = matrix![
        1.0,   0.0,    0.0,    6.0,    0.0;
        0.0,  10.5,    0.0,    0.0,    0.0;
        0.0,   0.0,  0.015,    0.0,    0.0;
@ -101,7 +110,7 @@ fn test_matrixmarket_sparse_real_general() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_sparse_int_symmetric() {
+fn test_matrixmarket_load_sparse_int_symmetric() {
    let file_str = r#"
 %%MatrixMarket matrix coordinate integer symmetric
 %
@ -117,7 +126,7 @@ fn test_matrixmarket_sparse_int_symmetric() {
    5  5  55
 "#;
    let sparse_mat = load_coo_from_matrix_market_str::<i128>(file_str).unwrap();
-    let expected = dmatrix![
+    let expected = matrix![
         11,  0,  0,   0, -15;
          0, 22, 23,  24,   0;
          0, 23, 33,   0,  35;
@ -129,7 +138,7 @@ fn test_matrixmarket_sparse_int_symmetric() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_sparse_complex_hermitian() {
+fn test_matrixmarket_load_sparse_complex_hermitian() {
    let file_str = r#"
 %%MatrixMarket matrix coordinate complex hermitian
 %
@ -144,19 +153,19 @@ fn test_matrixmarket_sparse_complex_hermitian() {

 "#;
    let sparse_mat = load_coo_from_matrix_market_str::<Complex<f64>>(file_str).unwrap();
-    let expected = dmatrix![
-        Complex::<f64>{re:1.0,im:0.0}, Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:0.0,im:0.0},Complex::<f64>{re:0.0,im:0.0};
-        Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:10.5,im:0.0}, Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:250.5,im:-22.22},Complex::<f64>{re:0.0,im:0.0};
-        Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:0.015,im:0.0}, Complex::<f64>{re:0.0,im:0.0},Complex::<f64>{re:0.0,im:0.0};
-        Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:250.5,im:22.22}, Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:-280.0,im:0.0},Complex::<f64>{re:0.0,im:-33.32};
-        Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:0.0,im:0.0}, Complex::<f64>{re:0.0,im:33.32},Complex::<f64>{re:12.0,im:0.0};
+    let expected = matrix![
+        C64{re:1.0,im:0.0}, C64{re:0.0,im:0.0}, C64{re:0.0,im:0.0}, C64{re:0.0,im:0.0},C64{re:0.0,im:0.0};
+        C64{re:0.0,im:0.0}, C64{re:10.5,im:0.0}, C64{re:0.0,im:0.0}, C64{re:250.5,im:-22.22},C64{re:0.0,im:0.0};
+        C64{re:0.0,im:0.0}, C64{re:0.0,im:0.0}, C64{re:0.015,im:0.0}, C64{re:0.0,im:0.0},C64{re:0.0,im:0.0};
+        C64{re:0.0,im:0.0}, C64{re:250.5,im:22.22}, C64{re:0.0,im:0.0}, C64{re:-280.0,im:0.0},C64{re:0.0,im:-33.32};
+        C64{re:0.0,im:0.0}, C64{re:0.0,im:0.0}, C64{re:0.0,im:0.0}, C64{re:0.0,im:33.32},C64{re:12.0,im:0.0};
    ];
    assert_matrix_eq!(sparse_mat, expected);
 }

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_sparse_real_skew() {
+fn test_matrixmarket_load_sparse_real_skew() {
    let file_str = r#"
 %%MatrixMarket matrix coordinate real skew-symmetric
 %
@ -167,7 +176,7 @@ fn test_matrixmarket_sparse_real_skew() {
    5  3  -35.0
 "#;
    let sparse_mat = load_coo_from_matrix_market_str::<f64>(file_str).unwrap();
-    let expected = dmatrix![
+    let expected = matrix![
      0.0,    0.0,   0.0,   0.0,  15.0;
      0.0,    0.0,  23.0,  24.0,   0.0;
      0.0,  -23.0,   0.0,   0.0,  35.0;
@ -179,7 +188,7 @@ fn test_matrixmarket_sparse_real_skew() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_sparse_pattern_general() {
+fn test_matrixmarket_load_sparse_pattern_general() {
    let file_str = r#"
 %%MatrixMarket matrix coordinate pattern general
 %
@ -198,10 +207,10 @@ fn test_matrixmarket_sparse_pattern_general() {
    let pattern_matrix = load_coo_from_matrix_market_str::<()>(file_str).unwrap();
    let nrows = pattern_matrix.nrows();
    let ncols = pattern_matrix.ncols();
-    let (row_idx, col_idx, val) = pattern_matrix.disassemble();
+    let (row_idx, col_idx, val) = pattern_matrix.clone().disassemble();
    let values = vec![1; val.len()];
    let sparse_mat = CooMatrix::try_from_triplets(nrows, ncols, row_idx, col_idx, values).unwrap();
-    let expected = dmatrix![
+    let expected = matrix![
        1, 0, 0, 0, 1;
        0, 0, 1, 1, 0;
        0, 1, 0, 0, 1;
@ -213,7 +222,7 @@ fn test_matrixmarket_sparse_pattern_general() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_dense_real_general() {
+fn test_matrixmarket_load_dense_real_general() {
    let file_str = r#"
 %%MatrixMarket matrix array real general
 %
@ -233,7 +242,7 @@ fn test_matrixmarket_dense_real_general() {

 "#;
    let sparse_mat = load_coo_from_matrix_market_str::<f32>(file_str).unwrap();
-    let expected = dmatrix![
+    let expected = matrix![
        1.0, 5.0,  9.0;
        2.0, 6.0, 10.0;
        3.0, 7.0, 11.0;
@ -244,7 +253,7 @@ fn test_matrixmarket_dense_real_general() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_dense_real_symmetric() {
+fn test_matrixmarket_load_dense_real_symmetric() {
    let file_str = r#"
 %%MatrixMarket matrix array real symmetric
 %
@ -262,7 +271,7 @@ fn test_matrixmarket_dense_real_symmetric() {

 "#;
    let sparse_mat = load_coo_from_matrix_market_str::<f32>(file_str).unwrap();
-    let expected = dmatrix![
+    let expected = matrix![
        1.0, 2.0, 3.0,  4.0;
        2.0, 5.0, 6.0,  7.0;
        3.0, 6.0, 8.0,  9.0;
@ -273,7 +282,7 @@ fn test_matrixmarket_dense_real_symmetric() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_dense_complex_hermitian() {
+fn test_matrixmarket_load_dense_complex_hermitian() {
    let file_str = r#"
 %%MatrixMarket matrix array complex hermitian
 %
@ -290,19 +299,19 @@ fn test_matrixmarket_dense_complex_hermitian() {
 10.0 0.0

 "#;
-    let sparse_mat = load_coo_from_matrix_market_str::<Complex<f64>>(file_str).unwrap();
-    let expected = dmatrix![
-        Complex::<f64>{re:1.0,im:0.0}, Complex::<f64>{re:2.0,im:-2.0} ,Complex::<f64>{re:3.0,im:-3.0} ,Complex::<f64>{re:4.0,im:-4.0};
-        Complex::<f64>{re:2.0,im:2.0}, Complex::<f64>{re:5.0,im:0.0} ,Complex::<f64>{re:6.0,im:-6.0} ,Complex::<f64>{re:7.0,im:-7.0};
-        Complex::<f64>{re:3.0,im:3.0}, Complex::<f64>{re:6.0,im:6.0} ,Complex::<f64>{re:8.0,im:0.0} ,Complex::<f64>{re:9.0,im:-9.0};
-        Complex::<f64>{re:4.0,im:4.0}, Complex::<f64>{re:7.0,im:7.0} ,Complex::<f64>{re:9.0,im:9.0} ,Complex::<f64>{re:10.0,im:0.0};
+    let sparse_mat = load_coo_from_matrix_market_str::<C64>(file_str).unwrap();
+    let expected = matrix![
+        C64{re:1.0,im:0.0}, C64{re:2.0,im:-2.0} ,C64{re:3.0,im:-3.0} ,C64{re:4.0,im:-4.0};
+        C64{re:2.0,im:2.0}, C64{re:5.0,im:0.0} ,C64{re:6.0,im:-6.0} ,C64{re:7.0,im:-7.0};
+        C64{re:3.0,im:3.0}, C64{re:6.0,im:6.0} ,C64{re:8.0,im:0.0} ,C64{re:9.0,im:-9.0};
+        C64{re:4.0,im:4.0}, C64{re:7.0,im:7.0} ,C64{re:9.0,im:9.0} ,C64{re:10.0,im:0.0};
    ];
    assert_matrix_eq!(sparse_mat, expected);
 }

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_dense_int_skew() {
+fn test_matrixmarket_load_dense_int_skew() {
    let file_str = r#"
 %%MatrixMarket matrix array integer skew-symmetric
 %
@ -315,7 +324,7 @@ fn test_matrixmarket_dense_int_skew() {
 6
 "#;
    let sparse_mat = load_coo_from_matrix_market_str::<i32>(file_str).unwrap();
-    let expected = dmatrix![
+    let expected = matrix![
        0,-1,-2,-3;
        1, 0,-4,-5;
        2, 4, 0,-6;
@ -326,7 +335,7 @@ fn test_matrixmarket_dense_int_skew() {

 #[test]
 #[rustfmt::skip]
-fn test_matrixmarket_dense_complex_general() {
+fn test_matrixmarket_load_dense_complex_general() {
    let file_str = r#"
 %%MatrixMarket matrix array complex general
 %
@ -336,10 +345,124 @@ fn test_matrixmarket_dense_complex_general() {
 1 0
 1 0
 "#;
-    let sparse_mat = load_coo_from_matrix_market_str::<Complex<f32>>(file_str).unwrap();
-    let expected = dmatrix![
-        Complex::<f32>{re:1.0,im:0.0},Complex::<f32>{re:1.0,im:0.0};
-        Complex::<f32>{re:1.0,im:0.0},Complex::<f32>{re:1.0,im:0.0};
+    let sparse_mat = load_coo_from_matrix_market_str::<C32>(file_str).unwrap();
+    let expected = matrix![
+        C32{re:1.0,im:0.0},C32{re:1.0,im:0.0};
+        C32{re:1.0,im:0.0},C32{re:1.0,im:0.0};
    ];
    assert_matrix_eq!(sparse_mat, expected);
 }
+
+#[test]
+#[rustfmt::skip]
+fn test_matrixmarket_write_real(){
+    let dense_matrix = matrix![
+        1.0, 2.0, 3.0;
+        2.0, 0.0, 3.0;
+        ];
+    let row_indices = vec![0,1,0,0,1];
+    let col_indices = vec![0,0,1,2,2];
+    let values = vec![1.0,2.0,2.0,3.0,3.0];
+    let coo_matrix = CooMatrix::try_from_triplets(2, 3, row_indices, col_indices, values).unwrap();
+    assert_matrix_eq!(dense_matrix,coo_matrix);
+    let expected = r#"%%matrixmarket matrix coordinate real general
+% matrixmarket file generated by nalgebra-sparse.
+2 3 5
+1 1 1
+2 1 2
+1 2 2
+1 3 3
+2 3 3
+"#;
+    let matrixmarket_str = save_to_matrix_market_str(&coo_matrix);
+    assert_eq!(matrixmarket_str,expected);
+}
+
+#[test]
+fn test_matrixmarket_write_int() {
+    let dense_matrix = matrix![
+    1,2,3;
+    2,0,3;
+    ];
+    let row_indices = vec![0, 1, 0, 0, 1];
+    let col_indices = vec![0, 0, 1, 2, 2];
+    let values = vec![1, 2, 2, 3, 3];
+    let coo_matrix = CooMatrix::try_from_triplets(2, 3, row_indices, col_indices, values).unwrap();
+    assert_matrix_eq!(dense_matrix, coo_matrix);
+    let expected = r#"%%matrixmarket matrix coordinate integer general
+% matrixmarket file generated by nalgebra-sparse.
+2 3 5
+1 1 1
+2 1 2
+1 2 2
+1 3 3
+2 3 3
+"#;
+    let matrixmarket_str = save_to_matrix_market_str(&coo_matrix);
+    assert_eq!(matrixmarket_str, expected);
+}
+
+#[test]
+fn test_matrixmarket_write_pattern() {
+    let row_indices = vec![0, 1, 0, 0, 1];
+    let col_indices = vec![0, 0, 1, 2, 2];
+    let values = vec![(), (), (), (), ()];
+    let coo_matrix = CooMatrix::try_from_triplets(2, 3, row_indices, col_indices, values).unwrap();
+    let expected = r#"%%matrixmarket matrix coordinate pattern general
+% matrixmarket file generated by nalgebra-sparse.
+2 3 5
+1 1 
+2 1 
+1 2 
+1 3 
+2 3 
+"#;
+    let matrixmarket_str = save_to_matrix_market_str(&coo_matrix);
+    assert_eq!(matrixmarket_str, expected);
+}
+
+#[test]
+fn test_matrixmarket_write_complex() {
+    let row_indices = vec![0, 1, 0, 0, 1];
+    let col_indices = vec![0, 0, 1, 2, 2];
+    let values = vec![
+        C64 { re: 1.0, im: 2.0 },
+        C64 { re: 2.0, im: 3.0 },
+        C64 { re: 3.0, im: 4.0 },
+        C64 { re: 4.0, im: 5.0 },
+        C64 { re: 5.0, im: 6.0 },
+    ];
+    let coo_matrix = CooMatrix::try_from_triplets(2, 3, row_indices, col_indices, values).unwrap();
+    let expected = r#"%%matrixmarket matrix coordinate complex general
+% matrixmarket file generated by nalgebra-sparse.
+2 3 5
+1 1 1 2
+2 1 2 3
+1 2 3 4
+1 3 4 5
+2 3 5 6
+"#;
+    let matrixmarket_str = save_to_matrix_market_str(&coo_matrix);
+    assert_eq!(matrixmarket_str, expected);
+}
+
+proptest! {
+    #[test]
+    fn coo_matrix_market_roundtrip_str(coo in coo_no_duplicates(-10 ..= 10, 0 ..= 10, 0..= 10, 100)) {
+        let generated_matrixmarket_string = save_to_matrix_market_str(&coo);
+        let generated_matrix = load_coo_from_matrix_market_str(&generated_matrixmarket_string).unwrap();
+        assert_matrix_eq!(generated_matrix, coo);
+    }
+}
+
+proptest! {
+    #[test]
+    fn coo_matrix_market_roundtrip_file(coo in coo_no_duplicates(-10 ..= 10, 0 ..= 10, 0..= 10, 100)) {
+        let temp_dir = tempdir().expect("Unable to create temporary directory");
+        let file_path = temp_dir.path().join("temp.mtx");
+        save_to_matrix_market_file(&coo,&file_path).unwrap();
+        let generated_matrix = load_coo_from_matrix_market_file(file_path).unwrap();
+        assert_matrix_eq!(generated_matrix, coo);
+        temp_dir.close().expect("Unable to delete temporary directory");
+    }
+}
--- a/nalgebra-sparse/tests/unit_tests/ops.rs
+++ b/nalgebra-sparse/tests/unit_tests/ops.rs
@ -6,7 +6,8 @@ use nalgebra_sparse::csc::CscMatrix;
 use nalgebra_sparse::csr::CsrMatrix;
 use nalgebra_sparse::ops::serial::{
    spadd_csc_prealloc, spadd_csr_prealloc, spadd_pattern, spmm_csc_dense, spmm_csc_prealloc,
-    spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc, spsolve_csc_lower_triangular,
+    spmm_csc_prealloc_unchecked, spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc,
+    spmm_csr_prealloc_unchecked, spsolve_csc_lower_triangular,
 };
 use nalgebra_sparse::ops::Op;
 use nalgebra_sparse::pattern::SparsityPattern;
@ -543,6 +544,29 @@ proptest! {
        prop_assert_eq!(&c_pattern, c_csr.pattern());
    }

+    #[test]
+    fn spmm_csr_prealloc_unchecked_test(SpmmCsrArgs { c, beta, alpha, a, b }
+        in spmm_csr_prealloc_args_strategy()
+    ) {
+        // Test that we get the expected result by comparing to an equivalent dense operation
+        // (here we give in the C matrix, so the sparsity pattern is essentially fixed)
+        let mut c_sparse = c.clone();
+        spmm_csr_prealloc_unchecked(beta, &mut c_sparse, alpha, a.as_ref(), b.as_ref()).unwrap();
+
+        let mut c_dense = DMatrix::from(&c);
+        let op_a_dense = match a {
+            Op::NoOp(ref a) => DMatrix::from(a),
+            Op::Transpose(ref a) => DMatrix::from(a).transpose(),
+        };
+        let op_b_dense = match b {
+            Op::NoOp(ref b) => DMatrix::from(b),
+            Op::Transpose(ref b) => DMatrix::from(b).transpose(),
+        };
+        c_dense = beta * c_dense + alpha * &op_a_dense * op_b_dense;
+
+        prop_assert_eq!(&DMatrix::from(&c_sparse), &c_dense);
+    }
+
    #[test]
    fn spmm_csr_prealloc_test(SpmmCsrArgs { c, beta, alpha, a, b }
        in spmm_csr_prealloc_args_strategy()
@ -705,6 +729,29 @@ proptest! {
        prop_assert_eq!(&DMatrix::from(&c_sparse), &c_dense);
    }

+    #[test]
+    fn spmm_csc_prealloc_unchecked_test(SpmmCscArgs { c, beta, alpha, a, b }
+        in spmm_csc_prealloc_args_strategy()
+    ) {
+        // Test that we get the expected result by comparing to an equivalent dense operation
+        // (here we give in the C matrix, so the sparsity pattern is essentially fixed)
+        let mut c_sparse = c.clone();
+        spmm_csc_prealloc_unchecked(beta, &mut c_sparse, alpha, a.as_ref(), b.as_ref()).unwrap();
+
+        let mut c_dense = DMatrix::from(&c);
+        let op_a_dense = match a {
+            Op::NoOp(ref a) => DMatrix::from(a),
+            Op::Transpose(ref a) => DMatrix::from(a).transpose(),
+        };
+        let op_b_dense = match b {
+            Op::NoOp(ref b) => DMatrix::from(b),
+            Op::Transpose(ref b) => DMatrix::from(b).transpose(),
+        };
+        c_dense = beta * c_dense + alpha * &op_a_dense * op_b_dense;
+
+        prop_assert_eq!(&DMatrix::from(&c_sparse), &c_dense);
+    }
+
    #[test]
    fn spmm_csc_prealloc_panics_on_dim_mismatch(
        (alpha, beta, c, a, b)
--- a/src/base/allocator.rs
+++ b/src/base/allocator.rs
@ -41,6 +41,41 @@ pub trait Allocator<T, R: Dim, C: Dim = U1>: Any + Sized {
        ncols: C,
        iter: I,
    ) -> Self::Buffer;
+
+    #[inline]
+    /// Allocates a buffer initialized with the content of the given row-major order iterator.
+    fn allocate_from_row_iterator<I: IntoIterator<Item = T>>(
+        nrows: R,
+        ncols: C,
+        iter: I,
+    ) -> Self::Buffer {
+        let mut res = Self::allocate_uninit(nrows, ncols);
+        let mut count = 0;
+
+        unsafe {
+            // OK because the allocated buffer is guaranteed to be contiguous.
+            let res_ptr = res.as_mut_slice_unchecked();
+
+            for (k, e) in iter
+                .into_iter()
+                .take(ncols.value() * nrows.value())
+                .enumerate()
+            {
+                let i = k / ncols.value();
+                let j = k % ncols.value();
+                // result[(i, j)] = e;
+                *res_ptr.get_unchecked_mut(i + j * nrows.value()) = MaybeUninit::new(e);
+                count += 1;
+            }
+
+            assert!(
+                count == nrows.value() * ncols.value(),
+                "Matrix init. from row iterator: iterator not long enough."
+            );
+
+            <Self as Allocator<T, R, C>>::assume_init(res)
+        }
+    }
 }

 /// A matrix reallocator. Changes the size of the memory buffer that initially contains (`RFrom` ×
--- a/src/base/construction.rs
+++ b/src/base/construction.rs
@ -86,6 +86,17 @@ where
        Self::from_data(DefaultAllocator::allocate_from_iterator(nrows, ncols, iter))
    }

+    /// Creates a matrix with all its elements filled by an row-major order iterator.
+    #[inline]
+    pub fn from_row_iterator_generic<I>(nrows: R, ncols: C, iter: I) -> Self
+    where
+        I: IntoIterator<Item = T>,
+    {
+        Self::from_data(DefaultAllocator::allocate_from_row_iterator(
+            nrows, ncols, iter,
+        ))
+    }
+
    /// Creates a matrix with its elements filled with the components provided by a slice in
    /// row-major order.
    ///
@ -479,6 +490,36 @@ macro_rules! impl_constructors(
            Self::from_iterator_generic($($gargs, )* iter)
        }

+        /// Creates a matrix or vector with all its elements filled by a row-major iterator.
+        ///
+        /// The output matrix is filled row-by-row.
+        ///
+        /// ## Example
+        /// ```
+        /// # use nalgebra::{Matrix2x3, Vector3, DVector, DMatrix};
+        /// # use std::iter;
+        ///
+        /// let v = Vector3::from_row_iterator((0..3).into_iter());
+        /// // The additional argument represents the vector dimension.
+        /// let dv = DVector::from_row_iterator(3, (0..3).into_iter());
+        /// let m = Matrix2x3::from_row_iterator((0..6).into_iter());
+        /// // The two additional arguments represent the matrix dimensions.
+        /// let dm = DMatrix::from_row_iterator(2, 3, (0..6).into_iter());
+        ///
+        /// // For Vectors from_row_iterator is identical to from_iterator
+        /// assert!(v.x == 0 && v.y == 1 && v.z == 2);
+        /// assert!(dv[0] == 0 && dv[1] == 1 && dv[2] == 2);
+        /// assert!(m.m11 == 0 && m.m12 == 1 && m.m13 == 2 &&
+        ///         m.m21 == 3 && m.m22 == 4 && m.m23 == 5);
+        /// assert!(dm[(0, 0)] == 0 && dm[(0, 1)] == 1 && dm[(0, 2)] == 2 &&
+        ///         dm[(1, 0)] == 3 && dm[(1, 1)] == 4 && dm[(1, 2)] == 5);
+        /// ```
+        #[inline]
+        pub fn from_row_iterator<I>($($args: usize,)* iter: I) -> Self
+            where I: IntoIterator<Item = T> {
+            Self::from_row_iterator_generic($($gargs, )* iter)
+        }
+
        /// Creates a matrix or vector filled with the results of a function applied to each of its
        /// component coordinates.
        ///
--- a/src/base/default_allocator.rs
+++ b/src/base/default_allocator.rs
@ -26,7 +26,7 @@ use std::mem::{ManuallyDrop, MaybeUninit};
 * Allocator.
 *
 */
-/// An allocator based on `GenericArray` and `VecStorage` for statically-sized and dynamically-sized
+/// An allocator based on [`ArrayStorage`] and [`VecStorage`] for statically-sized and dynamically-sized
 /// matrices respectively.
 #[derive(Copy, Clone, Debug)]
 pub struct DefaultAllocator;
--- a/src/base/dimension.rs
+++ b/src/base/dimension.rs
@ -252,14 +252,17 @@ pub trait ToTypenum {
 }

 unsafe impl<const T: usize> Dim for Const<T> {
+    #[inline]
    fn try_to_usize() -> Option<usize> {
        Some(T)
    }

+    #[inline]
    fn value(&self) -> usize {
        T
    }

+    #[inline]
    fn from_usize(dim: usize) -> Self {
        assert_eq!(dim, T);
        Self
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@ -2186,3 +2186,28 @@ where
        }
    }
 }
+
+impl<T, D, S> Unit<Vector<T, D, S>>
+where
+    T: Scalar,
+    D: Dim,
+    S: RawStorage<T, D, U1>,
+{
+    /// Cast the components of `self` to another type.
+    ///
+    /// # Example
+    /// ```
+    /// # use nalgebra::Vector3;
+    /// let v = Vector3::<f64>::y_axis();
+    /// let v2 = v.cast::<f32>();
+    /// assert_eq!(v2, Vector3::<f32>::y_axis());
+    /// ```
+    pub fn cast<T2: Scalar>(self) -> Unit<OVector<T2, D>>
+    where
+        T: Scalar,
+        OVector<T2, D>: SupersetOf<Vector<T, D, S>>,
+        DefaultAllocator: Allocator<T2, D, U1>,
+    {
+        Unit::new_unchecked(crate::convert_ref(self.as_ref()))
+    }
+}
--- a/src/geometry/point_construction.rs
+++ b/src/geometry/point_construction.rs
@ -202,7 +202,7 @@ impl<T: Scalar> Point1<T> {
    /// assert_eq!(p.x, 1.0);
    /// ```
    #[inline]
-    pub fn new(x: T) -> Self {
+    pub const fn new(x: T) -> Self {
        Point {
            coords: Vector1::new(x),
        }
@ -216,7 +216,7 @@ macro_rules! componentwise_constructors_impl(
            #[doc = $doc]
            #[doc = "```"]
            #[inline]
-            pub fn new($($args: T),*) -> Self {
+            pub const fn new($($args: T),*) -> Self {
                Point { coords: $Vector::new($($args),*) }
            }
        }
--- a/src/geometry/rotation_specialization.rs
+++ b/src/geometry/rotation_specialization.rs
@ -17,7 +17,9 @@ use std::ops::Neg;

 use crate::base::dimension::{U1, U2, U3};
 use crate::base::storage::Storage;
-use crate::base::{Matrix2, Matrix3, SMatrix, SVector, Unit, Vector, Vector1, Vector2, Vector3};
+use crate::base::{
+    Matrix2, Matrix3, SMatrix, SVector, Unit, UnitVector3, Vector, Vector1, Vector2, Vector3,
+};

 use crate::geometry::{Rotation2, Rotation3, UnitComplex, UnitQuaternion};

@ -730,9 +732,12 @@ where
        T: RealField,
    {
        if max_iter == 0 {
-            max_iter = usize::max_value();
+            max_iter = usize::MAX;
        }

+        // Using sqrt(eps) ensures we perturb with something larger than eps; clamp to eps to handle the case of eps > 1.0
+        let eps_disturbance = eps.clone().sqrt().max(eps.clone() * eps.clone());
+        let mut perturbation_axes = Vector3::x_axis();
        let mut rot = guess.into_inner();

        for _ in 0..max_iter {
@ -748,10 +753,36 @@ where
            if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, eps.clone()) {
                rot = Rotation3::from_axis_angle(&axis, angle) * rot;
            } else {
+                // Check if stuck in a maximum w.r.t. the norm (m - rot).norm()
+                let mut perturbed = rot.clone();
+                let norm_squared = (m - &rot).norm_squared();
+                let mut new_norm_squared: T;
+
+                // Perturb until the new norm is significantly different
+                loop {
+                    perturbed *=
+                        Rotation3::from_axis_angle(&perturbation_axes, eps_disturbance.clone());
+                    new_norm_squared = (m - &perturbed).norm_squared();
+                    if abs_diff_ne!(
+                        norm_squared,
+                        new_norm_squared,
+                        epsilon = T::default_epsilon()
+                    ) {
                        break;
                    }
                }

+                // If new norm is larger, it's a minimum
+                if norm_squared < new_norm_squared {
+                    break;
+                }
+
+                // If not, continue from perturbed rotation, but use a different axes for the next perturbation
+                perturbation_axes = UnitVector3::new_unchecked(perturbation_axes.yzx());
+                rot = perturbed;
+            }
+        }
+
        Self::from_matrix_unchecked(rot)
    }

--- a/src/geometry/unit_complex_construction.rs
+++ b/src/geometry/unit_complex_construction.rs
@ -131,10 +131,11 @@ where
    ///
    /// # Example
    /// ```
+    /// #[macro_use] extern crate approx;
    /// # use nalgebra::UnitComplex;
    /// let c = UnitComplex::new(1.0f64);
    /// let c2 = c.cast::<f32>();
-    /// assert_eq!(c2, UnitComplex::new(1.0f32));
+    /// assert_relative_eq!(c2, UnitComplex::new(1.0f32));
    /// ```
    pub fn cast<To: Scalar>(self) -> UnitComplex<To>
    where
--- a/src/lib.rs
+++ b/src/lib.rs
@ -46,28 +46,34 @@ fn main() {
 **nalgebra** is meant to be a general-purpose, low-dimensional, linear algebra library, with
 an optimized set of tools for computer graphics and physics. Those features include:

-* A single parametrizable type `Matrix` for vectors, (square or rectangular) matrices, and slices
-  with dimensions known either at compile-time (using type-level integers) or at runtime.
+* A single parametrizable type [`Matrix`](Matrix) for vectors, (square or rectangular) matrices, and
+  slices with dimensions known either at compile-time (using type-level integers) or at runtime.
 * Matrices and vectors with compile-time sizes are statically allocated while dynamic ones are
  allocated on the heap.
-* Convenient aliases for low-dimensional matrices and vectors: `Vector1` to `Vector6` and
-  `Matrix1x1` to `Matrix6x6`, including rectangular matrices like `Matrix2x5`.
-* Points sizes known at compile time, and convenience aliases: `Point1` to `Point6`.
-* Translation (seen as a transformation that composes by multiplication): `Translation2`,
-  `Translation3`.
-* Rotation matrices: `Rotation2`, `Rotation3`.
-* Quaternions: `Quaternion`, `UnitQuaternion` (for 3D rotation).
-* Unit complex numbers can be used for 2D rotation: `UnitComplex`.
-* Algebraic entities with a norm equal to one: `Unit<T>`, e.g., `Unit<Vector3<f32>>`.
-* Isometries (translation ⨯ rotation): `Isometry2`, `Isometry3`
-* Similarity transformations (translation ⨯ rotation ⨯ uniform scale): `Similarity2`, `Similarity3`.
-* Affine transformations stored as a homogeneous matrix: `Affine2`, `Affine3`.
-* Projective (i.e. invertible) transformations stored as a homogeneous matrix: `Projective2`,
-  `Projective3`.
+* Convenient aliases for low-dimensional matrices and vectors: [`Vector1`](Vector1) to
+  [`Vector6`](Vector6) and [`Matrix1x1`](Matrix1) to [`Matrix6x6`](Matrix6), including rectangular
+  matrices like [`Matrix2x5`](Matrix2x5).
+* Points sizes known at compile time, and convenience aliases: [`Point1`](Point1) to
+  [`Point6`](Point6).
+* Translation (seen as a transformation that composes by multiplication):
+  [`Translation2`](Translation2), [`Translation3`](Translation3).
+* Rotation matrices: [`Rotation2`](Rotation2), [`Rotation3`](Rotation3).
+* Quaternions: [`Quaternion`](Quaternion), [`UnitQuaternion`](UnitQuaternion) (for 3D rotation).
+* Unit complex numbers can be used for 2D rotation: [`UnitComplex`](UnitComplex).
+* Algebraic entities with a norm equal to one: [`Unit<T>`](Unit), e.g., `Unit<Vector3<f32>>`.
+* Isometries (translation ⨯ rotation): [`Isometry2`](Isometry2), [`Isometry3`](Isometry3)
+* Similarity transformations (translation ⨯ rotation ⨯ uniform scale):
+  [`Similarity2`](Similarity2), [`Similarity3`](Similarity3).
+* Affine transformations stored as a homogeneous matrix:
+  [`Affine2`](Affine2), [`Affine3`](Affine3).
+* Projective (i.e. invertible) transformations stored as a homogeneous matrix:
+  [`Projective2`](Projective2), [`Projective3`](Projective3).
 * General transformations that does not have to be invertible, stored as a homogeneous matrix:
-  `Transform2`, `Transform3`.
-* 3D projections for computer graphics: `Perspective3`, `Orthographic3`.
-* Matrix factorizations: `Cholesky`, `QR`, `LU`, `FullPivLU`, `SVD`, `Schur`, `Hessenberg`, `SymmetricEigen`.
+  [`Transform2`](Transform2), [`Transform3`](Transform3).
+* 3D projections for computer graphics: [`Perspective3`](Perspective3),
+  [`Orthographic3`](Orthographic3).
+* Matrix factorizations: [`Cholesky`](Cholesky), [`QR`](QR), [`LU`](LU), [`FullPivLU`](FullPivLU),
+  [`SVD`](SVD), [`Schur`](Schur), [`Hessenberg`](Hessenberg), [`SymmetricEigen`](SymmetricEigen).
 * Insertion and removal of rows of columns of a matrix.
 */

@ -109,8 +115,6 @@ extern crate alloc;
 #[cfg(not(feature = "std"))]
 extern crate core as std;

-#[cfg(feature = "io")]
-extern crate pest;
 #[macro_use]
 #[cfg(feature = "io")]
 extern crate pest_derive;
--- a/src/third_party/glam/common/glam_matrix.rs
+++ b/src/third_party/glam/common/glam_matrix.rs
@ -14,7 +14,7 @@ macro_rules! impl_vec_conversion(
        impl From<$Vec2> for Vector2<$N> {
            #[inline]
            fn from(e: $Vec2) -> Vector2<$N> {
-                (*e.as_ref()).into()
+                <[$N;2]>::from(e).into()
            }
        }

@ -31,7 +31,7 @@ macro_rules! impl_vec_conversion(
        impl From<$Vec3> for Vector3<$N> {
            #[inline]
            fn from(e: $Vec3) -> Vector3<$N> {
-                (*e.as_ref()).into()
+                <[$N;3]>::from(e).into()
            }
        }

@ -48,7 +48,7 @@ macro_rules! impl_vec_conversion(
        impl From<$Vec4> for Vector4<$N> {
            #[inline]
            fn from(e: $Vec4) -> Vector4<$N> {
-                (*e.as_ref()).into()
+                <[$N;4]>::from(e).into()
            }
        }

--- a/src/third_party/glam/common/glam_point.rs
+++ b/src/third_party/glam/common/glam_point.rs
@ -9,7 +9,7 @@ macro_rules! impl_point_conversion(
        impl From<$Vec2> for Point2<$N> {
            #[inline]
            fn from(e: $Vec2) -> Point2<$N> {
-                (*e.as_ref()).into()
+                <[$N;2]>::from(e).into()
            }
        }

@ -23,7 +23,7 @@ macro_rules! impl_point_conversion(
        impl From<$Vec3> for Point3<$N> {
            #[inline]
            fn from(e: $Vec3) -> Point3<$N> {
-                (*e.as_ref()).into()
+                <[$N;3]>::from(e).into()
            }
        }

@ -37,7 +37,7 @@ macro_rules! impl_point_conversion(
        impl From<$Vec4> for Point4<$N> {
            #[inline]
            fn from(e: $Vec4) -> Point4<$N> {
-                (*e.as_ref()).into()
+                <[$N;4]>::from(e).into()
            }
        }

--- a/src/third_party/glam/mod.rs
+++ b/src/third_party/glam/mod.rs
@ -12,3 +12,5 @@ mod v018;
 mod v019;
 #[cfg(feature = "glam020")]
 mod v020;
+#[cfg(feature = "glam021")]
+mod v021;
--- a/src/third_party/glam/v021/mod.rs
+++ b/src/third_party/glam/v021/mod.rs
@ -0,0 +1,18 @@
+#[path = "../common/glam_isometry.rs"]
+mod glam_isometry;
+#[path = "../common/glam_matrix.rs"]
+mod glam_matrix;
+#[path = "../common/glam_point.rs"]
+mod glam_point;
+#[path = "../common/glam_quaternion.rs"]
+mod glam_quaternion;
+#[path = "../common/glam_rotation.rs"]
+mod glam_rotation;
+#[path = "../common/glam_similarity.rs"]
+mod glam_similarity;
+#[path = "../common/glam_translation.rs"]
+mod glam_translation;
+#[path = "../common/glam_unit_complex.rs"]
+mod glam_unit_complex;
+
+pub(self) use glam021 as glam;
--- a/tests/geometry/rotation.rs
+++ b/tests/geometry/rotation.rs
@ -1,4 +1,7 @@
-use na::{Quaternion, RealField, UnitQuaternion, Vector2, Vector3};
+use na::{
+    Matrix3, Quaternion, RealField, Rotation3, UnitQuaternion, UnitVector3, Vector2, Vector3,
+};
+use std::f64::consts::PI;

 #[test]
 fn angle_2() {
@ -16,6 +19,58 @@ fn angle_3() {
    assert_eq!(a.angle(&b), 0.0);
 }

+#[test]
+fn from_rotation_matrix() {
+    // Test degenerate case when from_matrix gets stuck in Identity rotation
+    let identity =
+        Rotation3::from_matrix(&Matrix3::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0));
+    assert_relative_eq!(identity, &Rotation3::identity(), epsilon = 0.001);
+    let rotated_z =
+        Rotation3::from_matrix(&Matrix3::new(1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0));
+    assert_relative_eq!(
+        rotated_z,
+        &Rotation3::from_axis_angle(&UnitVector3::new_unchecked(Vector3::new(1.0, 0.0, 0.0)), PI),
+        epsilon = 0.001
+    );
+    // Test that issue 627 is fixed
+    let m_627 = Matrix3::<f64>::new(-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0);
+    assert_relative_ne!(identity, Rotation3::from_matrix(&m_627), epsilon = 0.01);
+    assert_relative_eq!(
+        Rotation3::from_matrix_unchecked(m_627.clone()),
+        Rotation3::from_matrix(&m_627),
+        epsilon = 0.001
+    );
+    // Test that issue 1078 is fixed
+    let m_1078 = Matrix3::<f64>::new(0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0);
+    assert_relative_ne!(identity, Rotation3::from_matrix(&m_1078), epsilon = 0.01);
+    assert_relative_eq!(
+        Rotation3::from_matrix_unchecked(m_1078.clone()),
+        Rotation3::from_matrix(&m_1078),
+        epsilon = 0.001
+    );
+    // Additional test cases for eps >= 1.0
+    assert_relative_ne!(
+        identity,
+        Rotation3::from_matrix_eps(&m_627, 1.2, 0, Rotation3::identity()),
+        epsilon = 0.6
+    );
+    assert_relative_eq!(
+        Rotation3::from_matrix_unchecked(m_627.clone()),
+        Rotation3::from_matrix_eps(&m_627, 1.2, 0, Rotation3::identity()),
+        epsilon = 0.6
+    );
+    assert_relative_ne!(
+        identity,
+        Rotation3::from_matrix_eps(&m_1078, 1.0, 0, Rotation3::identity()),
+        epsilon = 0.1
+    );
+    assert_relative_eq!(
+        Rotation3::from_matrix_unchecked(m_1078.clone()),
+        Rotation3::from_matrix_eps(&m_1078, 1.0, 0, Rotation3::identity()),
+        epsilon = 0.1
+    );
+}
+
 #[test]
 fn quaternion_euler_angles_issue_494() {
    let quat = UnitQuaternion::from_quaternion(Quaternion::new(