diff --git a/CHANGELOG.md b/CHANGELOG.md
index eebda429..54457072 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,9 @@ This project adheres to [Semantic Versioning](http://semver.org/).
 
 ## [0.13.0]
 
+The **nalgebra-lapack** crate has been updated. This now includes a broad range
+matrix decompositions using LAPACK bindings.
+
 ### Breaking semantic change
   * The implementation of slicing with steps now matches the documentation.
     Before, step identified the number to add to pass from one column/row index
@@ -19,10 +22,13 @@ This project adheres to [Semantic Versioning](http://semver.org/).
   * The trait `Axpy` takes one additional parameter for the type of `x`.
   * The alias `MatrixNM` is now deprecated. Use `MatrixMN` instead (we
     reordered M and N to be in alphabetical order).
+  * In-place componentwise multiplication and division
+    `.component_mul_mut(...)` and `.component_div_mut(...)` have bee deprecated
+    for a future renaming. Use `.component_mul_assign(...)` and
+    `.component_div_assign(...)` instead.
 
 ### Added
   * `alga::general::Real` is now re-exported by nalgebra.
-  * `.trace()` that computes the trace of a matrix (the sum of its diagonal
     elements.)
   * `::zeros(...)` that creates a matrix filled with zeroes.
   * `::from_partial_diagonal(...)` that creates a matrix from diagonal elements.
@@ -39,16 +45,20 @@ This project adheres to [Semantic Versioning](http://semver.org/).
   * Implements `IntoIterator` for `&Matrix`, `&mut Matrix` and `Matrix`.
   * `.mul_to(...)` multiplies two matrices and stores the result to the given buffer.
   * `.tr_mul_to(...)` left-multiplies `self.transpose()` to another matrix and stores the result to the given buffer.
-  * `.rows_range(...)` that retrieves a reference to a range of rows.
-  * `.rows_range_mut(...)` that retrieves a mutable reference to a range of rows.
-  * `.columns_range(...)` that retrieves a reference to a range of columns.
-  * `.columns_range_mut(...)` that retrieves a mutable reference to a range of columns.
   * `.add_scalar(...)` that adds a scalar to each component of a matrix.
   * `.add_scalar_mut(...)` that adds in-place a scalar to each component of a matrix.
   * `.kronecker(a, b)` computes the kronecker product (i.e. matrix tensor
     product) of two matrices.
+  * `.apply(f)` replaces each component of a matrix with the results of the
+    closure `f` called on each of them.
 
-Matrix decompositions:
+New slicing methods:
+  * `.rows_range(...)` that retrieves a reference to a range of rows.
+  * `.rows_range_mut(...)` that retrieves a mutable reference to a range of rows.
+  * `.columns_range(...)` that retrieves a reference to a range of columns.
+  * `.columns_range_mut(...)` that retrieves a mutable reference to a range of columns.
+
+Matrix decompositions implemented in pure Rust:
   * Cholesky, SVD, LU, QR, Hessenberg, Schur, Symmetric eigendecompositions,
     Bidiagonal, Symmetric tridiagonal
   * Computation of householder reflectors and givens rotations.
diff --git a/nalgebra-lapack/src/cholesky.rs b/nalgebra-lapack/src/cholesky.rs
index 065a6bb4..317b67a4 100644
--- a/nalgebra-lapack/src/cholesky.rs
+++ b/nalgebra-lapack/src/cholesky.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num::Zero;
 use num_complex::Complex;
 
@@ -9,11 +12,25 @@ use na::allocator::Allocator;
 use lapack::fortran as interface;
 
 /// The cholesky decomposion of a symmetric-definite-positive matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D>,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D>,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct Cholesky<N: Scalar, D: Dim>
     where DefaultAllocator: Allocator<N, D, D> {
     l: MatrixN<N, D>
 }
 
+impl<N: Scalar, D: Dim> Copy for Cholesky<N, D>
+    where DefaultAllocator: Allocator<N, D, D>,
+          MatrixN<N, D>: Copy { }
+
 impl<N: CholeskyScalar + Zero, D: Dim> Cholesky<N, D>
     where DefaultAllocator: Allocator<N, D, D> {
 
diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs
index 4121fbc3..4e0cc26e 100644
--- a/nalgebra-lapack/src/eigen.rs
+++ b/nalgebra-lapack/src/eigen.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num::Zero;
 use num_complex::Complex;
 
@@ -12,9 +15,21 @@ use na::allocator::Allocator;
 use lapack::fortran as interface;
 
 /// Eigendecomposition of a real square matrix with real eigenvalues.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct Eigen<N: Scalar, D: Dim>
     where DefaultAllocator: Allocator<N, D> +
-             Allocator<N, D, D> {
+                            Allocator<N, D, D> {
     /// The eigenvalues of the decomposed matrix.
     pub eigenvalues:       VectorN<N, D>,
     /// The (right) eigenvectors of the decomposed matrix.
@@ -23,6 +38,12 @@ pub struct Eigen<N: Scalar, D: Dim>
     pub left_eigenvectors: Option<MatrixN<N, D>>
 }
 
+impl<N: Scalar, D: Dim> Copy for Eigen<N, D>
+    where DefaultAllocator: Allocator<N, D> +
+                            Allocator<N, D, D>,
+          VectorN<N, D>: Copy,
+          MatrixN<N, D>: Copy { }
+
 
 impl<N: EigenScalar + Real, D: Dim> Eigen<N, D>
     where DefaultAllocator: Allocator<N, D, D> +
diff --git a/nalgebra-lapack/src/hessenberg.rs b/nalgebra-lapack/src/hessenberg.rs
index b7576022..ae417445 100644
--- a/nalgebra-lapack/src/hessenberg.rs
+++ b/nalgebra-lapack/src/hessenberg.rs
@@ -11,6 +11,20 @@ use lapack::fortran as interface;
 
 
 /// The Hessenberg decomposition of a general matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, DimDiff<D, U1>>,
+         MatrixN<N, D>: serde::Serialize,
+         VectorN<N, DimDiff<D, U1>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, DimDiff<D, U1>>,
+         MatrixN<N, D>: serde::Deserialize<'de>,
+         VectorN<N, DimDiff<D, U1>>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct Hessenberg<N: Scalar, D: DimSub<U1>>
     where DefaultAllocator: Allocator<N, D, D> +
                             Allocator<N, DimDiff<D, U1>> {
@@ -19,6 +33,12 @@ pub struct Hessenberg<N: Scalar, D: DimSub<U1>>
 }
 
 
+impl<N: Scalar, D: DimSub<U1>> Copy for Hessenberg<N, D>
+    where DefaultAllocator: Allocator<N, D, D> +
+                            Allocator<N, DimDiff<D, U1>>,
+          MatrixN<N, D>: Copy,
+          VectorN<N, DimDiff<D, U1>>: Copy { }
+
 impl<N: HessenbergScalar + Zero, D: DimSub<U1>> Hessenberg<N, D>
         where DefaultAllocator: Allocator<N, D, D> +
                                 Allocator<N, DimDiff<D, U1>> {
diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs
index 066c1cea..c36e7387 100644
--- a/nalgebra-lapack/src/lib.rs
+++ b/nalgebra-lapack/src/lib.rs
@@ -30,20 +30,24 @@
 //! the system installation of netlib without LAPACKE (note the E) or
 //! CBLAS:
 //! 
-//!     sudo apt-get install gfortran libblas3gf liblapack3gf
-//!     export CARGO_FEATURE_SYSTEM_NETLIB=1
-//!     export CARGO_FEATURE_EXCLUDE_LAPACKE=1
-//!     export CARGO_FEATURE_EXCLUDE_CBLAS=1
+//! ```.ignore
+//! sudo apt-get install gfortran libblas3gf liblapack3gf
+//! export CARGO_FEATURE_SYSTEM_NETLIB=1
+//! export CARGO_FEATURE_EXCLUDE_LAPACKE=1
+//! export CARGO_FEATURE_EXCLUDE_CBLAS=1
 //! 
-//!     export CARGO_FEATURES='--no-default-features --features netlib'
-//!     cargo build ${CARGO_FEATURES}
+//! export CARGO_FEATURES='--no-default-features --features netlib'
+//! cargo build ${CARGO_FEATURES}
+//! ```
 //! 
 //! ### Mac OS X
 //! 
 //! On Mac OS X, do this to use Apple's Accelerate framework:
 //! 
-//!     export CARGO_FEATURES='--no-default-features --features accelerate'
-//!     cargo build ${CARGO_FEATURES}
+//! ```.ignore
+//! export CARGO_FEATURES='--no-default-features --features accelerate'
+//! cargo build ${CARGO_FEATURES}
+//! ```
 //! 
 //! [version-img]: https://img.shields.io/crates/v/nalgebra-lapack.svg
 //! [version-url]: https://crates.io/crates/nalgebra-lapack
diff --git a/nalgebra-lapack/src/lu.rs b/nalgebra-lapack/src/lu.rs
index f3b5bdb0..07de2d92 100644
--- a/nalgebra-lapack/src/lu.rs
+++ b/nalgebra-lapack/src/lu.rs
@@ -17,6 +17,20 @@ use lapack::fortran as interface;
 /// * `P` which is a `m * m` permutation matrix.
 ///
 /// Those are such that `M == P * L * U`.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<i32, DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Serialize,
+         PermutationSequence<DimMinimum<R, C>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<i32, DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Deserialize<'de>,
+         PermutationSequence<DimMinimum<R, C>>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct LU<N: Scalar, R: DimMin<C>, C: Dim>
     where DefaultAllocator: Allocator<i32, DimMinimum<R, C>> +
                             Allocator<N, R, C> {
@@ -24,6 +38,12 @@ pub struct LU<N: Scalar, R: DimMin<C>, C: Dim>
     p:  VectorN<i32, DimMinimum<R, C>>
 }
 
+impl<N: Scalar, R: DimMin<C>, C: Dim> Copy for LU<N, R, C>
+    where DefaultAllocator: Allocator<N, R, C> +
+                            Allocator<i32, DimMinimum<R, C>>,
+          MatrixMN<N, R, C>: Copy,
+          VectorN<i32, DimMinimum<R, C>>: Copy { }
+
 impl<N: LUScalar, R: Dim, C: Dim> LU<N, R, C>
     where N: Zero + One,
           R: DimMin<C>,
diff --git a/nalgebra-lapack/src/qr.rs b/nalgebra-lapack/src/qr.rs
index b69452f8..d221b86b 100644
--- a/nalgebra-lapack/src/qr.rs
+++ b/nalgebra-lapack/src/qr.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num_complex::Complex;
 use num::Zero;
 
@@ -11,6 +14,20 @@ use lapack::fortran as interface;
 
 
 /// The QR decomposition of a general matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<N, DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Serialize,
+         VectorN<N, DimMinimum<R, C>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<N, DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Deserialize<'de>,
+         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct QR<N: Scalar, R: DimMin<C>, C: Dim>
     where DefaultAllocator: Allocator<N, R, C> +
                             Allocator<N, DimMinimum<R, C>> {
@@ -18,6 +35,12 @@ pub struct QR<N: Scalar, R: DimMin<C>, C: Dim>
     tau: VectorN<N, DimMinimum<R, C>>
 }
 
+impl<N: Scalar, R: DimMin<C>, C: Dim> Copy for QR<N, R, C>
+    where DefaultAllocator: Allocator<N, R, C> +
+                            Allocator<N, DimMinimum<R, C>>,
+          MatrixMN<N, R, C>: Copy,
+          VectorN<N, DimMinimum<R, C>>: Copy { }
+
 impl<N: QRScalar + Zero, R: DimMin<C>, C: Dim> QR<N, R, C>
     where DefaultAllocator: Allocator<N, R, C> +
                             Allocator<N, R, DimMinimum<R, C>> +
diff --git a/nalgebra-lapack/src/schur.rs b/nalgebra-lapack/src/schur.rs
index 225c3241..8576b336 100644
--- a/nalgebra-lapack/src/schur.rs
+++ b/nalgebra-lapack/src/schur.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num::Zero;
 use num_complex::Complex;
 
@@ -12,9 +15,21 @@ use na::allocator::Allocator;
 use lapack::fortran as interface;
 
 /// Eigendecomposition of a real square matrix with real eigenvalues.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct RealSchur<N: Scalar, D: Dim>
     where DefaultAllocator: Allocator<N, D> +
-    Allocator<N, D, D> {
+                            Allocator<N, D, D> {
 
     re: VectorN<N, D>,
     im: VectorN<N, D>,
@@ -22,6 +37,11 @@ pub struct RealSchur<N: Scalar, D: Dim>
     q:  MatrixN<N, D>
 }
 
+impl<N: Scalar, D: Dim> Copy for RealSchur<N, D>
+    where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+                            MatrixN<N, D>: Copy,
+                            VectorN<N, D>: Copy { }
+
 
 impl<N: RealSchurScalar + Real, D: Dim> RealSchur<N, D>
     where DefaultAllocator: Allocator<N, D, D> +
diff --git a/nalgebra-lapack/src/svd.rs b/nalgebra-lapack/src/svd.rs
index 7f13c5f3..a6587051 100644
--- a/nalgebra-lapack/src/svd.rs
+++ b/nalgebra-lapack/src/svd.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use std::cmp;
 use num::Signed;
 
@@ -11,18 +14,43 @@ use lapack::fortran as interface;
 
 
 /// The SVD decomposition of a general matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
+                           Allocator<N, R, R> +
+                           Allocator<N, C, C>,
+         MatrixN<N, R>: serde::Serialize,
+         MatrixN<N, C>: serde::Serialize,
+         VectorN<N, DimMinimum<R, C>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
+                           Allocator<N, R, R> +
+                           Allocator<N, C, C>,
+         MatrixN<N, R>: serde::Deserialize<'de>,
+         MatrixN<N, C>: serde::Deserialize<'de>,
+         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct SVD<N: Scalar, R: DimMin<C>, C: Dim>
     where DefaultAllocator: Allocator<N, R, R>                 +
                             Allocator<N, DimMinimum<R, C>> +
                             Allocator<N, C, C> {
     /// The left-singular vectors `U` of this SVD.
-    pub u:  MatrixN<N, R>,
+    pub u:  MatrixN<N, R>, // FIXME: should be MatrixMN<N, R, DimMinimum<R, C>>
     /// The right-singular vectors `V^t` of this SVD.
-    pub vt: MatrixN<N, C>,
+    pub vt: MatrixN<N, C>, // FIXME: should be MatrixMN<N, DimMinimum<R, C>, C>
     /// The singular values of this SVD.
     pub singular_values: VectorN<N, DimMinimum<R, C>>
 }
 
+impl<N: Scalar, R: DimMin<C>, C: Dim> Copy for SVD<N, R, C>
+    where DefaultAllocator: Allocator<N, C, C> +
+                            Allocator<N, R, R> +
+                            Allocator<N, DimMinimum<R, C>>,
+          MatrixMN<N, R, R>: Copy,
+          MatrixMN<N, C, C>: Copy,
+          VectorN<N, DimMinimum<R, C>>: Copy { }
 
 /// Trait implemented by floats (`f32`, `f64`) and complex floats (`Complex<f32>`, `Complex<f64>`)
 /// supported by the Singular Value Decompotition.
diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs
index 881fe2f5..16d3f47b 100644
--- a/nalgebra-lapack/src/symmetric_eigen.rs
+++ b/nalgebra-lapack/src/symmetric_eigen.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num::Zero;
 use std::ops::MulAssign;
 
@@ -11,7 +14,21 @@ use na::allocator::Allocator;
 
 use lapack::fortran as interface;
 
-/// SymmetricEigendecomposition of a real square matrix with real eigenvalues.
+/// Eigendecomposition of a real square symmetric matrix with real eigenvalues.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, D>,
+         VectorN<N, D>: serde::Deserialize<'de>,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
+#[derive(Clone, Debug)]
 pub struct SymmetricEigen<N: Scalar, D: Dim>
     where DefaultAllocator: Allocator<N, D> +
     Allocator<N, D, D> {
@@ -23,6 +40,12 @@ pub struct SymmetricEigen<N: Scalar, D: Dim>
 }
 
 
+impl<N: Scalar, D: Dim> Copy for SymmetricEigen<N, D>
+    where DefaultAllocator: Allocator<N, D, D> +
+                            Allocator<N, D>,
+          MatrixN<N, D>: Copy,
+          VectorN<N, D>: Copy { }
+
 impl<N: SymmetricEigenScalar + Real, D: Dim> SymmetricEigen<N, D>
     where DefaultAllocator: Allocator<N, D, D> +
                             Allocator<N, D> {
diff --git a/nalgebra-lapack/tests/linalg/svd.rs b/nalgebra-lapack/tests/linalg/svd.rs
index 7a76d445..9ab7a99e 100644
--- a/nalgebra-lapack/tests/linalg/svd.rs
+++ b/nalgebra-lapack/tests/linalg/svd.rs
@@ -5,10 +5,10 @@ quickcheck!{
     fn svd(m: DMatrix<f64>) -> bool {
         if m.nrows() != 0 && m.ncols() != 0 {
             let svd = SVD::new(m.clone()).unwrap();
-            let sm  = DMatrix::from_partial_diagonal(m.nrows(), m.ncols(), svd.s.as_slice());
+            let sm  = DMatrix::from_partial_diagonal(m.nrows(), m.ncols(), svd.singular_values.as_slice());
 
             let reconstructed_m = &svd.u * sm * &svd.vt;
-            let reconstructed_m2 = svd.matrix();
+            let reconstructed_m2 = svd.recompose();
 
             relative_eq!(reconstructed_m, m, epsilon = 1.0e-7) &&
             relative_eq!(reconstructed_m2, reconstructed_m, epsilon = 1.0e-7)
@@ -20,10 +20,10 @@ quickcheck!{
 
     fn svd_static(m: Matrix3x4<f64>) -> bool {
         let svd = SVD::new(m).unwrap();
-        let sm  = Matrix3x4::from_partial_diagonal(svd.s.as_slice());
+        let sm  = Matrix3x4::from_partial_diagonal(svd.singular_values.as_slice());
 
         let reconstructed_m  = &svd.u * &sm * &svd.vt;
-        let reconstructed_m2 = svd.matrix();
+        let reconstructed_m2 = svd.recompose();
 
         relative_eq!(reconstructed_m, m, epsilon = 1.0e-7) &&
         relative_eq!(reconstructed_m2, m, epsilon = 1.0e-7)
diff --git a/src/core/componentwise.rs b/src/core/componentwise.rs
index 031af3de..65a4d1a9 100644
--- a/src/core/componentwise.rs
+++ b/src/core/componentwise.rs
@@ -33,7 +33,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
 }
 
 macro_rules! component_binop_impl(
-    ($($binop: ident, $binop_mut: ident, $Trait: ident . $binop_assign: ident, $desc:expr, $desc_mut:expr);* $(;)*) => {$(
+    ($($binop: ident, $binop_mut: ident, $binop_assign: ident, $Trait: ident . $op_assign: ident, $desc:expr, $desc_mut:expr);* $(;)*) => {$(
         impl<N: Scalar, R1: Dim, C1: Dim, SA: Storage<N, R1, C1>> Matrix<N, R1, C1, SA> {
             #[doc = $desc]
             #[inline]
@@ -50,7 +50,7 @@ macro_rules! component_binop_impl(
                 for j in 0 .. res.ncols() {
                     for i in 0 .. res.nrows() {
                         unsafe {
-                            res.get_unchecked_mut(i, j).$binop_assign(*rhs.get_unchecked(i, j));
+                            res.get_unchecked_mut(i, j).$op_assign(*rhs.get_unchecked(i, j));
                         }
                     }
                 }
@@ -62,7 +62,7 @@ macro_rules! component_binop_impl(
         impl<N: Scalar, R1: Dim, C1: Dim, SA: StorageMut<N, R1, C1>> Matrix<N, R1, C1, SA> {
             #[doc = $desc_mut]
             #[inline]
-            pub fn $binop_mut<R2, C2, SB>(&mut self, rhs: &Matrix<N, R2, C2, SB>)
+            pub fn $binop_assign<R2, C2, SB>(&mut self, rhs: &Matrix<N, R2, C2, SB>)
                 where N: $Trait,
                       R2: Dim,
                       C2: Dim,
@@ -74,19 +74,31 @@ macro_rules! component_binop_impl(
                 for j in 0 .. self.ncols() {
                     for i in 0 .. self.nrows() {
                         unsafe {
-                            self.get_unchecked_mut(i, j).$binop_assign(*rhs.get_unchecked(i, j));
+                            self.get_unchecked_mut(i, j).$op_assign(*rhs.get_unchecked(i, j));
                         }
                     }
                 }
             }
+
+            #[doc = $desc_mut]
+            #[inline]
+            #[deprecated(note = "This is renamed using the `_assign` sufix instead of the `_mut` suffix.")]
+            pub fn $binop_mut<R2, C2, SB>(&mut self, rhs: &Matrix<N, R2, C2, SB>)
+                where N: $Trait,
+                      R2: Dim,
+                      C2: Dim,
+                      SB: Storage<N, R2, C2>,
+                      ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2> {
+                self.$binop_assign(rhs)
+            }
         }
     )*}
 );
 
 component_binop_impl!(
-    component_mul, component_mul_mut, ClosedMul.mul_assign,
+    component_mul, component_mul_mut, component_mul_assign, ClosedMul.mul_assign,
     "Componentwise matrix multiplication.", "Mutable, componentwise matrix multiplication.";
-    component_div, component_div_mut, ClosedDiv.div_assign,
+    component_div, component_div_mut, component_div_assign, ClosedDiv.div_assign,
     "Componentwise matrix division.", "Mutable, componentwise matrix division.";
     // FIXME: add other operators like bitshift, etc. ?
 );
diff --git a/src/core/edition.rs b/src/core/edition.rs
index baf33d54..da07555e 100644
--- a/src/core/edition.rs
+++ b/src/core/edition.rs
@@ -206,8 +206,8 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         self.remove_fixed_columns::<U1>(i)
     }
 
-    /// Removes a fixed number − `D::dim()` − of consecutive columns from this matrix, starting
-    /// with the `i`-th (included).
+    /// Removes `D::dim()` consecutive columns from this matrix, starting with the `i`-th
+    /// (included).
     #[inline]
     pub fn remove_fixed_columns<D>(self, i: usize) -> MatrixMN<N, R, DimDiff<C, D>>
         where D: DimName,
@@ -271,8 +271,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         self.remove_fixed_rows::<U1>(i)
     }
 
-    /// Removes a fixed number − `D::dim()` − of consecutive rows from this matrix, starting
-    /// with the `i`-th (included).
+    /// Removes `D::dim()` consecutive rows from this matrix, starting with the `i`-th (included).
     #[inline]
     pub fn remove_fixed_rows<D>(self, i: usize) -> MatrixMN<N, DimDiff<R, D>, C>
         where D: DimName,
@@ -328,7 +327,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         self.insert_fixed_columns::<U1>(i, val)
     }
 
-    /// Inserts `D` column filled with `val` at the `i-th` position.
+    /// Inserts `D::dim()` columns filled with `val` starting at the `i-th` position.
     #[inline]
     pub fn insert_fixed_columns<D>(self, i: usize, val: N) -> MatrixMN<N, R, DimSum<C, D>>
         where D: DimName,
@@ -339,7 +338,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         res
     }
 
-    /// Inserts `n` column filled with `val` at the `i-th` position.
+    /// Inserts `n` columns filled with `val` starting at the `i-th` position.
     #[inline]
     pub fn insert_columns(self, i: usize, n: usize, val: N) -> MatrixMN<N, R, Dynamic>
         where C: DimAdd<Dynamic, Output = Dynamic>,
@@ -349,7 +348,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         res
     }
 
-    /// Inserts `ninsert.value()` columns at the `i-th` place of this matrix.
+    /// Inserts `ninsert.value()` columns starting at the `i-th` place of this matrix.
     ///
     /// The added column values are not initialized.
     #[inline]
@@ -388,7 +387,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         self.insert_fixed_rows::<U1>(i, val)
     }
 
-    /// Inserts `D` row filled with `val` at the `i-th` position.
+    /// Inserts `D::dim()` rows filled with `val` starting at the `i-th` position.
     #[inline]
     pub fn insert_fixed_rows<D>(self, i: usize, val: N) -> MatrixMN<N, DimSum<R, D>, C>
         where D: DimName,
@@ -399,7 +398,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         res
     }
 
-    /// Inserts `n` rows filled with `val` at the `i-th` position.
+    /// Inserts `n` rows filled with `val` starting at the `i-th` position.
     #[inline]
     pub fn insert_rows(self, i: usize, n: usize, val: N) -> MatrixMN<N, Dynamic, C>
         where R: DimAdd<Dynamic, Output = Dynamic>,
@@ -412,6 +411,8 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     /// Inserts `ninsert.value()` rows at the `i-th` place of this matrix.
     ///
     /// The added rows values are not initialized.
+    /// This is the generic implementation of `.insert_rows(...)` and
+    /// `.insert_fixed_rows(...)` which have nicer API interfaces.
     #[inline]
     pub unsafe fn insert_rows_generic_uninitialized<D>(self, i: usize, ninsert: D)
         -> MatrixMN<N, DimSum<R, D>, C>
diff --git a/src/core/matrix.rs b/src/core/matrix.rs
index e767488b..769007d4 100644
--- a/src/core/matrix.rs
+++ b/src/core/matrix.rs
@@ -79,11 +79,9 @@ impl<N, R, C, S> Serialize for Matrix<N, R, C, S>
     where N: Scalar,
           R: Dim,
           C: Dim,
-          S: Serialize,
-{
+          S: Serialize, {
     fn serialize<T>(&self, serializer: T) -> Result<T::Ok, T::Error>
-        where T: Serializer
-    {
+        where T: Serializer {
         self.data.serialize(serializer)
     }
 }
@@ -93,8 +91,7 @@ impl<'de, N, R, C, S> Deserialize<'de> for Matrix<N, R, C, S>
     where N: Scalar,
           R: Dim,
           C: Dim,
-          S: Deserialize<'de>,
-{
+          S: Deserialize<'de> {
     fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
         where D: Deserializer<'de>
     {
diff --git a/src/lib.rs b/src/lib.rs
index 8b1ee29b..095c21f1 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -50,7 +50,7 @@ an optimized set of tools for computer graphics and physics. Those features incl
 * 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`.
+  `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`.
@@ -66,7 +66,8 @@ an optimized set of tools for computer graphics and physics. Those features incl
 * General transformations that does not have to be invertible, stored as an homogeneous matrix:
   `Transform2`, `Transform3`.
 * 3D projections for computer graphics: `Perspective3`, `Orthographic3`.
-* Linear algebra decompositions: Cholesky, QR, LU, SVD, Schur, symmetric-eigendecomposition,.
+* Matrix factorizations: `Cholesky`, `QR`, `LU`, `FullPivLU`, `SVD`, `RealSchur`, `Hessenberg`, `SymmetricEigen`.
+* Insertion and removal of rows of columns of a matrix.
 * Implements traits from the [alga](https://crates.io/crates/alga) crate for
   generic programming.
 */
diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs
index 565dee65..2ccbb24a 100644
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use alga::general::Real;
 use core::{Unit, Matrix, MatrixN, MatrixMN, VectorN, DefaultAllocator};
 use dimension::{Dim, DimMin, DimMinimum, DimSub, DimDiff, Dynamic, U1};
@@ -10,6 +13,25 @@ use geometry::Reflection;
 
 
 /// The bidiagonalization of a general matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DimMinimum<R, C>: DimSub<U1>,
+         DefaultAllocator: Allocator<N, R, C>             +
+                           Allocator<N, DimMinimum<R, C>> +
+                           Allocator<N, DimDiff<DimMinimum<R, C>, U1>>,
+         MatrixMN<N, R, C>: serde::Serialize,
+         VectorN<N, DimMinimum<R, C>>: serde::Serialize,
+         VectorN<N, DimDiff<DimMinimum<R, C>, U1>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DimMinimum<R, C>: DimSub<U1>,
+         DefaultAllocator: Allocator<N, R, C>             +
+                           Allocator<N, DimMinimum<R, C>> +
+                           Allocator<N, DimDiff<DimMinimum<R, C>, U1>>,
+         MatrixMN<N, R, C>: serde::Deserialize<'de>,
+         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>,
+         VectorN<N, DimDiff<DimMinimum<R, C>, U1>>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct Bidiagonal<N: Real, R: DimMin<C>, C: Dim>
     where DimMinimum<R, C>: DimSub<U1>,
@@ -35,6 +57,7 @@ impl<N: Real, R: DimMin<C>, C: Dim> Copy for Bidiagonal<N, R, C>
           VectorN<N, DimMinimum<R, C>>: Copy,
           VectorN<N, DimDiff<DimMinimum<R, C>, U1>>: Copy { }
 
+
 impl<N: Real, R: DimMin<C>, C: Dim> Bidiagonal<N, R, C>
     where DimMinimum<R, C>: DimSub<U1>,
           DefaultAllocator: Allocator<N, R, C>             +
@@ -93,7 +116,7 @@ impl<N: Real, R: DimMin<C>, C: Dim> Bidiagonal<N, R, C>
         }
     }
 
-    /// Unpacks thi decomposition into its thrme matrix factors `(U, D, V^t)`.
+    /// Unpacks this decomposition into its three matrix factors `(U, D, V^t)`.
     ///
     /// The decomposed matrix `M` is equal to `U * D * V^t`.
     #[inline]
diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs
index 6037cc9f..7f9f6c0c 100644
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use alga::general::Real;
 
 use core::{DefaultAllocator, MatrixN, MatrixMN, Matrix, SquareMatrix};
@@ -6,7 +9,16 @@ use storage::{Storage, StorageMut};
 use allocator::Allocator;
 use dimension::{Dim, Dynamic, DimSub};
 
-/// The cholesky decomposion of a symmetric-definite-positive matrix.
+/// The Cholesky decomposion of a symmetric-definite-positive matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D>,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D>,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct Cholesky<N: Real, D: Dim>
     where DefaultAllocator: Allocator<N, D, D> {
@@ -19,7 +31,7 @@ impl<N: Real, D: Dim> Copy for Cholesky<N, D>
 
 impl<N: Real, D: DimSub<Dynamic>> Cholesky<N, D>
     where DefaultAllocator: Allocator<N, D, D> {
-    /// Attempts to compute the sholesky decomposition of `matrix`.
+    /// Attempts to compute the Cholesky decomposition of `matrix`.
     ///
     /// Returns `None` if the input matrix is not definite-positive. The intput matrix is assumed
     /// to be symmetric and only the lower-triangular part is read.
@@ -55,27 +67,29 @@ impl<N: Real, D: DimSub<Dynamic>> Cholesky<N, D>
         Some(Cholesky { chol: matrix })
     }
 
-    /// Retrieves the lower-triangular factor of the cholesky decomposition.
+    /// Retrieves the lower-triangular factor of the Cholesky decomposition with its strictly
+    /// upper-triangular part filled with zeros.
     pub fn unpack(mut self) -> MatrixN<N, D> {
         self.chol.fill_upper_triangle(N::zero(), 1);
         self.chol
     }
 
-    /// Retrieves the lower-triangular factor of che cholesky decomposition, without zeroing-out
+    /// Retrieves the lower-triangular factor of the Cholesky decomposition, without zeroing-out
     /// its strict upper-triangular part.
     ///
-    /// This is an allocation-less version of `self.l()`. The values of the strict upper-triangular
-    /// part are garbage and should be ignored by further computations.
+    /// The values of the strict upper-triangular part are garbage and should be ignored by further
+    /// computations.
     pub fn unpack_dirty(self) -> MatrixN<N, D> {
         self.chol
     }
 
-    /// Retrieves the lower-triangular factor of the cholesky decomposition.
+    /// Retrieves the lower-triangular factor of the Cholesky decomposition with its strictly
+    /// uppen-triangular part filled with zeros.
     pub fn l(&self) -> MatrixN<N, D> {
         self.chol.lower_triangle()
     }
 
-    /// Retrieves the lower-triangular factor of the cholesky decomposition, without zeroing-out
+    /// Retrieves the lower-triangular factor of the Cholesky decomposition, without zeroing-out
     /// its strict upper-triangular part.
     ///
     /// This is an allocation-less version of `self.l()`. The values of the strict upper-triangular
@@ -94,9 +108,8 @@ impl<N: Real, D: DimSub<Dynamic>> Cholesky<N, D>
         self.chol.tr_solve_lower_triangular_mut(b);
     }
 
-    /// Solves the system `self * x = b` where `self` is the decomposed matrix and `x` the unknown.
-    ///
-    /// The result is stored on `b`.
+    /// Returns the solution of the system `self * x = b` where `self` is the decomposed matrix and
+    /// `x` the unknown.
     pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> MatrixMN<N, R2, C2>
         where S2: StorageMut<N, R2, C2>,
               DefaultAllocator: Allocator<N, R2, C2>,
@@ -119,7 +132,7 @@ impl<N: Real, D: DimSub<Dynamic>> Cholesky<N, D>
 impl<N: Real, D: DimSub<Dynamic>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
     where DefaultAllocator: Allocator<N, D, D> {
 
-    /// Attempts to compute the sholesky decomposition of this matrix.
+    /// Attempts to compute the Cholesky decomposition of this matrix.
     ///
     /// Returns `None` if the input matrix is not definite-positive. The intput matrix is assumed
     /// to be symmetric and only the lower-triangular part is read.
diff --git a/src/linalg/eigen.rs b/src/linalg/eigen.rs
index 9d0ec50c..26c2b312 100644
--- a/src/linalg/eigen.rs
+++ b/src/linalg/eigen.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use std::fmt::Display;
 use std::cmp;
 use num_complex::Complex;
@@ -15,7 +18,18 @@ use linalg::RealSchur;
 use geometry::{Reflection, UnitComplex};
 
 
-/// The eigendecomposition of a matrix with real eigenvalues.
+/// Eigendecomposition of a matrix with real eigenvalues.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct RealEigen<N: Real, D: Dim>
     where DefaultAllocator: Allocator<N, D, D> + 
diff --git a/src/linalg/full_piv_lu.rs b/src/linalg/full_piv_lu.rs
index a972efff..bf7f084c 100644
--- a/src/linalg/full_piv_lu.rs
+++ b/src/linalg/full_piv_lu.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use alga::general::Real;
 use core::{Matrix, MatrixN, MatrixMN, DefaultAllocator};
 use dimension::{Dim, DimMin, DimMinimum};
@@ -10,7 +13,20 @@ use linalg::PermutationSequence;
 
 
 
-/// LU decomposition with full pivoting.
+/// LU decomposition with full row and column pivoting.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<(usize, usize), DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Serialize,
+         PermutationSequence<DimMinimum<R, C>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<(usize, usize), DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Deserialize<'de>,
+         PermutationSequence<DimMinimum<R, C>>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct FullPivLU<N: Real, R: DimMin<C>, C: Dim>
     where DefaultAllocator: Allocator<N, R, C> +
@@ -31,7 +47,7 @@ impl<N: Real, R: DimMin<C>, C: Dim> Copy for FullPivLU<N, R, C>
 impl<N: Real, R: DimMin<C>, C: Dim> FullPivLU<N, R, C>
     where DefaultAllocator: Allocator<N, R, C> +
                             Allocator<(usize, usize), DimMinimum<R, C>> {
-    /// Computes the LU decomposition with full-pivoting of `matrix`.
+    /// Computes the LU decomposition with full pivoting of `matrix`.
     ///
     /// This effectively computes `P, L, U, Q` such that `P * matrix * Q = LU`.
     pub fn new(mut matrix: MatrixMN<N, R, C>) -> Self {
@@ -103,13 +119,13 @@ impl<N: Real, R: DimMin<C>, C: Dim> FullPivLU<N, R, C>
         &self.p
     }
 
-    /// The q permutations of this decomposition.
+    /// The column permutations of this decomposition.
     #[inline]
     pub fn q(&self) -> &PermutationSequence<DimMinimum<R, C>> {
         &self.q
     }
 
-    /// The two matrix of this decomposition and the permutation matrix: `(P, L, U, Q)`.
+    /// The two matrices of this decomposition and the row and column permutations: `(P, L, U, Q)`.
     #[inline]
     pub fn unpack(self) -> (PermutationSequence<DimMinimum<R, C>>,
                             MatrixMN<N, R, DimMinimum<R, C>>,
@@ -148,8 +164,8 @@ impl<N: Real, D: DimMin<D, Output = D>> FullPivLU<N, D, D>
 
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
     ///
-    /// If the decomposed matrix is not invertible, this returns `false` and its input `b` is
-    /// left unchanged.
+    /// If the decomposed matrix is not invertible, this returns `false` and its input `b` may
+    /// be overwritten with garbage.
     pub fn solve_mut<R2: Dim, C2: Dim, S2>(&self, b: &mut Matrix<N, R2, C2, S2>) -> bool
         where S2: StorageMut<N, R2, C2>,
               ShapeConstraint: SameNumberOfRows<R2, D> {
@@ -217,7 +233,7 @@ impl<N: Real, D: DimMin<D, Output = D>> FullPivLU<N, D, D>
 impl<N: Real, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
     where DefaultAllocator: Allocator<N, R, C> +
                             Allocator<(usize, usize), DimMinimum<R, C>> {
-    /// Computes the LU decomposition with full-pivoting of `matrix`.
+    /// Computes the LU decomposition with full pivoting of `matrix`.
     ///
     /// This effectively computes `P, L, U, Q` such that `P * matrix * Q = LU`.
     pub fn full_piv_lu(self) -> FullPivLU<N, R, C> {
diff --git a/src/linalg/hessenberg.rs b/src/linalg/hessenberg.rs
index f8514df6..f3f1dc3f 100644
--- a/src/linalg/hessenberg.rs
+++ b/src/linalg/hessenberg.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use alga::general::Real;
 use core::{SquareMatrix, MatrixN, MatrixMN, VectorN, DefaultAllocator};
 use dimension::{DimSub, DimDiff, Dynamic, U1};
@@ -7,7 +10,20 @@ use constraint::{ShapeConstraint, DimEq};
 
 use linalg::householder;
 
-/// The Hessenberg decomposition of a general matrix.
+/// Hessenberg decomposition of a general matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, DimDiff<D, U1>>,
+         MatrixN<N, D>: serde::Serialize,
+         VectorN<N, DimDiff<D, U1>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, DimDiff<D, U1>>,
+         MatrixN<N, D>: serde::Deserialize<'de>,
+         VectorN<N, DimDiff<D, U1>>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct Hessenberg<N: Real, D: DimSub<U1>>
     where DefaultAllocator: Allocator<N, D, D> +
@@ -36,7 +52,8 @@ impl<N: Real, D: DimSub<U1>> Hessenberg<N, D>
 
     /// Computes the Hessenberg decomposition using householder reflections.
     ///
-    /// The workspace containing `D` elements must be provided and may be uninitialized.
+    /// The workspace containing `D` elements must be provided but its content does not have to be
+    /// initialized.
     pub fn new_with_workspace(mut hess: MatrixN<N, D>, work: &mut VectorN<N, D>) -> Self {
         assert!(hess.is_square(), "Cannot compute the hessenberg decomposition of a non-square matrix.");
 
@@ -81,7 +98,7 @@ impl<N: Real, D: DimSub<U1>> Hessenberg<N, D>
     // FIXME: add a h that moves out of self.
     /// Retrieves the upper trapezoidal submatrix `H` of this decomposition.
     ///
-    /// This is less efficient than `.unpack_h()`.
+    /// This is less efficient than `.unpack_h()` as it allocates a new matrix.
     #[inline]
     pub fn h(&self) -> MatrixN<N, D>
         where ShapeConstraint: DimEq<Dynamic, DimDiff<D, U1>> {
diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs
index 50941108..959cff23 100644
--- a/src/linalg/householder.rs
+++ b/src/linalg/householder.rs
@@ -1,4 +1,4 @@
-//! Construction of transformations.
+//! Construction of householder elementary reflections.
 
 use alga::general::Real;
 use core::{Unit, MatrixN, MatrixMN, Vector, VectorN, DefaultAllocator};
diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs
index 6236251a..bdf66935 100644
--- a/src/linalg/lu.rs
+++ b/src/linalg/lu.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use std::mem;
 use alga::general::{Field, Real};
 use core::{Scalar, Matrix, MatrixN, MatrixMN, DefaultAllocator};
@@ -11,6 +14,19 @@ use linalg::PermutationSequence;
 
 
 /// LU decomposition with partial (row) pivoting.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<(usize, usize), DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Serialize,
+         PermutationSequence<DimMinimum<R, C>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<(usize, usize), DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Deserialize<'de>,
+         PermutationSequence<DimMinimum<R, C>>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct LU<N: Real, R: DimMin<C>, C: Dim>
     where DefaultAllocator: Allocator<N, R, C> +
@@ -147,13 +163,13 @@ impl<N: Real, R: DimMin<C>, C: Dim> LU<N, R, C>
         self.lu.rows_generic(0, nrows.min(ncols)).upper_triangle()
     }
 
-    /// The row permutation matrix of this decomposition.
+    /// The row permutations of this decomposition.
     #[inline]
     pub fn p(&self) -> &PermutationSequence<DimMinimum<R, C>> {
         &self.p
     }
 
-    /// The two matrix of this decomposition and the permutation matrix: `(P, L, U)`.
+    /// The row permutations and two triangular matrices of this decomposition: `(P, L, U)`.
     #[inline]
     pub fn unpack(self) -> (PermutationSequence<DimMinimum<R, C>>,
                             MatrixMN<N, R, DimMinimum<R, C>>,
@@ -190,8 +206,8 @@ impl<N: Real, D: DimMin<D, Output = D>> LU<N, D, D>
 
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
     ///
-    /// If the decomposed matrix is not invertible, this returns `false` and its input `b` is
-    /// overwritten with meaningless informations.
+    /// If the decomposed matrix is not invertible, this returns `false` and its input `b` may
+    /// be overwritten with garbage.
     pub fn solve_mut<R2: Dim, C2: Dim, S2>(&self, b: &mut Matrix<N, R2, C2, S2>) -> bool
         where S2: StorageMut<N, R2, C2>,
               ShapeConstraint: SameNumberOfRows<R2, D> {
@@ -206,7 +222,7 @@ impl<N: Real, D: DimMin<D, Output = D>> LU<N, D, D>
 
     /// Computes the inverse of the decomposed matrix.
     ///
-    /// Returnrs `None` if the matrix is not invertible.
+    /// Returns `None` if the matrix is not invertible.
     pub fn try_inverse(&self) -> Option<MatrixN<N, D>> {
         assert!(self.lu.is_square(), "LU inverse: unable to compute the inverse of a non-square matrix.");
 
@@ -222,8 +238,8 @@ impl<N: Real, D: DimMin<D, Output = D>> LU<N, D, D>
 
     /// Computes the inverse of the decomposed matrix and outputs the result to `out`.
     ///
-    /// If the decomposed matrix is not invertible, this returns `false` and `out` may contain
-    /// meaninless informations.
+    /// If the decomposed matrix is not invertible, this returns `false` and `out` may be
+    /// overwritten with garbage.
     pub fn try_inverse_to<S2: StorageMut<N, D, D>>(&self, out: &mut Matrix<N, D, D, S2>) -> bool {
         assert!(self.lu.is_square(), "LU inverse: unable to compute the inverse of a non-square matrix.");
         assert!(self.lu.shape() == out.shape(), "LU inverse: mismatched output shape.");
diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs
index c0b81792..53fed898 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -1,4 +1,4 @@
-//! Linear algebra operators: decompositions, inversion, linear system resolution, etc.
+//! Factorization of real matrices.
 
 mod solve;
 mod determinant;
diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs
index 6d38be31..84b23354 100644
--- a/src/linalg/permutation_sequence.rs
+++ b/src/linalg/permutation_sequence.rs
@@ -1,13 +1,25 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num::One;
 use alga::general::ClosedNeg;
 
 use core::{Scalar, Matrix, VectorN, DefaultAllocator};
-use dimension::{Dim, U1};
+use dimension::{Dim, DimName, Dynamic, U1};
 use storage::StorageMut;
 use allocator::Allocator;
 
 
-/// A sequence of permutations.
+/// A sequence of row or column permutations.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<(usize, usize), D>,
+         VectorN<(usize, usize), D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<(usize, usize), D>,
+         VectorN<(usize, usize), D>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct PermutationSequence<D: Dim>
     where DefaultAllocator: Allocator<(usize, usize), D> {
@@ -19,10 +31,28 @@ impl<D: Dim> Copy for PermutationSequence<D>
     where DefaultAllocator: Allocator<(usize, usize), D>,
           VectorN<(usize, usize), D>: Copy { }
 
-impl<D: Dim> PermutationSequence<D>
+impl<D: DimName> PermutationSequence<D>
     where DefaultAllocator: Allocator<(usize, usize), D> {
 
-    // XXX: Add non-generic constructors.
+    /// Creates a new statically-allocated sequence of `D` identity permutations.
+    #[inline]
+    pub fn identity() -> Self {
+        Self::identity_generic(D::name())
+    }
+}
+
+impl PermutationSequence<Dynamic>
+    where DefaultAllocator: Allocator<(usize, usize), Dynamic> {
+
+    /// Creates a new dynamically-allocated sequence of `n` identity permutations.
+    #[inline]
+    pub fn identity(n: usize) -> Self {
+        Self::identity_generic(Dynamic::new(n))
+    }
+}
+
+impl<D: Dim> PermutationSequence<D>
+    where DefaultAllocator: Allocator<(usize, usize), D> {
     /// Creates a new sequence of D identity permutations.
     #[inline]
     pub fn identity_generic(dim: D) -> Self {
@@ -34,7 +64,8 @@ impl<D: Dim> PermutationSequence<D>
         }
     }
 
-    /// Adds the interchange of the row `i` with the row `i2` to this sequence of permutations.
+    /// Adds the interchange of the row (or column) `i` with the row (or column) `i2` to this
+    /// sequence of permutations.
     #[inline]
     pub fn append_permutation(&mut self, i: usize, i2: usize) {
         if i != i2 {
@@ -44,8 +75,7 @@ impl<D: Dim> PermutationSequence<D>
         }
     }
 
-    /// Permutes the rows of `rhs`, applying a sequence of permutations from the 0-th row to the
-    /// last.
+    /// Applies this sequence of permutations to the rows of `rhs`.
     #[inline]
     pub fn permute_rows<N: Scalar, R2: Dim, C2: Dim, S2>(&self, rhs: &mut Matrix<N, R2, C2, S2>)
         where S2: StorageMut<N, R2, C2> {
@@ -55,7 +85,7 @@ impl<D: Dim> PermutationSequence<D>
         }
     }
 
-    /// Permutes the rows of `rhs` using the inverse permutation matrix of this LUP decomposition.
+    /// Applies this sequence of permutations in reverse to the rows of `rhs`.
     #[inline]
     pub fn inv_permute_rows<N: Scalar, R2: Dim, C2: Dim, S2>(&self, rhs: &mut Matrix<N, R2, C2, S2>)
         where S2: StorageMut<N, R2, C2> {
@@ -66,8 +96,7 @@ impl<D: Dim> PermutationSequence<D>
         }
     }
 
-    /// Permutes the columns of `rhs`, applying a sequence of permutations from the 0-th row to the
-    /// last.
+    /// Applies this sequence of permutations to the columns of `rhs`.
     #[inline]
     pub fn permute_columns<N: Scalar, R2: Dim, C2: Dim, S2>(&self, rhs: &mut Matrix<N, R2, C2, S2>)
         where S2: StorageMut<N, R2, C2> {
@@ -77,7 +106,7 @@ impl<D: Dim> PermutationSequence<D>
         }
     }
 
-    /// Permutes the columns of `rhs` using the inverse permutation matrix of this LUP decomposition.
+    /// Applies this sequence of permutations in reverse to the columns of `rhs`.
     #[inline]
     pub fn inv_permute_columns<N: Scalar, R2: Dim, C2: Dim, S2>(&self, rhs: &mut Matrix<N, R2, C2, S2>)
         where S2: StorageMut<N, R2, C2> {
@@ -88,7 +117,7 @@ impl<D: Dim> PermutationSequence<D>
         }
     }
 
-    /// The number of permutations applied by this matrix.
+    /// The number of non-identity permutations applied by this sequence.
     pub fn len(&self) -> usize {
         self.len
     }
diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs
index 50b4fc19..c9130dfe 100644
--- a/src/linalg/qr.rs
+++ b/src/linalg/qr.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use alga::general::Real;
 use core::{Unit, Matrix, MatrixN, MatrixMN, VectorN, DefaultAllocator};
 use dimension::{Dim, DimMin, DimMinimum, U1};
@@ -10,6 +13,19 @@ use geometry::Reflection;
 
 
 /// The QR decomposition of a general matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<N, DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Serialize,
+         VectorN<N, DimMinimum<R, C>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, R, C> +
+                           Allocator<N, DimMinimum<R, C>>,
+         MatrixMN<N, R, C>: serde::Deserialize<'de>,
+         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct QR<N: Real, R: DimMin<C>, C: Dim>
     where DefaultAllocator: Allocator<N, R, C> +
@@ -150,7 +166,7 @@ impl<N: Real, D: DimMin<D, Output = D>> QR<N, D, D>
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
     ///
     /// If the decomposed matrix is not invertible, this returns `false` and its input `b` is
-    /// overwritten with meaningless informations.
+    /// overwritten with garbage.
     pub fn solve_mut<R2: Dim, C2: Dim, S2>(&self, b: &mut Matrix<N, R2, C2, S2>) -> bool
         where S2: StorageMut<N, R2, C2>,
               ShapeConstraint: SameNumberOfRows<R2, D> {
diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs
index 7b2eb290..15e59211 100644
--- a/src/linalg/schur.rs
+++ b/src/linalg/schur.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use std::cmp;
 use num_complex::Complex;
 use alga::general::Real;
@@ -14,19 +17,26 @@ use geometry::{Reflection, UnitComplex};
 
 
 
-/// Real RealSchur decomposition of a square matrix.
+/// Real Schur decomposition of a square matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D>,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D>,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct RealSchur<N: Real, D: Dim>
-    where DefaultAllocator: Allocator<N, D, D> + 
-                            Allocator<N, D> {
+    where DefaultAllocator: Allocator<N, D, D> {
     q: MatrixN<N, D>,
     t: MatrixN<N, D>
 }
 
 
 impl<N: Real, D: Dim> Copy for RealSchur<N, D>
-    where DefaultAllocator: Allocator<N, D, D> + 
-                            Allocator<N, D>,
+    where DefaultAllocator: Allocator<N, D, D>,
           MatrixN<N, D>: Copy { }
 
 impl<N: Real, D: Dim> RealSchur<N, D>
@@ -37,12 +47,12 @@ impl<N: Real, D: Dim> RealSchur<N, D>
                             Allocator<N, D, D>              + 
                             Allocator<N, D> {
 
-    /// Computes the schur decomposition of a square matrix.
+    /// Computes the Schur decomposition of a square matrix.
     pub fn new(m: MatrixN<N, D>) -> RealSchur<N, D> {
         Self::try_new(m, N::default_epsilon(), 0).unwrap()
     }
 
-    /// Attempts to compute the schur decomposition of a square matrix.
+    /// Attempts to compute the Schur decomposition of a square matrix.
     ///
     /// If only eigenvalues are needed, it is more efficient to call the matrix method
     /// `.eigenvalues()` instead.
@@ -451,12 +461,12 @@ impl<N: Real, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S>
                             Allocator<N, DimDiff<D, U1>>    + // For Hessenberg.
                             Allocator<N, D, D>              + 
                             Allocator<N, D> {
-    /// Computes the schur decomposition of a square matrix.
+    /// Computes the Schur decomposition of a square matrix.
     pub fn real_schur(self) -> RealSchur<N, D> {
         RealSchur::new(self.into_owned())
     }
 
-    /// Attempts to compute the schur decomposition of a square matrix.
+    /// Attempts to compute the Schur decomposition of a square matrix.
     ///
     /// If only eigenvalues are needed, it is more efficient to call the matrix method
     /// `.eigenvalues()` instead.
diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs
index 710d547d..c22ba6af 100644
--- a/src/linalg/svd.rs
+++ b/src/linalg/svd.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num_complex::Complex;
 use std::ops::MulAssign;
 
@@ -15,11 +18,29 @@ use geometry::UnitComplex;
 
 
 
-/// The Singular Value Decomposition of a real matrix.
+/// Singular Value Decomposition of a general matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, R, C>                +
+                           Allocator<N, DimMinimum<R, C>>    +
+                           Allocator<N, DimMinimum<R, C>, C> +
+                           Allocator<N, R, DimMinimum<R, C>>,
+         MatrixMN<N, R, DimMinimum<R, C>>: serde::Serialize,
+         MatrixMN<N, DimMinimum<R, C>, C>: serde::Serialize,
+         VectorN<N, DimMinimum<R, C>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, R, C>                +
+                           Allocator<N, DimMinimum<R, C>>    +
+                           Allocator<N, DimMinimum<R, C>, C> +
+                           Allocator<N, R, DimMinimum<R, C>>,
+         MatrixMN<N, R, DimMinimum<R, C>>: serde::Deserialize<'de>,
+         MatrixMN<N, DimMinimum<R, C>, C>: serde::Deserialize<'de>,
+         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct SVD<N: Real, R: DimMin<C>, C: Dim>
-    where DefaultAllocator: Allocator<N, R, C>                +
-                            Allocator<N, DimMinimum<R, C>, C> +
+    where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C> +
                             Allocator<N, R, DimMinimum<R, C>> +
                             Allocator<N, DimMinimum<R, C>> {
     /// The left-singular vectors `U` of this SVD.
@@ -31,9 +52,8 @@ pub struct SVD<N: Real, R: DimMin<C>, C: Dim>
 }
 
 
-impl<N: Real, R: DimMin<C>, C: Dim> SVD<N, R, C>
-    where DefaultAllocator: Allocator<N, R, C>                +
-                            Allocator<N, DimMinimum<R, C>, C> +
+impl<N: Real, R: DimMin<C>, C: Dim> Copy for SVD<N, R, C>
+    where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C> +
                             Allocator<N, R, DimMinimum<R, C>> +
                             Allocator<N, DimMinimum<R, C>>,
           MatrixMN<N, R, DimMinimum<R, C>>: Copy,
diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs
index b66154dd..1d417dd1 100644
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use num_complex::Complex;
 use std::ops::MulAssign;
 
@@ -12,7 +15,20 @@ use linalg::SymmetricTridiagonal;
 use geometry::UnitComplex;
 
 
-/// The eigendecomposition of a symmetric matrix.
+/// Eigendecomposition of a symmetric matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, D>,
+         VectorN<N, D>: serde::Serialize,
+         MatrixN<N, D>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, D>,
+         VectorN<N, D>: serde::Deserialize<'de>,
+         MatrixN<N, D>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct SymmetricEigen<N: Real, D: Dim>
     where DefaultAllocator: Allocator<N, D, D> +
@@ -24,7 +40,7 @@ pub struct SymmetricEigen<N: Real, D: Dim>
     pub eigenvalues: VectorN<N, D>
 }
 
-impl<N: Real, D: Dim> SymmetricEigen<N, D>
+impl<N: Real, D: Dim> Copy for SymmetricEigen<N, D>
     where DefaultAllocator: Allocator<N, D, D> +
                             Allocator<N, D>,
           MatrixN<N, D>: Copy,
@@ -35,7 +51,7 @@ impl<N: Real, D: Dim> SymmetricEigen<N, D>
                             Allocator<N, D> {
     /// Computes the eigendecomposition of the given symmetric matrix.
     ///
-    /// Only the lower-triangular and diagonal parts of `m` are read.
+    /// Only the lower-triangular parts (including its diagonal) of `m` is read.
     pub fn new(m: MatrixN<N, D>) -> Self
         where D: DimSub<U1>,
               DefaultAllocator: Allocator<N, DimDiff<D, U1>> {
@@ -46,11 +62,11 @@ impl<N: Real, D: Dim> SymmetricEigen<N, D>
     /// Computes the eigendecomposition of the given symmetric matrix with user-specified
     /// convergence parameters.
     ///
-    /// Only the lower-triangular and diagonal parts of `m` are read.
+    /// Only the lower-triangular part (including its diagonal) of `m` is read.
     ///
     /// # Arguments
     ///
-    /// * `eps`       − tolerence used to determine when a value converged to 0.
+    /// * `eps`       − tolerance used to determine when a value converged to 0.
     /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
     /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
     /// continues indefinitely until convergence.
@@ -277,7 +293,7 @@ impl<N: Real, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
 
     /// Computes the eigendecomposition of this symmetric matrix.
     ///
-    /// Only the lower-triangular part (including the diagonal) of `m` are read.
+    /// Only the lower-triangular part (including the diagonal) of `m` is read.
     pub fn symmetric_eigen(self) -> SymmetricEigen<N, D> {
         SymmetricEigen::new(self.into_owned())
     }
@@ -285,11 +301,11 @@ impl<N: Real, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
     /// Computes the eigendecomposition of the given symmetric matrix with user-specified
     /// convergence parameters.
     ///
-    /// Only the lower-triangular and diagonal parts of `m` are read.
+    /// Only the lower-triangular part (including the diagonal) of `m` is read.
     ///
     /// # Arguments
     ///
-    /// * `eps`       − tolerence used to determine when a value converged to 0.
+    /// * `eps`       − tolerance used to determine when a value converged to 0.
     /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
     /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
     /// continues indefinitely until convergence.
diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs
index bb922928..63dd35d8 100644
--- a/src/linalg/symmetric_tridiagonal.rs
+++ b/src/linalg/symmetric_tridiagonal.rs
@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde;
+
 use alga::general::Real;
 use core::{SquareMatrix, MatrixN, MatrixMN, VectorN, DefaultAllocator};
 use dimension::{DimSub, DimDiff, U1};
@@ -7,7 +10,20 @@ use allocator::Allocator;
 use linalg::householder;
 
 
-/// The tridiagonalization of a symmetric matrix.
+/// Tridiagonalization of a symmetric matrix.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(serialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, DimDiff<D, U1>>,
+         MatrixN<N, D>: serde::Serialize,
+         VectorN<N, DimDiff<D, U1>>: serde::Serialize")))]
+#[cfg_attr(feature = "serde-serialize",
+    serde(bound(deserialize =
+        "DefaultAllocator: Allocator<N, D, D> +
+                           Allocator<N, DimDiff<D, U1>>,
+         MatrixN<N, D>: serde::Deserialize<'de>,
+         VectorN<N, DimDiff<D, U1>>: serde::Deserialize<'de>")))]
 #[derive(Clone, Debug)]
 pub struct SymmetricTridiagonal<N: Real, D: DimSub<U1>>
     where DefaultAllocator: Allocator<N, D, D> +
@@ -16,7 +32,7 @@ pub struct SymmetricTridiagonal<N: Real, D: DimSub<U1>>
     off_diagonal: VectorN<N, DimDiff<D, U1>>
 }
 
-impl<N: Real, D: DimSub<U1>> SymmetricTridiagonal<N, D>
+impl<N: Real, D: DimSub<U1>> Copy for SymmetricTridiagonal<N, D>
     where DefaultAllocator: Allocator<N, D, D> +
                             Allocator<N, DimDiff<D, U1>>,
           MatrixN<N, D>: Copy,
@@ -28,7 +44,7 @@ impl<N: Real, D: DimSub<U1>> SymmetricTridiagonal<N, D>
 
     /// Computes the tridiagonalization of the symmetric matrix `m`.
     ///
-    /// Only the lower-triangular and diagonal parts of `m` are read.
+    /// Only the lower-triangular part (including the diagonal) of `m` is read.
     pub fn new(mut m: MatrixN<N, D>) -> Self {
         let dim = m.data.shape().0;
 
@@ -124,7 +140,7 @@ impl<N: Real, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
 
     /// Computes the tridiagonalization of this symmetric matrix.
     ///
-    /// Only the lower-triangular and diagonal parts of `self` are read.
+    /// Only the lower-triangular part (including the diagonal) of `m` is read.
     pub fn symmetric_tridiagonalize(self) -> SymmetricTridiagonal<N, D> {
         SymmetricTridiagonal::new(self.into_owned())
     }