diff --git a/CHANGELOG.md b/CHANGELOG.md
index c00c01fc..2950a1df 100644
--- 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
diff --git a/Cargo.toml b/Cargo.toml
index 732676ec..5b9940df 100644
--- 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
 
diff --git a/nalgebra-sparse/Cargo.toml b/nalgebra-sparse/Cargo.toml
index 7220373e..b9c343e4 100644
--- 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"]
diff --git a/nalgebra-sparse/src/io/matrix_market.rs b/nalgebra-sparse/src/io/matrix_market.rs
index dea284ee..975c586e 100644
--- 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(())
+}
diff --git a/nalgebra-sparse/src/io/mod.rs b/nalgebra-sparse/src/io/mod.rs
index 89b21ffb..1f4cb081 100644
--- 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;
diff --git a/nalgebra-sparse/src/lib.rs b/nalgebra-sparse/src/lib.rs
index 8567261a..c62677c3 100644
--- 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;
diff --git a/nalgebra-sparse/src/ops/impl_std_ops.rs b/nalgebra-sparse/src/ops/impl_std_ops.rs
index 107c38ba..91b6574f 100644
--- 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.
diff --git a/nalgebra-sparse/src/ops/serial/cs.rs b/nalgebra-sparse/src/ops/serial/cs.rs
index 86484053..cc13c168 100644
--- 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>,
diff --git a/nalgebra-sparse/src/ops/serial/csc.rs b/nalgebra-sparse/src/ops/serial/csc.rs
index e5c9ae4e..6a691ef9 100644
--- 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,35 +85,81 @@ 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)
         }
-        _ => {
-            // Currently we handle transposition by explicitly precomputing transposed matrices
-            // and calling the operation again without transposition
-            let a_ref: &CscMatrix<T> = a.inner_ref();
-            let b_ref: &CscMatrix<T> = b.inner_ref();
-            let (a, b) = {
-                use Cow::*;
-                match (&a, &b) {
-                    (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()))
-                    }
-                }
-            };
-
-            spmm_csc_prealloc(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
-        }
+        _ => 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();
+    let b_ref: &CscMatrix<T> = b.inner_ref();
+    let (a, b) = {
+        use Cow::*;
+        match (&a, &b) {
+            (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())),
+        }
+    };
+    spmm_kernel(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
+}
+
 /// Solve the lower triangular system `op(L) X = B`.
 ///
 /// Only the lower triangular part of L is read, and the result is stored in B.
diff --git a/nalgebra-sparse/src/ops/serial/csr.rs b/nalgebra-sparse/src/ops/serial/csr.rs
index fa317bbf..708c81a3 100644
--- 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,30 +79,73 @@ 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),
-        _ => {
-            // 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) = {
-                use Cow::*;
-                match (&a, &b) {
-                    (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()))
-                    }
-                }
-            };
-
-            spmm_csr_prealloc(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
-        }
+        _ => 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
+    let a_ref: &CsrMatrix<T> = a.inner_ref();
+    let b_ref: &CsrMatrix<T> = b.inner_ref();
+    let (a, b) = {
+        use Cow::*;
+        match (&a, &b) {
+            (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())),
+        }
+    };
+    spmm_kernel(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref()))
+}
diff --git a/nalgebra-sparse/tests/unit_tests/matrix_market.rs b/nalgebra-sparse/tests/unit_tests/matrix_market.rs
index 48ff1a78..44e27644 100644
--- 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");
+    }
+}
diff --git a/nalgebra-sparse/tests/unit_tests/ops.rs b/nalgebra-sparse/tests/unit_tests/ops.rs
index f2a02fd8..0a335567 100644
--- 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)
diff --git a/src/base/allocator.rs b/src/base/allocator.rs
index 29286420..bb227d3b 100644
--- 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` ×
diff --git a/src/base/construction.rs b/src/base/construction.rs
index fe4e4b08..009d8c76 100644
--- 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.
         ///
diff --git a/src/geometry/rotation_specialization.rs b/src/geometry/rotation_specialization.rs
index 41405c87..c9197fd6 100644
--- 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,7 +753,33 @@ where
             if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, eps.clone()) {
                 rot = Rotation3::from_axis_angle(&axis, angle) * rot;
             } else {
-                break;
+                // 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;
             }
         }
 
diff --git a/src/geometry/unit_complex_construction.rs b/src/geometry/unit_complex_construction.rs
index ebf4e81d..9536e87f 100644
--- 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
diff --git a/src/lib.rs b/src/lib.rs
index 1ee1a3ba..8da86531 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -109,8 +109,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;
diff --git a/src/third_party/glam/common/glam_matrix.rs b/src/third_party/glam/common/glam_matrix.rs
index fa9f713f..26f58803 100644
--- 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()
             }
         }
 
diff --git a/src/third_party/glam/common/glam_point.rs b/src/third_party/glam/common/glam_point.rs
index b15a6c6d..205986b8 100644
--- 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()
             }
         }
 
diff --git a/src/third_party/glam/mod.rs b/src/third_party/glam/mod.rs
index ae2c4514..06bdcfb7 100644
--- 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;
diff --git a/src/third_party/glam/v021/mod.rs b/src/third_party/glam/v021/mod.rs
new file mode 100644
index 00000000..d0168923
--- /dev/null
+++ 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;
diff --git a/tests/geometry/rotation.rs b/tests/geometry/rotation.rs
index 84bba676..2e8d2482 100644
--- 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(