From e6265da98033c6ebf4e939737b8de2cd4d38cbd6 Mon Sep 17 00:00:00 2001
From: Mara Bos <m-ou.se@m-ou.se>
Date: Tue, 16 Apr 2019 10:11:27 +0200
Subject: [PATCH 01/20] Implement Default for MatrixMN.

---
 src/base/array_storage.rs | 13 +++++++++++++
 src/base/matrix.rs        | 15 +++++++++++++++
 2 files changed, 28 insertions(+)

diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs
index bebb8740..450127ba 100644
--- a/src/base/array_storage.rs
+++ b/src/base/array_storage.rs
@@ -48,6 +48,19 @@ where
 /// Renamed to [ArrayStorage].
 pub type MatrixArray<N, R, C> = ArrayStorage<N, R, C>;
 
+impl<N, R, C> Default for ArrayStorage<N, R, C>
+where
+    R: DimName,
+    C: DimName,
+    R::Value: Mul<C::Value>,
+    Prod<R::Value, C::Value>: ArrayLength<N>,
+    N: Default,
+{
+    fn default() -> Self {
+        ArrayStorage { data: Default::default() }
+    }
+}
+
 impl<N, R, C> Hash for ArrayStorage<N, R, C>
 where
     N: Hash,
diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index 7119a8ad..1aa171e4 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -90,6 +90,21 @@ impl<N: Scalar, R: Dim, C: Dim, S: fmt::Debug> fmt::Debug for Matrix<N, R, C, S>
     }
 }
 
+impl<N, R, C, S> Default for Matrix<N, R, C, S>
+where
+    N: Scalar,
+    R: Dim,
+    C: Dim,
+    S: Default,
+{
+    fn default() -> Self {
+        Matrix {
+            data: Default::default(),
+            _phantoms: PhantomData,
+        }
+    }
+}
+
 #[cfg(feature = "serde-serialize")]
 impl<N, R, C, S> Serialize for Matrix<N, R, C, S>
 where

From 85c931520abc0fd0d416327062b35fcb5ba56320 Mon Sep 17 00:00:00 2001
From: Mara Bos <m-ou.se@m-ou.se>
Date: Tue, 16 Apr 2019 10:12:09 +0200
Subject: [PATCH 02/20] Implement Default for Quaternion and UnitQuaternion.

---
 src/geometry/quaternion.rs | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs
index 40c085f2..5b126983 100755
--- a/src/geometry/quaternion.rs
+++ b/src/geometry/quaternion.rs
@@ -30,6 +30,14 @@ pub struct Quaternion<N: RealField> {
     pub coords: Vector4<N>,
 }
 
+impl<N: RealField> Default for Quaternion<N> {
+    fn default() -> Self {
+        Quaternion {
+            coords: Vector4::zeros()
+        }
+    }
+}
+
 #[cfg(feature = "abomonation-serialize")]
 impl<N: RealField> Abomonation for Quaternion<N>
 where Vector4<N>: Abomonation
@@ -1422,6 +1430,12 @@ impl<N: RealField> UnitQuaternion<N> {
     }
 }
 
+impl<N: RealField> Default for UnitQuaternion<N> {
+    fn default() -> Self {
+        Self::identity()
+    }
+}
+
 impl<N: RealField + fmt::Display> fmt::Display for UnitQuaternion<N> {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         if let Some(axis) = self.axis() {

From c0a6df55b17b767340e2bd65e98df136b356ff2c Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Wed, 1 Apr 2020 22:36:05 +0200
Subject: [PATCH 03/20] Addition of matrix exponent for static size matrices.

---
 src/linalg/exp.rs   | 445 ++++++++++++++++++++++++++++++++++++++++++++
 src/linalg/mod.rs   |   1 +
 tests/linalg/exp.rs |  97 ++++++++++
 tests/linalg/mod.rs |   2 +
 4 files changed, 545 insertions(+)
 create mode 100644 src/linalg/exp.rs
 create mode 100644 tests/linalg/exp.rs

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
new file mode 100644
index 00000000..8e8227d5
--- /dev/null
+++ b/src/linalg/exp.rs
@@ -0,0 +1,445 @@
+use crate::{
+    base::{
+        allocator::Allocator,
+        dimension::{DimMin, DimMinimum, DimName},
+        DefaultAllocator,
+    },
+    try_convert, ComplexField, MatrixN, RealField,
+};
+
+// https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/matfuncs.py
+struct ExpmPadeHelper<N, R>
+where
+    N: RealField,
+    R: DimName + DimMin<R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    use_exact_norm: bool,
+    ident: MatrixN<N, R>,
+
+    a: MatrixN<N, R>,
+    a2: Option<MatrixN<N, R>>,
+    a4: Option<MatrixN<N, R>>,
+    a6: Option<MatrixN<N, R>>,
+    a8: Option<MatrixN<N, R>>,
+    a10: Option<MatrixN<N, R>>,
+
+    d4_exact: Option<N>,
+    d6_exact: Option<N>,
+    d8_exact: Option<N>,
+    d10_exact: Option<N>,
+
+    d4_approx: Option<N>,
+    d6_approx: Option<N>,
+    d8_approx: Option<N>,
+    d10_approx: Option<N>,
+}
+
+impl<N, R> ExpmPadeHelper<N, R>
+where
+    N: RealField,
+    R: DimName + DimMin<R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    fn new(a: MatrixN<N, R>, use_exact_norm: bool) -> Self {
+        ExpmPadeHelper {
+            use_exact_norm,
+            ident: MatrixN::<N, R>::identity(),
+            a,
+            a2: None,
+            a4: None,
+            a6: None,
+            a8: None,
+            a10: None,
+            d4_exact: None,
+            d6_exact: None,
+            d8_exact: None,
+            d10_exact: None,
+            d4_approx: None,
+            d6_approx: None,
+            d8_approx: None,
+            d10_approx: None,
+        }
+    }
+
+    fn a2(&self) -> &MatrixN<N, R> {
+        if self.a2.is_none() {
+            let ap = &self.a2 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(&self.a * &self.a);
+            }
+        }
+        self.a2.as_ref().unwrap()
+    }
+
+    fn a4(&self) -> &MatrixN<N, R> {
+        if self.a4.is_none() {
+            let ap = &self.a4 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let a2 = self.a2();
+            unsafe {
+                *ap = Some(a2 * a2);
+            }
+        }
+        self.a4.as_ref().unwrap()
+    }
+
+    fn a6(&self) -> &MatrixN<N, R> {
+        if self.a6.is_none() {
+            let a2 = self.a2();
+            let a4 = self.a4();
+            let ap = &self.a6 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(a4 * a2);
+            }
+        }
+        self.a6.as_ref().unwrap()
+    }
+
+    fn a8(&self) -> &MatrixN<N, R> {
+        if self.a8.is_none() {
+            let a2 = self.a2();
+            let a6 = self.a6();
+            let ap = &self.a8 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(a6 * a2);
+            }
+        }
+        self.a8.as_ref().unwrap()
+    }
+
+    fn a10(&mut self) -> &MatrixN<N, R> {
+        if self.a10.is_none() {
+            let a4 = self.a4();
+            let a6 = self.a6();
+            let ap = &self.a10 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(a6 * a4);
+            }
+        }
+        self.a10.as_ref().unwrap()
+    }
+
+    fn d4_tight(&mut self) -> N {
+        if self.d4_exact.is_none() {
+            self.d4_exact = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+        }
+        self.d4_exact.unwrap()
+    }
+
+    fn d6_tight(&mut self) -> N {
+        if self.d6_exact.is_none() {
+            self.d6_exact = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+        }
+        self.d6_exact.unwrap()
+    }
+
+    fn d8_tight(&mut self) -> N {
+        if self.d8_exact.is_none() {
+            self.d8_exact = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+        }
+        self.d8_exact.unwrap()
+    }
+
+    fn d10_tight(&mut self) -> N {
+        if self.d10_exact.is_none() {
+            self.d10_exact = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+        }
+        self.d10_exact.unwrap()
+    }
+
+    fn d4_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d4_tight();
+        }
+
+        if self.d4_exact.is_some() {
+            return self.d4_exact.unwrap();
+        }
+
+        if self.d4_approx.is_none() {
+            self.d4_approx = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+        }
+
+        self.d4_approx.unwrap()
+    }
+
+    fn d6_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d6_tight();
+        }
+
+        if self.d6_exact.is_some() {
+            return self.d6_exact.unwrap();
+        }
+
+        if self.d6_approx.is_none() {
+            self.d6_approx = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+        }
+
+        self.d6_approx.unwrap()
+    }
+
+    fn d8_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d8_tight();
+        }
+
+        if self.d8_exact.is_some() {
+            return self.d8_exact.unwrap();
+        }
+
+        if self.d8_approx.is_none() {
+            self.d8_approx = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+        }
+
+        self.d8_approx.unwrap()
+    }
+
+    fn d10_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d10_tight();
+        }
+
+        if self.d10_exact.is_some() {
+            return self.d10_exact.unwrap();
+        }
+
+        if self.d10_approx.is_none() {
+            self.d10_approx = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+        }
+
+        self.d10_approx.unwrap()
+    }
+
+    fn pade3(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(120.0).unwrap(),
+            N::from_f64(60.0).unwrap(),
+            N::from_f64(12.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u = &self.a * (self.a2() * b[3] + &self.ident * b[1]);
+        let v = self.a2() * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade5(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(30240.0).unwrap(),
+            N::from_f64(15120.0).unwrap(),
+            N::from_f64(3360.0).unwrap(),
+            N::from_f64(420.0).unwrap(),
+            N::from_f64(30.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u = &self.a * (self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
+        let v = self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade7(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(17297280.0).unwrap(),
+            N::from_f64(8648640.0).unwrap(),
+            N::from_f64(1995840.0).unwrap(),
+            N::from_f64(277200.0).unwrap(),
+            N::from_f64(25200.0).unwrap(),
+            N::from_f64(1512.0).unwrap(),
+            N::from_f64(56.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u =
+            &self.a * (self.a6() * b[7] + self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
+        let v = self.a6() * b[6] + self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade9(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(17643225600.0).unwrap(),
+            N::from_f64(8821612800.0).unwrap(),
+            N::from_f64(2075673600.0).unwrap(),
+            N::from_f64(302702400.0).unwrap(),
+            N::from_f64(30270240.0).unwrap(),
+            N::from_f64(2162160.0).unwrap(),
+            N::from_f64(110880.0).unwrap(),
+            N::from_f64(3960.0).unwrap(),
+            N::from_f64(90.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u = &self.a
+            * (self.a8() * b[9]
+                + self.a6() * b[7]
+                + self.a4() * b[5]
+                + self.a2() * b[3]
+                + &self.ident * b[1]);
+        let v = self.a8() * b[8]
+            + self.a6() * b[6]
+            + self.a4() * b[4]
+            + self.a2() * b[2]
+            + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade13_scaled(&mut self, s: u64) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(64764752532480000.0).unwrap(),
+            N::from_f64(32382376266240000.0).unwrap(),
+            N::from_f64(7771770303897600.0).unwrap(),
+            N::from_f64(1187353796428800.0).unwrap(),
+            N::from_f64(129060195264000.0).unwrap(),
+            N::from_f64(10559470521600.0).unwrap(),
+            N::from_f64(670442572800.0).unwrap(),
+            N::from_f64(33522128640.0).unwrap(),
+            N::from_f64(1323241920.0).unwrap(),
+            N::from_f64(40840800.0).unwrap(),
+            N::from_f64(960960.0).unwrap(),
+            N::from_f64(16380.0).unwrap(),
+            N::from_f64(182.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let s = s as f64;
+
+        let mb = &self.a * N::from_f64(2.0.powf(-s)).unwrap();
+        let mb2 = self.a2() * N::from_f64(2.0.powf(-2.0 * s)).unwrap();
+        let mb4 = self.a4() * N::from_f64(2.0.powf(-4.0 * s)).unwrap();
+        let mb6 = self.a6() * N::from_f64(2.0.powf(-6.0 * s)).unwrap();
+
+        let u2 = &mb6 * (&mb6 * b[13] + &mb4 * b[11] + &mb2 * b[9]);
+        let u = &mb * (&u2 + &mb6 * b[7] + &mb4 * b[5] + &mb2 * b[3] + &self.ident * b[1]);
+        let v2 = &mb6 * (&mb6 * b[12] + &mb4 * b[10] + &mb2 * b[8]);
+        let v = v2 + &mb6 * b[6] + &mb4 * b[4] + &mb2 * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+}
+
+fn factorial(n: u128) -> u128 {
+    if n == 1 {
+        return 1;
+    }
+    n * factorial(n - 1)
+}
+
+fn onenorm_matrix_power_nnm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
+where
+    N: RealField,
+    R: DimName,
+    DefaultAllocator: Allocator<N, R, R>,
+{
+    let mut v = MatrixN::<N, R>::repeat(N::from_f64(1.0).unwrap());
+    let m = a.transpose();
+
+    for _ in 0..p {
+        v = &m * v;
+    }
+
+    v.amax()
+}
+
+fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
+where
+    N: RealField,
+    R: DimName,
+    DefaultAllocator: Allocator<N, R, R>,
+{
+    // 2m choose m = (2m)!/(m! * (2m-m)!)
+
+    let a_abs_onenorm = onenorm_matrix_power_nnm(&a.abs(), 2 * m + 1);
+
+    if a_abs_onenorm == N::zero() {
+        return 0;
+    }
+
+    let choose_2m_m =
+        factorial(2 * m as u128) / (factorial(m as u128) * factorial(2 * m as u128 - m as u128));
+    let abs_c_recip = choose_2m_m * factorial(2 * m as u128 + 1);
+    let alpha = a_abs_onenorm / a.amax();
+    let alpha = alpha / N::from_u128(abs_c_recip).unwrap();
+
+    let u = N::from_f64(2_f64.powf(-53.0)).unwrap();
+    let log2_alpha_div_u = try_convert((alpha / u).log2()).unwrap();
+    let value = (log2_alpha_div_u / (2.0 * m as f64)).ceil();
+    if value > 0.0 {
+        value as u64
+    } else {
+        0
+    }
+}
+
+fn solve_p_q<N, R>(u: MatrixN<N, R>, v: MatrixN<N, R>) -> MatrixN<N, R>
+where
+    N: ComplexField,
+    R: DimMin<R, Output = R> + DimName,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    let p = &u + &v;
+    let q = &v - &u;
+
+    q.lu().solve(&p).unwrap()
+}
+
+impl<N: RealField, R: DimMin<R, Output = R> + DimName> MatrixN<N, R>
+where
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    /// Computes exp of this matrix
+    pub fn exp(&self) -> Self {
+        // Simple case
+        if self.nrows() == 1 {
+            return self.clone().map(|v| v.exp());
+        }
+
+        let mut h = ExpmPadeHelper::new(self.clone(), true);
+
+        let eta_1 = N::max(h.d4_loose(), h.d6_loose());
+        if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
+            let (u, v) = h.pade3();
+            return solve_p_q(u, v);
+        }
+
+        let eta_2 = N::max(h.d4_tight(), h.d6_loose());
+        if eta_2 < N::from_f64(2.539398330063230e-001).unwrap() && ell(&h.a, 5) == 0 {
+            let (u, v) = h.pade5();
+            return solve_p_q(u, v);
+        }
+
+        let eta_3 = N::max(h.d6_tight(), h.d8_loose());
+        if eta_3 < N::from_f64(9.504178996162932e-001).unwrap() && ell(&h.a, 7) == 0 {
+            let (u, v) = h.pade7();
+            return solve_p_q(u, v);
+        }
+        if eta_3 < N::from_f64(2.097847961257068e+000).unwrap() && ell(&h.a, 9) == 0 {
+            let (u, v) = h.pade9();
+            return solve_p_q(u, v);
+        }
+
+        let eta_4 = N::max(h.d8_loose(), h.d10_loose());
+        let eta_5 = N::min(eta_3, eta_4);
+        let theta_13 = N::from_f64(4.25).unwrap();
+
+        let mut s = if eta_5 == N::zero() {
+            0
+        } else {
+            let l2 = try_convert((eta_5 / theta_13).log2().ceil()).unwrap();
+
+            if l2 < 0.0 {
+                0
+            } else {
+                l2 as u64
+            }
+        };
+
+        s += ell(
+            &(&h.a * N::from_f64(2.0_f64.powf(-(s as f64))).unwrap()),
+            13,
+        );
+
+        let (u, v) = h.pade13_scaled(s);
+        let mut x = solve_p_q(u, v);
+
+        for _ in 0..s {
+            x = &x * &x;
+        }
+        x
+    }
+}
diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs
index de1108f7..42bfde63 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -5,6 +5,7 @@ mod bidiagonal;
 mod cholesky;
 mod convolution;
 mod determinant;
+mod exp;
 mod full_piv_lu;
 pub mod givens;
 mod hessenberg;
diff --git a/tests/linalg/exp.rs b/tests/linalg/exp.rs
new file mode 100644
index 00000000..ac3f8dd7
--- /dev/null
+++ b/tests/linalg/exp.rs
@@ -0,0 +1,97 @@
+#[cfg(test)]
+mod tests {
+    //https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/tests/test_matfuncs.py
+    #[test]
+    fn exp() {
+        use nalgebra::{Matrix1, Matrix2, Matrix3};
+
+        {
+            let m = Matrix1::new(1.0);
+
+            let f = m.exp();
+
+            assert_eq!(f, Matrix1::new(1_f64.exp()));
+        }
+
+        {
+            let m = Matrix2::new(0.0, 1.0, 0.0, 0.0);
+
+            assert_eq!(m.exp(), Matrix2::new(1.0, 1.0, 0.0, 1.0));
+        }
+
+        {
+            let a: f64 = 1.0;
+            let b: f64 = 2.0;
+            let c: f64 = 3.0;
+            let d: f64 = 4.0;
+            let m = Matrix2::new(a, b, c, d);
+
+            let delta = ((a - d).powf(2.0) + 4.0 * b * c).sqrt();
+            let delta_2 = delta / 2.0;
+            let ad_2 = (a + d) / 2.0;
+            let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh());
+            let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh();
+            let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh();
+            let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
+
+            let f = Matrix2::new(m11, m12, m21, m22) / delta;
+            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+        }
+
+        {
+            // https://mathworld.wolfram.com/MatrixExponential.html
+            use rand::{
+                distributions::{Distribution, Uniform},
+                thread_rng,
+            };
+            let mut rng = thread_rng();
+            let dist = Uniform::new(-10.0, 10.0);
+            loop {
+                let a: f64 = dist.sample(&mut rng);
+                let b: f64 = dist.sample(&mut rng);
+                let c: f64 = dist.sample(&mut rng);
+                let d: f64 = dist.sample(&mut rng);
+                let m = Matrix2::new(a, b, c, d);
+
+                let delta_sq = (a - d).powf(2.0) + 4.0 * b * c;
+                if delta_sq < 0.0 {
+                    continue;
+                }
+
+                let delta = delta_sq.sqrt();
+                let delta_2 = delta / 2.0;
+                let ad_2 = (a + d) / 2.0;
+                let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh());
+                let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh();
+                let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh();
+                let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
+
+                let f = Matrix2::new(m11, m12, m21, m22) / delta;
+                println!("a: {}", m);
+                assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+                break;
+            }
+        }
+
+        {
+            let m = Matrix3::new(1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0);
+
+            let e1 = 1.0_f64.exp();
+            let e2 = 2.0_f64.exp();
+
+            let f = Matrix3::new(
+                e1,
+                3.0 * e1,
+                15.0 * (e2 - 2.0 * e1),
+                0.0,
+                e1,
+                5.0 * (e2 - e1),
+                0.0,
+                0.0,
+                e2,
+            );
+
+            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+        }
+    }
+}
diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs
index 234cac39..c810df1f 100644
--- a/tests/linalg/mod.rs
+++ b/tests/linalg/mod.rs
@@ -1,7 +1,9 @@
 mod balancing;
 mod bidiagonal;
 mod cholesky;
+mod convolution;
 mod eigen;
+mod exp;
 mod full_piv_lu;
 mod hessenberg;
 mod inverse;

From 2696ffd7a1837b73756b9512b32781532b4c66ca Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Thu, 2 Apr 2020 09:38:18 +0200
Subject: [PATCH 04/20] Fixed bug by introducing one norm

---
 src/linalg/exp.rs | 51 +++++++++++++++++++++++++++++++++++++----------
 1 file changed, 40 insertions(+), 11 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 8e8227d5..877ce0d3 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -121,28 +121,28 @@ where
 
     fn d4_tight(&mut self) -> N {
         if self.d4_exact.is_none() {
-            self.d4_exact = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+            self.d4_exact = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
         }
         self.d4_exact.unwrap()
     }
 
     fn d6_tight(&mut self) -> N {
         if self.d6_exact.is_none() {
-            self.d6_exact = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_exact = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
         }
         self.d6_exact.unwrap()
     }
 
     fn d8_tight(&mut self) -> N {
         if self.d8_exact.is_none() {
-            self.d8_exact = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_exact = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
         }
         self.d8_exact.unwrap()
     }
 
     fn d10_tight(&mut self) -> N {
         if self.d10_exact.is_none() {
-            self.d10_exact = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_exact = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
         }
         self.d10_exact.unwrap()
     }
@@ -157,7 +157,7 @@ where
         }
 
         if self.d4_approx.is_none() {
-            self.d4_approx = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+            self.d4_approx = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
         }
 
         self.d4_approx.unwrap()
@@ -173,7 +173,7 @@ where
         }
 
         if self.d6_approx.is_none() {
-            self.d6_approx = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_approx = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
         }
 
         self.d6_approx.unwrap()
@@ -189,7 +189,7 @@ where
         }
 
         if self.d8_approx.is_none() {
-            self.d8_approx = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_approx = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
         }
 
         self.d8_approx.unwrap()
@@ -205,7 +205,7 @@ where
         }
 
         if self.d10_approx.is_none() {
-            self.d10_approx = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_approx = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
         }
 
         self.d10_approx.unwrap()
@@ -333,7 +333,7 @@ where
         v = &m * v;
     }
 
-    v.amax()
+    one_norm(&v)
 }
 
 fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
@@ -353,7 +353,7 @@ where
     let choose_2m_m =
         factorial(2 * m as u128) / (factorial(m as u128) * factorial(2 * m as u128 - m as u128));
     let abs_c_recip = choose_2m_m * factorial(2 * m as u128 + 1);
-    let alpha = a_abs_onenorm / a.amax();
+    let alpha = a_abs_onenorm / one_norm(a);
     let alpha = alpha / N::from_u128(abs_c_recip).unwrap();
 
     let u = N::from_f64(2_f64.powf(-53.0)).unwrap();
@@ -378,6 +378,24 @@ where
     q.lu().solve(&p).unwrap()
 }
 
+pub fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
+where
+    N: RealField,
+    R: DimName,
+    DefaultAllocator: Allocator<N, R, R>,
+{
+    let mut col_sums = vec![N::zero(); m.ncols()];
+    for i in 0..m.ncols() {
+        let col = m.column(i);
+        col.iter().for_each(|v| col_sums[i] += v.abs());
+    }
+    let mut max = col_sums[0];
+    for i in 1..col_sums.len() {
+        max = N::max(max, col_sums[i]);
+    }
+    max
+}
+
 impl<N: RealField, R: DimMin<R, Output = R> + DimName> MatrixN<N, R>
 where
     DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
@@ -389,7 +407,7 @@ where
             return self.clone().map(|v| v.exp());
         }
 
-        let mut h = ExpmPadeHelper::new(self.clone(), true);
+        let mut h = ExpmPadeHelper::new(self.clone(), false);
 
         let eta_1 = N::max(h.d4_loose(), h.d6_loose());
         if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
@@ -443,3 +461,14 @@ where
         x
     }
 }
+
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn one_norm() {
+        use crate::Matrix3;
+        let m = Matrix3::new(-3.0, 5.0, 7.0, 2.0, 6.0, 4.0, 0.0, 2.0, 8.0);
+
+        assert_eq!(super::one_norm(&m), 19.0);
+    }
+}

From e1565616778f4635088d47f88bb46c1ff8bc4c4d Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Thu, 2 Apr 2020 09:39:50 +0200
Subject: [PATCH 05/20] one_norm is not a public function

---
 src/linalg/exp.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 877ce0d3..5b76dcaf 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -378,7 +378,7 @@ where
     q.lu().solve(&p).unwrap()
 }
 
-pub fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
+fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
 where
     N: RealField,
     R: DimName,

From b3ef66fd14b17ffb54f8cb560a5866c81a6aff4a Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Thu, 2 Apr 2020 10:55:00 +0200
Subject: [PATCH 06/20] Fixed bug in onenorm_matrix_power_norm

---
 src/linalg/exp.rs | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 5b76dcaf..7c087159 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -320,31 +320,32 @@ fn factorial(n: u128) -> u128 {
     n * factorial(n - 1)
 }
 
-fn onenorm_matrix_power_nnm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
+/// Compute the 1-norm of a non-negative integer power of a non-negative matrix.
+fn onenorm_matrix_power_nonm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
 where
     N: RealField,
     R: DimName,
-    DefaultAllocator: Allocator<N, R, R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
 {
-    let mut v = MatrixN::<N, R>::repeat(N::from_f64(1.0).unwrap());
+    let mut v = crate::VectorN::<N, R>::repeat(N::from_f64(1.0).unwrap());
     let m = a.transpose();
 
     for _ in 0..p {
         v = &m * v;
     }
 
-    one_norm(&v)
+    v.max()
 }
 
 fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
 where
     N: RealField,
     R: DimName,
-    DefaultAllocator: Allocator<N, R, R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
 {
     // 2m choose m = (2m)!/(m! * (2m-m)!)
 
-    let a_abs_onenorm = onenorm_matrix_power_nnm(&a.abs(), 2 * m + 1);
+    let a_abs_onenorm = onenorm_matrix_power_nonm(&a.abs(), 2 * m + 1);
 
     if a_abs_onenorm == N::zero() {
         return 0;
@@ -396,9 +397,11 @@ where
     max
 }
 
-impl<N: RealField, R: DimMin<R, Output = R> + DimName> MatrixN<N, R>
+impl<N: RealField, R> MatrixN<N, R>
 where
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    R: DimMin<R, Output = R> + DimName,
+    DefaultAllocator:
+        Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>> + Allocator<N, R>,
 {
     /// Computes exp of this matrix
     pub fn exp(&self) -> Self {

From 4ec94408b5652e12b779df462479e15df8bc7f33 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Fri, 3 Apr 2020 11:11:14 +0200
Subject: [PATCH 07/20] Pub use of exp in the linalg module

---
 src/linalg/exp.rs | 6 ++++--
 src/linalg/mod.rs | 1 +
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 7c087159..8822027b 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -1,3 +1,5 @@
+//! This module provides the matrix exponent (exp) function to square matrices.
+//!
 use crate::{
     base::{
         allocator::Allocator,
@@ -403,14 +405,14 @@ where
     DefaultAllocator:
         Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>> + Allocator<N, R>,
 {
-    /// Computes exp of this matrix
+    /// Computes exponential of this matrix
     pub fn exp(&self) -> Self {
         // Simple case
         if self.nrows() == 1 {
             return self.clone().map(|v| v.exp());
         }
 
-        let mut h = ExpmPadeHelper::new(self.clone(), false);
+        let mut h = ExpmPadeHelper::new(self.clone(), true);
 
         let eta_1 = N::max(h.d4_loose(), h.d6_loose());
         if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs
index 42bfde63..f96cef0c 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -27,6 +27,7 @@ mod symmetric_tridiagonal;
 pub use self::bidiagonal::*;
 pub use self::cholesky::*;
 pub use self::convolution::*;
+pub use self::exp::*;
 pub use self::full_piv_lu::*;
 pub use self::hessenberg::*;
 pub use self::lu::*;

From dbbf87a3dd475c454a5f69cfb355f9fd6b4da5c3 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Tue, 7 Apr 2020 09:44:06 +0200
Subject: [PATCH 08/20] Rebased against dev

---
 tests/linalg/mod.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs
index c810df1f..7fc01396 100644
--- a/tests/linalg/mod.rs
+++ b/tests/linalg/mod.rs
@@ -12,5 +12,4 @@ mod qr;
 mod schur;
 mod solve;
 mod svd;
-mod convolution;
 mod tridiagonal;

From c7d9e415ce3531b2ffb8399c2c4e3682c54068bf Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Tue, 7 Apr 2020 09:55:58 +0200
Subject: [PATCH 09/20] Converted tests to use relative_eq

---
 tests/linalg/exp.rs | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/tests/linalg/exp.rs b/tests/linalg/exp.rs
index ac3f8dd7..37abc4d2 100644
--- a/tests/linalg/exp.rs
+++ b/tests/linalg/exp.rs
@@ -10,13 +10,17 @@ mod tests {
 
             let f = m.exp();
 
-            assert_eq!(f, Matrix1::new(1_f64.exp()));
+            assert!(relative_eq!(f, Matrix1::new(1_f64.exp()), epsilon = 1.0e-7));
         }
 
         {
             let m = Matrix2::new(0.0, 1.0, 0.0, 0.0);
 
-            assert_eq!(m.exp(), Matrix2::new(1.0, 1.0, 0.0, 1.0));
+            assert!(relative_eq!(
+                m.exp(),
+                Matrix2::new(1.0, 1.0, 0.0, 1.0),
+                epsilon = 1.0e-7
+            ));
         }
 
         {
@@ -35,7 +39,7 @@ mod tests {
             let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
 
             let f = Matrix2::new(m11, m12, m21, m22) / delta;
-            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+            assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
         }
 
         {
@@ -68,7 +72,7 @@ mod tests {
 
                 let f = Matrix2::new(m11, m12, m21, m22) / delta;
                 println!("a: {}", m);
-                assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+                assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
                 break;
             }
         }
@@ -91,7 +95,7 @@ mod tests {
                 e2,
             );
 
-            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+            assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
         }
     }
 }

From 0a3ee99cdb0dfa9b9202ced71ae5c277dc18b0f4 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Sun, 12 Apr 2020 11:46:00 +0200
Subject: [PATCH 10/20] Changed dimension name R to D Changed N::from_x to
 crate::convert

---
 src/linalg/exp.rs | 230 ++++++++++++++++++++++------------------------
 1 file changed, 111 insertions(+), 119 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 8822027b..f68cbe6c 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -6,25 +6,25 @@ use crate::{
         dimension::{DimMin, DimMinimum, DimName},
         DefaultAllocator,
     },
-    try_convert, ComplexField, MatrixN, RealField,
+    convert, try_convert, ComplexField, MatrixN, RealField,
 };
 
 // https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/matfuncs.py
-struct ExpmPadeHelper<N, R>
+struct ExpmPadeHelper<N, D>
 where
     N: RealField,
-    R: DimName + DimMin<R>,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    D: DimName + DimMin<D>,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     use_exact_norm: bool,
-    ident: MatrixN<N, R>,
+    ident: MatrixN<N, D>,
 
-    a: MatrixN<N, R>,
-    a2: Option<MatrixN<N, R>>,
-    a4: Option<MatrixN<N, R>>,
-    a6: Option<MatrixN<N, R>>,
-    a8: Option<MatrixN<N, R>>,
-    a10: Option<MatrixN<N, R>>,
+    a: MatrixN<N, D>,
+    a2: Option<MatrixN<N, D>>,
+    a4: Option<MatrixN<N, D>>,
+    a6: Option<MatrixN<N, D>>,
+    a8: Option<MatrixN<N, D>>,
+    a10: Option<MatrixN<N, D>>,
 
     d4_exact: Option<N>,
     d6_exact: Option<N>,
@@ -37,16 +37,16 @@ where
     d10_approx: Option<N>,
 }
 
-impl<N, R> ExpmPadeHelper<N, R>
+impl<N, D> ExpmPadeHelper<N, D>
 where
     N: RealField,
-    R: DimName + DimMin<R>,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    D: DimName + DimMin<D>,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
-    fn new(a: MatrixN<N, R>, use_exact_norm: bool) -> Self {
+    fn new(a: MatrixN<N, D>, use_exact_norm: bool) -> Self {
         ExpmPadeHelper {
             use_exact_norm,
-            ident: MatrixN::<N, R>::identity(),
+            ident: MatrixN::<N, D>::identity(),
             a,
             a2: None,
             a4: None,
@@ -64,9 +64,9 @@ where
         }
     }
 
-    fn a2(&self) -> &MatrixN<N, R> {
+    fn a2(&self) -> &MatrixN<N, D> {
         if self.a2.is_none() {
-            let ap = &self.a2 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a2 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(&self.a * &self.a);
             }
@@ -74,9 +74,9 @@ where
         self.a2.as_ref().unwrap()
     }
 
-    fn a4(&self) -> &MatrixN<N, R> {
+    fn a4(&self) -> &MatrixN<N, D> {
         if self.a4.is_none() {
-            let ap = &self.a4 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a4 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             let a2 = self.a2();
             unsafe {
                 *ap = Some(a2 * a2);
@@ -85,11 +85,11 @@ where
         self.a4.as_ref().unwrap()
     }
 
-    fn a6(&self) -> &MatrixN<N, R> {
+    fn a6(&self) -> &MatrixN<N, D> {
         if self.a6.is_none() {
             let a2 = self.a2();
             let a4 = self.a4();
-            let ap = &self.a6 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a6 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(a4 * a2);
             }
@@ -97,11 +97,11 @@ where
         self.a6.as_ref().unwrap()
     }
 
-    fn a8(&self) -> &MatrixN<N, R> {
+    fn a8(&self) -> &MatrixN<N, D> {
         if self.a8.is_none() {
             let a2 = self.a2();
             let a6 = self.a6();
-            let ap = &self.a8 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a8 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(a6 * a2);
             }
@@ -109,11 +109,11 @@ where
         self.a8.as_ref().unwrap()
     }
 
-    fn a10(&mut self) -> &MatrixN<N, R> {
+    fn a10(&mut self) -> &MatrixN<N, D> {
         if self.a10.is_none() {
             let a4 = self.a4();
             let a6 = self.a6();
-            let ap = &self.a10 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a10 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(a6 * a4);
             }
@@ -123,28 +123,28 @@ where
 
     fn d4_tight(&mut self) -> N {
         if self.d4_exact.is_none() {
-            self.d4_exact = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
+            self.d4_exact = Some(one_norm(self.a4()).powf(convert(0.25)));
         }
         self.d4_exact.unwrap()
     }
 
     fn d6_tight(&mut self) -> N {
         if self.d6_exact.is_none() {
-            self.d6_exact = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_exact = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
         }
         self.d6_exact.unwrap()
     }
 
     fn d8_tight(&mut self) -> N {
         if self.d8_exact.is_none() {
-            self.d8_exact = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_exact = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
         }
         self.d8_exact.unwrap()
     }
 
     fn d10_tight(&mut self) -> N {
         if self.d10_exact.is_none() {
-            self.d10_exact = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_exact = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
         }
         self.d10_exact.unwrap()
     }
@@ -159,7 +159,7 @@ where
         }
 
         if self.d4_approx.is_none() {
-            self.d4_approx = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
+            self.d4_approx = Some(one_norm(self.a4()).powf(convert(0.25)));
         }
 
         self.d4_approx.unwrap()
@@ -175,7 +175,7 @@ where
         }
 
         if self.d6_approx.is_none() {
-            self.d6_approx = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_approx = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
         }
 
         self.d6_approx.unwrap()
@@ -191,7 +191,7 @@ where
         }
 
         if self.d8_approx.is_none() {
-            self.d8_approx = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_approx = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
         }
 
         self.d8_approx.unwrap()
@@ -207,48 +207,43 @@ where
         }
 
         if self.d10_approx.is_none() {
-            self.d10_approx = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_approx = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
         }
 
         self.d10_approx.unwrap()
     }
 
-    fn pade3(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(120.0).unwrap(),
-            N::from_f64(60.0).unwrap(),
-            N::from_f64(12.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
-        ];
+    fn pade3(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 4] = [convert(120.0), convert(60.0), convert(12.0), convert(1.0)];
         let u = &self.a * (self.a2() * b[3] + &self.ident * b[1]);
         let v = self.a2() * b[2] + &self.ident * b[0];
         (u, v)
     }
 
-    fn pade5(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(30240.0).unwrap(),
-            N::from_f64(15120.0).unwrap(),
-            N::from_f64(3360.0).unwrap(),
-            N::from_f64(420.0).unwrap(),
-            N::from_f64(30.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade5(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 6] = [
+            convert(30240.0),
+            convert(15120.0),
+            convert(3360.0),
+            convert(420.0),
+            convert(30.0),
+            convert(1.0),
         ];
         let u = &self.a * (self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
         let v = self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
         (u, v)
     }
 
-    fn pade7(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(17297280.0).unwrap(),
-            N::from_f64(8648640.0).unwrap(),
-            N::from_f64(1995840.0).unwrap(),
-            N::from_f64(277200.0).unwrap(),
-            N::from_f64(25200.0).unwrap(),
-            N::from_f64(1512.0).unwrap(),
-            N::from_f64(56.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade7(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 8] = [
+            convert(17297280.0),
+            convert(8648640.0),
+            convert(1995840.0),
+            convert(277200.0),
+            convert(25200.0),
+            convert(1512.0),
+            convert(56.0),
+            convert(1.0),
         ];
         let u =
             &self.a * (self.a6() * b[7] + self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
@@ -256,18 +251,18 @@ where
         (u, v)
     }
 
-    fn pade9(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(17643225600.0).unwrap(),
-            N::from_f64(8821612800.0).unwrap(),
-            N::from_f64(2075673600.0).unwrap(),
-            N::from_f64(302702400.0).unwrap(),
-            N::from_f64(30270240.0).unwrap(),
-            N::from_f64(2162160.0).unwrap(),
-            N::from_f64(110880.0).unwrap(),
-            N::from_f64(3960.0).unwrap(),
-            N::from_f64(90.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade9(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 10] = [
+            convert(17643225600.0),
+            convert(8821612800.0),
+            convert(2075673600.0),
+            convert(302702400.0),
+            convert(30270240.0),
+            convert(2162160.0),
+            convert(110880.0),
+            convert(3960.0),
+            convert(90.0),
+            convert(1.0),
         ];
         let u = &self.a
             * (self.a8() * b[9]
@@ -283,29 +278,29 @@ where
         (u, v)
     }
 
-    fn pade13_scaled(&mut self, s: u64) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(64764752532480000.0).unwrap(),
-            N::from_f64(32382376266240000.0).unwrap(),
-            N::from_f64(7771770303897600.0).unwrap(),
-            N::from_f64(1187353796428800.0).unwrap(),
-            N::from_f64(129060195264000.0).unwrap(),
-            N::from_f64(10559470521600.0).unwrap(),
-            N::from_f64(670442572800.0).unwrap(),
-            N::from_f64(33522128640.0).unwrap(),
-            N::from_f64(1323241920.0).unwrap(),
-            N::from_f64(40840800.0).unwrap(),
-            N::from_f64(960960.0).unwrap(),
-            N::from_f64(16380.0).unwrap(),
-            N::from_f64(182.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade13_scaled(&mut self, s: u64) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 14] = [
+            convert(64764752532480000.0),
+            convert(32382376266240000.0),
+            convert(7771770303897600.0),
+            convert(1187353796428800.0),
+            convert(129060195264000.0),
+            convert(10559470521600.0),
+            convert(670442572800.0),
+            convert(33522128640.0),
+            convert(1323241920.0),
+            convert(40840800.0),
+            convert(960960.0),
+            convert(16380.0),
+            convert(182.0),
+            convert(1.0),
         ];
         let s = s as f64;
 
-        let mb = &self.a * N::from_f64(2.0.powf(-s)).unwrap();
-        let mb2 = self.a2() * N::from_f64(2.0.powf(-2.0 * s)).unwrap();
-        let mb4 = self.a4() * N::from_f64(2.0.powf(-4.0 * s)).unwrap();
-        let mb6 = self.a6() * N::from_f64(2.0.powf(-6.0 * s)).unwrap();
+        let mb = &self.a * convert::<f64, N>(2.0_f64.powf(-s));
+        let mb2 = self.a2() * convert::<f64, N>(2.0_f64.powf(-2.0 * s));
+        let mb4 = self.a4() * convert::<f64, N>(2.0.powf(-4.0 * s));
+        let mb6 = self.a6() * convert::<f64, N>(2.0.powf(-6.0 * s));
 
         let u2 = &mb6 * (&mb6 * b[13] + &mb4 * b[11] + &mb2 * b[9]);
         let u = &mb * (&u2 + &mb6 * b[7] + &mb4 * b[5] + &mb2 * b[3] + &self.ident * b[1]);
@@ -323,13 +318,13 @@ fn factorial(n: u128) -> u128 {
 }
 
 /// Compute the 1-norm of a non-negative integer power of a non-negative matrix.
-fn onenorm_matrix_power_nonm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
+fn onenorm_matrix_power_nonm<N, D>(a: &MatrixN<N, D>, p: u64) -> N
 where
     N: RealField,
-    R: DimName,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
+    D: DimName,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
-    let mut v = crate::VectorN::<N, R>::repeat(N::from_f64(1.0).unwrap());
+    let mut v = crate::VectorN::<N, D>::repeat(convert(1.0));
     let m = a.transpose();
 
     for _ in 0..p {
@@ -339,11 +334,11 @@ where
     v.max()
 }
 
-fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
+fn ell<N, D>(a: &MatrixN<N, D>, m: u64) -> u64
 where
     N: RealField,
-    R: DimName,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
+    D: DimName,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
     // 2m choose m = (2m)!/(m! * (2m-m)!)
 
@@ -357,10 +352,10 @@ where
         factorial(2 * m as u128) / (factorial(m as u128) * factorial(2 * m as u128 - m as u128));
     let abs_c_recip = choose_2m_m * factorial(2 * m as u128 + 1);
     let alpha = a_abs_onenorm / one_norm(a);
-    let alpha = alpha / N::from_u128(abs_c_recip).unwrap();
+    let alpha: f64 = try_convert(alpha).unwrap() / abs_c_recip as f64;
 
-    let u = N::from_f64(2_f64.powf(-53.0)).unwrap();
-    let log2_alpha_div_u = try_convert((alpha / u).log2()).unwrap();
+    let u = 2_f64.powf(-53.0);
+    let log2_alpha_div_u = (alpha / u).log2();
     let value = (log2_alpha_div_u / (2.0 * m as f64)).ceil();
     if value > 0.0 {
         value as u64
@@ -369,11 +364,11 @@ where
     }
 }
 
-fn solve_p_q<N, R>(u: MatrixN<N, R>, v: MatrixN<N, R>) -> MatrixN<N, R>
+fn solve_p_q<N, D>(u: MatrixN<N, D>, v: MatrixN<N, D>) -> MatrixN<N, D>
 where
     N: ComplexField,
-    R: DimMin<R, Output = R> + DimName,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    D: DimMin<D, Output = D> + DimName,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     let p = &u + &v;
     let q = &v - &u;
@@ -381,11 +376,11 @@ where
     q.lu().solve(&p).unwrap()
 }
 
-fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
+fn one_norm<N, D>(m: &MatrixN<N, D>) -> N
 where
     N: RealField,
-    R: DimName,
-    DefaultAllocator: Allocator<N, R, R>,
+    D: DimName,
+    DefaultAllocator: Allocator<N, D, D>,
 {
     let mut col_sums = vec![N::zero(); m.ncols()];
     for i in 0..m.ncols() {
@@ -399,11 +394,11 @@ where
     max
 }
 
-impl<N: RealField, R> MatrixN<N, R>
+impl<N: RealField, D> MatrixN<N, D>
 where
-    R: DimMin<R, Output = R> + DimName,
+    D: DimMin<D, Output = D> + DimName,
     DefaultAllocator:
-        Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>> + Allocator<N, R>,
+        Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>> + Allocator<N, D>,
 {
     /// Computes exponential of this matrix
     pub fn exp(&self) -> Self {
@@ -415,30 +410,30 @@ where
         let mut h = ExpmPadeHelper::new(self.clone(), true);
 
         let eta_1 = N::max(h.d4_loose(), h.d6_loose());
-        if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
+        if eta_1 < convert(1.495585217958292e-002) && ell(&h.a, 3) == 0 {
             let (u, v) = h.pade3();
             return solve_p_q(u, v);
         }
 
         let eta_2 = N::max(h.d4_tight(), h.d6_loose());
-        if eta_2 < N::from_f64(2.539398330063230e-001).unwrap() && ell(&h.a, 5) == 0 {
+        if eta_2 < convert(2.539398330063230e-001) && ell(&h.a, 5) == 0 {
             let (u, v) = h.pade5();
             return solve_p_q(u, v);
         }
 
         let eta_3 = N::max(h.d6_tight(), h.d8_loose());
-        if eta_3 < N::from_f64(9.504178996162932e-001).unwrap() && ell(&h.a, 7) == 0 {
+        if eta_3 < convert(9.504178996162932e-001) && ell(&h.a, 7) == 0 {
             let (u, v) = h.pade7();
             return solve_p_q(u, v);
         }
-        if eta_3 < N::from_f64(2.097847961257068e+000).unwrap() && ell(&h.a, 9) == 0 {
+        if eta_3 < convert(2.097847961257068e+000) && ell(&h.a, 9) == 0 {
             let (u, v) = h.pade9();
             return solve_p_q(u, v);
         }
 
         let eta_4 = N::max(h.d8_loose(), h.d10_loose());
         let eta_5 = N::min(eta_3, eta_4);
-        let theta_13 = N::from_f64(4.25).unwrap();
+        let theta_13 = convert(4.25);
 
         let mut s = if eta_5 == N::zero() {
             0
@@ -452,10 +447,7 @@ where
             }
         };
 
-        s += ell(
-            &(&h.a * N::from_f64(2.0_f64.powf(-(s as f64))).unwrap()),
-            13,
-        );
+        s += ell(&(&h.a * convert::<f64, N>(2.0_f64.powf(-(s as f64)))), 13);
 
         let (u, v) = h.pade13_scaled(s);
         let mut x = solve_p_q(u, v);

From e914afe2afd243b568b2d8c41e09a621f465851b Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Sun, 12 Apr 2020 11:59:06 +0200
Subject: [PATCH 11/20] Added support for dynamic matrices

---
 src/linalg/exp.rs   | 23 +++++++++++++----------
 tests/linalg/exp.rs | 30 +++++++++++++++++++++++++++++-
 2 files changed, 42 insertions(+), 11 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index f68cbe6c..e07ddd2f 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -3,7 +3,8 @@
 use crate::{
     base::{
         allocator::Allocator,
-        dimension::{DimMin, DimMinimum, DimName},
+        dimension::{Dim, DimMin, DimMinimum, U1},
+        storage::Storage,
         DefaultAllocator,
     },
     convert, try_convert, ComplexField, MatrixN, RealField,
@@ -13,7 +14,7 @@ use crate::{
 struct ExpmPadeHelper<N, D>
 where
     N: RealField,
-    D: DimName + DimMin<D>,
+    D: DimMin<D>,
     DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     use_exact_norm: bool,
@@ -40,13 +41,14 @@ where
 impl<N, D> ExpmPadeHelper<N, D>
 where
     N: RealField,
-    D: DimName + DimMin<D>,
+    D: DimMin<D>,
     DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     fn new(a: MatrixN<N, D>, use_exact_norm: bool) -> Self {
+        let (nrows, ncols) = a.data.shape();
         ExpmPadeHelper {
             use_exact_norm,
-            ident: MatrixN::<N, D>::identity(),
+            ident: MatrixN::<N, D>::identity_generic(nrows, ncols),
             a,
             a2: None,
             a4: None,
@@ -321,10 +323,11 @@ fn factorial(n: u128) -> u128 {
 fn onenorm_matrix_power_nonm<N, D>(a: &MatrixN<N, D>, p: u64) -> N
 where
     N: RealField,
-    D: DimName,
+    D: Dim,
     DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
-    let mut v = crate::VectorN::<N, D>::repeat(convert(1.0));
+    let nrows = a.data.shape().0;
+    let mut v = crate::VectorN::<N, D>::repeat_generic(nrows, U1, convert(1.0));
     let m = a.transpose();
 
     for _ in 0..p {
@@ -337,7 +340,7 @@ where
 fn ell<N, D>(a: &MatrixN<N, D>, m: u64) -> u64
 where
     N: RealField,
-    D: DimName,
+    D: Dim,
     DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
     // 2m choose m = (2m)!/(m! * (2m-m)!)
@@ -367,7 +370,7 @@ where
 fn solve_p_q<N, D>(u: MatrixN<N, D>, v: MatrixN<N, D>) -> MatrixN<N, D>
 where
     N: ComplexField,
-    D: DimMin<D, Output = D> + DimName,
+    D: DimMin<D, Output = D>,
     DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     let p = &u + &v;
@@ -379,7 +382,7 @@ where
 fn one_norm<N, D>(m: &MatrixN<N, D>) -> N
 where
     N: RealField,
-    D: DimName,
+    D: Dim,
     DefaultAllocator: Allocator<N, D, D>,
 {
     let mut col_sums = vec![N::zero(); m.ncols()];
@@ -396,7 +399,7 @@ where
 
 impl<N: RealField, D> MatrixN<N, D>
 where
-    D: DimMin<D, Output = D> + DimName,
+    D: DimMin<D, Output = D>,
     DefaultAllocator:
         Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>> + Allocator<N, D>,
 {
diff --git a/tests/linalg/exp.rs b/tests/linalg/exp.rs
index 37abc4d2..75122107 100644
--- a/tests/linalg/exp.rs
+++ b/tests/linalg/exp.rs
@@ -2,7 +2,7 @@
 mod tests {
     //https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/tests/test_matfuncs.py
     #[test]
-    fn exp() {
+    fn exp_static() {
         use nalgebra::{Matrix1, Matrix2, Matrix3};
 
         {
@@ -98,4 +98,32 @@ mod tests {
             assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
         }
     }
+
+    #[test]
+    fn exp_dynamic() {
+        use nalgebra::DMatrix;
+
+        let m = DMatrix::from_row_slice(3, 3, &[1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0]);
+
+        let e1 = 1.0_f64.exp();
+        let e2 = 2.0_f64.exp();
+
+        let f = DMatrix::from_row_slice(
+            3,
+            3,
+            &[
+                e1,
+                3.0 * e1,
+                15.0 * (e2 - 2.0 * e1),
+                0.0,
+                e1,
+                5.0 * (e2 - e1),
+                0.0,
+                0.0,
+                e2,
+            ],
+        );
+
+        assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
+    }
 }

From 583a8fb1102f075e4457208f3c1cf68b21707384 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Sun, 12 Apr 2020 13:27:11 +0200
Subject: [PATCH 12/20] Removed unsafe immutable hack

---
 src/linalg/exp.rs | 146 ++++++++++++++++++++++++++--------------------
 1 file changed, 83 insertions(+), 63 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index e07ddd2f..86bd81a7 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -66,87 +66,78 @@ where
         }
     }
 
-    fn a2(&self) -> &MatrixN<N, D> {
+    fn calc_a2(&mut self) {
         if self.a2.is_none() {
-            let ap = &self.a2 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(&self.a * &self.a);
-            }
+            self.a2 = Some(&self.a * &self.a);
         }
-        self.a2.as_ref().unwrap()
     }
 
-    fn a4(&self) -> &MatrixN<N, D> {
+    fn calc_a4(&mut self) {
         if self.a4.is_none() {
-            let ap = &self.a4 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            let a2 = self.a2();
-            unsafe {
-                *ap = Some(a2 * a2);
-            }
+            self.calc_a2();
+            let a2 = self.a2.as_ref().unwrap();
+            self.a4 = Some(a2 * a2);
         }
-        self.a4.as_ref().unwrap()
     }
 
-    fn a6(&self) -> &MatrixN<N, D> {
+    fn calc_a6(&mut self) {
         if self.a6.is_none() {
-            let a2 = self.a2();
-            let a4 = self.a4();
-            let ap = &self.a6 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(a4 * a2);
-            }
+            self.calc_a2();
+            self.calc_a4();
+            let a2 = self.a2.as_ref().unwrap();
+            let a4 = self.a4.as_ref().unwrap();
+            self.a6 = Some(a4 * a2);
         }
-        self.a6.as_ref().unwrap()
     }
 
-    fn a8(&self) -> &MatrixN<N, D> {
+    fn calc_a8(&mut self) {
         if self.a8.is_none() {
-            let a2 = self.a2();
-            let a6 = self.a6();
-            let ap = &self.a8 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(a6 * a2);
-            }
+            self.calc_a2();
+            self.calc_a6();
+            let a2 = self.a2.as_ref().unwrap();
+            let a6 = self.a6.as_ref().unwrap();
+            self.a8 = Some(a6 * a2);
         }
-        self.a8.as_ref().unwrap()
     }
 
-    fn a10(&mut self) -> &MatrixN<N, D> {
+    fn calc_a10(&mut self) {
         if self.a10.is_none() {
-            let a4 = self.a4();
-            let a6 = self.a6();
-            let ap = &self.a10 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(a6 * a4);
-            }
+            self.calc_a4();
+            self.calc_a6();
+            let a4 = self.a4.as_ref().unwrap();
+            let a6 = self.a6.as_ref().unwrap();
+            self.a10 = Some(a6 * a4);
         }
-        self.a10.as_ref().unwrap()
     }
 
     fn d4_tight(&mut self) -> N {
         if self.d4_exact.is_none() {
-            self.d4_exact = Some(one_norm(self.a4()).powf(convert(0.25)));
+            self.calc_a4();
+            self.d4_exact = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25)));
         }
         self.d4_exact.unwrap()
     }
 
     fn d6_tight(&mut self) -> N {
         if self.d6_exact.is_none() {
-            self.d6_exact = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
+            self.calc_a6();
+            self.d6_exact = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0)));
         }
         self.d6_exact.unwrap()
     }
 
     fn d8_tight(&mut self) -> N {
         if self.d8_exact.is_none() {
-            self.d8_exact = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
+            self.calc_a8();
+            self.d8_exact = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0)));
         }
         self.d8_exact.unwrap()
     }
 
     fn d10_tight(&mut self) -> N {
         if self.d10_exact.is_none() {
-            self.d10_exact = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
+            self.calc_a10();
+            self.d10_exact = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0)));
         }
         self.d10_exact.unwrap()
     }
@@ -161,7 +152,8 @@ where
         }
 
         if self.d4_approx.is_none() {
-            self.d4_approx = Some(one_norm(self.a4()).powf(convert(0.25)));
+            self.calc_a4();
+            self.d4_approx = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25)));
         }
 
         self.d4_approx.unwrap()
@@ -177,7 +169,8 @@ where
         }
 
         if self.d6_approx.is_none() {
-            self.d6_approx = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
+            self.calc_a6();
+            self.d6_approx = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0)));
         }
 
         self.d6_approx.unwrap()
@@ -193,7 +186,8 @@ where
         }
 
         if self.d8_approx.is_none() {
-            self.d8_approx = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
+            self.calc_a8();
+            self.d8_approx = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0)));
         }
 
         self.d8_approx.unwrap()
@@ -209,7 +203,8 @@ where
         }
 
         if self.d10_approx.is_none() {
-            self.d10_approx = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
+            self.calc_a10();
+            self.d10_approx = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0)));
         }
 
         self.d10_approx.unwrap()
@@ -217,8 +212,10 @@ where
 
     fn pade3(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
         let b: [N; 4] = [convert(120.0), convert(60.0), convert(12.0), convert(1.0)];
-        let u = &self.a * (self.a2() * b[3] + &self.ident * b[1]);
-        let v = self.a2() * b[2] + &self.ident * b[0];
+        self.calc_a2();
+        let a2 = self.a2.as_ref().unwrap();
+        let u = &self.a * (a2 * b[3] + &self.ident * b[1]);
+        let v = a2 * b[2] + &self.ident * b[0];
         (u, v)
     }
 
@@ -231,8 +228,15 @@ where
             convert(30.0),
             convert(1.0),
         ];
-        let u = &self.a * (self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
-        let v = self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        self.calc_a2();
+        self.calc_a6();
+        let u = &self.a
+            * (self.a4.as_ref().unwrap() * b[5]
+                + self.a2.as_ref().unwrap() * b[3]
+                + &self.ident * b[1]);
+        let v = self.a4.as_ref().unwrap() * b[4]
+            + self.a2.as_ref().unwrap() * b[2]
+            + &self.ident * b[0];
         (u, v)
     }
 
@@ -247,9 +251,18 @@ where
             convert(56.0),
             convert(1.0),
         ];
-        let u =
-            &self.a * (self.a6() * b[7] + self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
-        let v = self.a6() * b[6] + self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        self.calc_a2();
+        self.calc_a4();
+        self.calc_a6();
+        let u = &self.a
+            * (self.a6.as_ref().unwrap() * b[7]
+                + self.a4.as_ref().unwrap() * b[5]
+                + self.a2.as_ref().unwrap() * b[3]
+                + &self.ident * b[1]);
+        let v = self.a6.as_ref().unwrap() * b[6]
+            + self.a4.as_ref().unwrap() * b[4]
+            + self.a2.as_ref().unwrap() * b[2]
+            + &self.ident * b[0];
         (u, v)
     }
 
@@ -266,16 +279,20 @@ where
             convert(90.0),
             convert(1.0),
         ];
+        self.calc_a2();
+        self.calc_a4();
+        self.calc_a6();
+        self.calc_a8();
         let u = &self.a
-            * (self.a8() * b[9]
-                + self.a6() * b[7]
-                + self.a4() * b[5]
-                + self.a2() * b[3]
+            * (self.a8.as_ref().unwrap() * b[9]
+                + self.a6.as_ref().unwrap() * b[7]
+                + self.a4.as_ref().unwrap() * b[5]
+                + self.a2.as_ref().unwrap() * b[3]
                 + &self.ident * b[1]);
-        let v = self.a8() * b[8]
-            + self.a6() * b[6]
-            + self.a4() * b[4]
-            + self.a2() * b[2]
+        let v = self.a8.as_ref().unwrap() * b[8]
+            + self.a6.as_ref().unwrap() * b[6]
+            + self.a4.as_ref().unwrap() * b[4]
+            + self.a2.as_ref().unwrap() * b[2]
             + &self.ident * b[0];
         (u, v)
     }
@@ -300,9 +317,12 @@ where
         let s = s as f64;
 
         let mb = &self.a * convert::<f64, N>(2.0_f64.powf(-s));
-        let mb2 = self.a2() * convert::<f64, N>(2.0_f64.powf(-2.0 * s));
-        let mb4 = self.a4() * convert::<f64, N>(2.0.powf(-4.0 * s));
-        let mb6 = self.a6() * convert::<f64, N>(2.0.powf(-6.0 * s));
+        self.calc_a2();
+        self.calc_a4();
+        self.calc_a6();
+        let mb2 = self.a2.as_ref().unwrap() * convert::<f64, N>(2.0_f64.powf(-2.0 * s));
+        let mb4 = self.a4.as_ref().unwrap() * convert::<f64, N>(2.0.powf(-4.0 * s));
+        let mb6 = self.a6.as_ref().unwrap() * convert::<f64, N>(2.0.powf(-6.0 * s));
 
         let u2 = &mb6 * (&mb6 * b[13] + &mb4 * b[11] + &mb2 * b[9]);
         let u = &mb * (&u2 + &mb6 * b[7] + &mb4 * b[5] + &mb2 * b[3] + &self.ident * b[1]);

From a50e567fcd07746285a62bc17fc84e45162bfaf6 Mon Sep 17 00:00:00 2001
From: Max Orok <maxwellorok@gmail.com>
Date: Mon, 20 Apr 2020 10:05:39 -0400
Subject: [PATCH 13/20] Fix alga cfg gate for some quickcheck tests.

---
 tests/core/conversion.rs | 2 +-
 tests/core/matrix.rs     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/core/conversion.rs b/tests/core/conversion.rs
index 74aa41e9..b7a8c5f8 100644
--- a/tests/core/conversion.rs
+++ b/tests/core/conversion.rs
@@ -1,4 +1,4 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(all(feature = "arbitrary", feature = "alga"))]
 use alga::linear::Transformation;
 use na::{
     self, Affine3, Isometry3, Matrix2, Matrix2x3, Matrix2x4, Matrix2x5, Matrix2x6, Matrix3,
diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs
index 9e25db58..6ade807f 100644
--- a/tests/core/matrix.rs
+++ b/tests/core/matrix.rs
@@ -945,7 +945,7 @@ mod normalization_tests {
     }
 }
 
-#[cfg(feature = "arbitrary")]
+#[cfg(all(feature = "arbitrary", feature = "alga"))]
 // FIXME: move this to alga ?
 mod finite_dim_inner_space_tests {
     use super::*;

From f4c089776438094e33be29405ffe0a28d3c69f0a Mon Sep 17 00:00:00 2001
From: sebcrozet <developer@crozet.re>
Date: Tue, 21 Apr 2020 10:19:03 +0200
Subject: [PATCH 14/20] Fix compilation of matrix exponential when targetting
 no-std.

---
 src/linalg/exp.rs | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 86bd81a7..df375ac9 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -42,7 +42,8 @@ impl<N, D> ExpmPadeHelper<N, D>
 where
     N: RealField,
     D: DimMin<D>,
-    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
+    DefaultAllocator:
+        Allocator<N, D, D> + Allocator<N, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     fn new(a: MatrixN<N, D>, use_exact_norm: bool) -> Self {
         let (nrows, ncols) = a.data.shape();
@@ -405,15 +406,13 @@ where
     D: Dim,
     DefaultAllocator: Allocator<N, D, D>,
 {
-    let mut col_sums = vec![N::zero(); m.ncols()];
+    let mut max = N::zero();
+
     for i in 0..m.ncols() {
         let col = m.column(i);
-        col.iter().for_each(|v| col_sums[i] += v.abs());
-    }
-    let mut max = col_sums[0];
-    for i in 1..col_sums.len() {
-        max = N::max(max, col_sums[i]);
+        max = max.max(col.iter().fold(N::zero(), |a, b| a + b.abs()));
     }
+
     max
 }
 
@@ -427,7 +426,7 @@ where
     pub fn exp(&self) -> Self {
         // Simple case
         if self.nrows() == 1 {
-            return self.clone().map(|v| v.exp());
+            return self.map(|v| v.exp());
         }
 
         let mut h = ExpmPadeHelper::new(self.clone(), true);

From 6dc739b70b44a0806d3b10bcc5646cea94696a95 Mon Sep 17 00:00:00 2001
From: sebcrozet <developer@crozet.re>
Date: Tue, 21 Apr 2020 10:37:14 +0200
Subject: [PATCH 15/20] Remove unused allocator bound.

---
 src/linalg/exp.rs | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index df375ac9..c11fd757 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -42,8 +42,7 @@ impl<N, D> ExpmPadeHelper<N, D>
 where
     N: RealField,
     D: DimMin<D>,
-    DefaultAllocator:
-        Allocator<N, D, D> + Allocator<N, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     fn new(a: MatrixN<N, D>, use_exact_norm: bool) -> Self {
         let (nrows, ncols) = a.data.shape();

From 7363caba172f1964d7741691470caeb3a43951fa Mon Sep 17 00:00:00 2001
From: Chris Hodapp <hodapp87@gmail.com>
Date: Tue, 21 Apr 2020 10:23:00 -0400
Subject: [PATCH 16/20] Fix doc typo (v is right-singular vectors, not left)

---
 src/linalg/svd.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs
index 6d74a438..4be0ecdb 100644
--- a/src/linalg/svd.rs
+++ b/src/linalg/svd.rs
@@ -95,7 +95,7 @@ where
     /// # Arguments
     ///
     /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
-    /// * `compute_v` − set this to `true` to enable the computation of left-singular vectors.
+    /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
     /// * `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
@@ -626,7 +626,7 @@ where
     /// # Arguments
     ///
     /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
-    /// * `compute_v` − set this to `true` to enable the computation of left-singular vectors.
+    /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
     /// * `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

From 3359e25435b7ecbd3632ab38769840413c4fbd7b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= <developer@crozet.re>
Date: Sat, 23 May 2020 15:13:13 +0200
Subject: [PATCH 17/20] Cholesky: add unchecked construction compatible with
 AoSoA SIMD.

---
 src/linalg/cholesky.rs |  94 ++++++++----
 src/linalg/solve.rs    | 334 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 398 insertions(+), 30 deletions(-)

diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs
index da3549d7..cf338a95 100644
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@@ -3,6 +3,7 @@ use serde::{Deserialize, Serialize};
 
 use num::One;
 use simba::scalar::ComplexField;
+use simba::simd::SimdComplexField;
 
 use crate::allocator::Allocator;
 use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix, Vector};
@@ -23,29 +24,26 @@ use crate::storage::{Storage, StorageMut};
          MatrixN<N, D>: Deserialize<'de>"))
 )]
 #[derive(Clone, Debug)]
-pub struct Cholesky<N: ComplexField, D: Dim>
-where
-    DefaultAllocator: Allocator<N, D, D>,
+pub struct Cholesky<N: SimdComplexField, D: Dim>
+where DefaultAllocator: Allocator<N, D, D>
 {
     chol: MatrixN<N, D>,
 }
 
-impl<N: ComplexField, D: Dim> Copy for Cholesky<N, D>
+impl<N: SimdComplexField, D: Dim> Copy for Cholesky<N, D>
 where
     DefaultAllocator: Allocator<N, D, D>,
     MatrixN<N, D>: Copy,
 {
 }
 
-impl<N: ComplexField, D: DimSub<Dynamic>> Cholesky<N, D>
-where
-    DefaultAllocator: Allocator<N, D, D>,
+impl<N: SimdComplexField, D: Dim> Cholesky<N, D>
+where DefaultAllocator: Allocator<N, D, D>
 {
-    /// Attempts to compute the Cholesky decomposition of `matrix`.
+    /// Computes the Cholesky decomposition of `matrix` without checking that the matrix is definite-positive.
     ///
-    /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed
-    /// to be symmetric and only the lower-triangular part is read.
-    pub fn new(mut matrix: MatrixN<N, D>) -> Option<Self> {
+    /// If the input matrix is not definite-positive, the decomposition may contain trash values (Inf, NaN, etc.)
+    pub fn new_unchecked(mut matrix: MatrixN<N, D>) -> Self {
         assert!(matrix.is_square(), "The input matrix must be square.");
 
         let n = matrix.nrows();
@@ -57,29 +55,21 @@ where
                 let (mut col_j, col_k) = matrix.columns_range_pair_mut(j, k);
                 let mut col_j = col_j.rows_range_mut(j..);
                 let col_k = col_k.rows_range(j..);
-
-                col_j.axpy(factor.conjugate(), &col_k, N::one());
+                col_j.axpy(factor.simd_conjugate(), &col_k, N::one());
             }
 
             let diag = unsafe { *matrix.get_unchecked((j, j)) };
-            if !diag.is_zero() {
-                if let Some(denom) = diag.try_sqrt() {
-                    unsafe {
-                        *matrix.get_unchecked_mut((j, j)) = denom;
-                    }
+            let denom = diag.simd_sqrt();
 
-                    let mut col = matrix.slice_range_mut(j + 1.., j);
-                    col /= denom;
-                    continue;
-                }
+            unsafe {
+                *matrix.get_unchecked_mut((j, j)) = denom;
             }
 
-            // The diagonal element is either zero or its square root could not
-            // be taken (e.g. for negative real numbers).
-            return None;
+            let mut col = matrix.slice_range_mut(j + 1.., j);
+            col /= denom;
         }
 
-        Some(Cholesky { chol: matrix })
+        Cholesky { chol: matrix }
     }
 
     /// Retrieves the lower-triangular factor of the Cholesky decomposition with its strictly
@@ -121,8 +111,8 @@ where
         S2: StorageMut<N, R2, C2>,
         ShapeConstraint: SameNumberOfRows<R2, D>,
     {
-        let _ = self.chol.solve_lower_triangular_mut(b);
-        let _ = self.chol.ad_solve_lower_triangular_mut(b);
+        self.chol.solve_lower_triangular_unchecked_mut(b);
+        self.chol.ad_solve_lower_triangular_unchecked_mut(b);
     }
 
     /// Returns the solution of the system `self * x = b` where `self` is the decomposed matrix and
@@ -146,6 +136,51 @@ where
         self.solve_mut(&mut res);
         res
     }
+}
+
+impl<N: ComplexField, D: Dim> Cholesky<N, D>
+where DefaultAllocator: Allocator<N, D, D>
+{
+    /// Attempts to compute the Cholesky decomposition of `matrix`.
+    ///
+    /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed
+    /// to be symmetric and only the lower-triangular part is read.
+    pub fn new(mut matrix: MatrixN<N, D>) -> Option<Self> {
+        assert!(matrix.is_square(), "The input matrix must be square.");
+
+        let n = matrix.nrows();
+
+        for j in 0..n {
+            for k in 0..j {
+                let factor = unsafe { -*matrix.get_unchecked((j, k)) };
+
+                let (mut col_j, col_k) = matrix.columns_range_pair_mut(j, k);
+                let mut col_j = col_j.rows_range_mut(j..);
+                let col_k = col_k.rows_range(j..);
+
+                col_j.axpy(factor.conjugate(), &col_k, N::one());
+            }
+
+            let diag = unsafe { *matrix.get_unchecked((j, j)) };
+            if !diag.is_zero() {
+                if let Some(denom) = diag.try_sqrt() {
+                    unsafe {
+                        *matrix.get_unchecked_mut((j, j)) = denom;
+                    }
+
+                    let mut col = matrix.slice_range_mut(j + 1.., j);
+                    col /= denom;
+                    continue;
+                }
+            }
+
+            // The diagonal element is either zero or its square root could not
+            // be taken (e.g. for negative real numbers).
+            return None;
+        }
+
+        Some(Cholesky { chol: matrix })
+    }
 
     /// Given the Cholesky decomposition of a matrix `M`, a scalar `sigma` and a vector `v`,
     /// performs a rank one update such that we end up with the decomposition of `M + sigma * (v * v.adjoint())`.
@@ -327,8 +362,7 @@ where
 }
 
 impl<N: ComplexField, D: DimSub<Dynamic>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
-where
-    DefaultAllocator: Allocator<N, D, D>,
+where DefaultAllocator: Allocator<N, D, D>
 {
     /// Attempts to compute the Cholesky decomposition of this matrix.
     ///
diff --git a/src/linalg/solve.rs b/src/linalg/solve.rs
index 56db4ade..ac5bff46 100644
--- a/src/linalg/solve.rs
+++ b/src/linalg/solve.rs
@@ -1,4 +1,5 @@
 use simba::scalar::ComplexField;
+use simba::simd::SimdComplexField;
 
 use crate::base::allocator::Allocator;
 use crate::base::constraint::{SameNumberOfRows, ShapeConstraint};
@@ -432,3 +433,336 @@ impl<N: ComplexField, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
         true
     }
 }
+
+/*
+ *
+ * SIMD-compatible unchecked versions.
+ *
+ */
+
+impl<N: SimdComplexField, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
+    /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only
+    /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[inline]
+    pub fn solve_lower_triangular_unchecked<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &Matrix<N, R2, C2, S2>,
+    ) -> MatrixMN<N, R2, C2>
+    where
+        S2: Storage<N, R2, C2>,
+        DefaultAllocator: Allocator<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let mut res = b.clone_owned();
+        self.solve_lower_triangular_unchecked_mut(&mut res);
+        res
+    }
+
+    /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only
+    /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[inline]
+    pub fn solve_upper_triangular_unchecked<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &Matrix<N, R2, C2, S2>,
+    ) -> MatrixMN<N, R2, C2>
+    where
+        S2: Storage<N, R2, C2>,
+        DefaultAllocator: Allocator<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let mut res = b.clone_owned();
+        self.solve_upper_triangular_unchecked_mut(&mut res);
+        res
+    }
+
+    /// Solves the linear system `self . x = b` where `x` is the unknown and only the
+    /// lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    pub fn solve_lower_triangular_unchecked_mut<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        S2: StorageMut<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        for i in 0..b.ncols() {
+            self.solve_lower_triangular_vector_unchecked_mut(&mut b.column_mut(i));
+        }
+    }
+
+    fn solve_lower_triangular_vector_unchecked_mut<R2: Dim, S2>(&self, b: &mut Vector<N, R2, S2>)
+    where
+        S2: StorageMut<N, R2, U1>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let dim = self.nrows();
+
+        for i in 0..dim {
+            let coeff;
+
+            unsafe {
+                let diag = *self.get_unchecked((i, i));
+                coeff = *b.vget_unchecked(i) / diag;
+                *b.vget_unchecked_mut(i) = coeff;
+            }
+
+            b.rows_range_mut(i + 1..)
+                .axpy(-coeff, &self.slice_range(i + 1.., i), N::one());
+        }
+    }
+
+    // FIXME: add the same but for solving upper-triangular.
+    /// Solves the linear system `self . x = b` where `x` is the unknown and only the
+    /// lower-triangular part of `self` is considered not-zero. The diagonal is never read as it is
+    /// assumed to be equal to `diag`. Returns `false` and does not modify its inputs if `diag` is zero.
+    pub fn solve_lower_triangular_with_diag_unchecked_mut<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &mut Matrix<N, R2, C2, S2>,
+        diag: N,
+    ) where
+        S2: StorageMut<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let dim = self.nrows();
+        let cols = b.ncols();
+
+        for k in 0..cols {
+            let mut bcol = b.column_mut(k);
+
+            for i in 0..dim - 1 {
+                let coeff = unsafe { *bcol.vget_unchecked(i) } / diag;
+                bcol.rows_range_mut(i + 1..)
+                    .axpy(-coeff, &self.slice_range(i + 1.., i), N::one());
+            }
+        }
+    }
+
+    /// Solves the linear system `self . x = b` where `x` is the unknown and only the
+    /// upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    pub fn solve_upper_triangular_unchecked_mut<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        S2: StorageMut<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        for i in 0..b.ncols() {
+            self.solve_upper_triangular_vector_unchecked_mut(&mut b.column_mut(i))
+        }
+    }
+
+    fn solve_upper_triangular_vector_unchecked_mut<R2: Dim, S2>(&self, b: &mut Vector<N, R2, S2>)
+    where
+        S2: StorageMut<N, R2, U1>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let dim = self.nrows();
+
+        for i in (0..dim).rev() {
+            let coeff;
+
+            unsafe {
+                let diag = *self.get_unchecked((i, i));
+                coeff = *b.vget_unchecked(i) / diag;
+                *b.vget_unchecked_mut(i) = coeff;
+            }
+
+            b.rows_range_mut(..i)
+                .axpy(-coeff, &self.slice_range(..i, i), N::one());
+        }
+    }
+
+    /*
+     *
+     * Transpose and adjoint versions
+     *
+     */
+    /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only
+    /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[inline]
+    pub fn tr_solve_lower_triangular_unchecked<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &Matrix<N, R2, C2, S2>,
+    ) -> MatrixMN<N, R2, C2>
+    where
+        S2: Storage<N, R2, C2>,
+        DefaultAllocator: Allocator<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let mut res = b.clone_owned();
+        self.tr_solve_lower_triangular_unchecked_mut(&mut res);
+        res
+    }
+
+    /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only
+    /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[inline]
+    pub fn tr_solve_upper_triangular_unchecked<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &Matrix<N, R2, C2, S2>,
+    ) -> MatrixMN<N, R2, C2>
+    where
+        S2: Storage<N, R2, C2>,
+        DefaultAllocator: Allocator<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let mut res = b.clone_owned();
+        self.tr_solve_upper_triangular_unchecked_mut(&mut res);
+        res
+    }
+
+    /// Solves the linear system `self.transpose() . x = b` where `x` is the unknown and only the
+    /// lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    pub fn tr_solve_lower_triangular_unchecked_mut<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        S2: StorageMut<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        for i in 0..b.ncols() {
+            self.xx_solve_lower_triangular_vector_unchecked_mut(
+                &mut b.column_mut(i),
+                |e| e,
+                |a, b| a.dot(b),
+            )
+        }
+    }
+
+    /// Solves the linear system `self.transpose() . x = b` where `x` is the unknown and only the
+    /// upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    pub fn tr_solve_upper_triangular_unchecked_mut<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        S2: StorageMut<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        for i in 0..b.ncols() {
+            self.xx_solve_upper_triangular_vector_unchecked_mut(
+                &mut b.column_mut(i),
+                |e| e,
+                |a, b| a.dot(b),
+            )
+        }
+    }
+
+    /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only
+    /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[inline]
+    pub fn ad_solve_lower_triangular_unchecked<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &Matrix<N, R2, C2, S2>,
+    ) -> MatrixMN<N, R2, C2>
+    where
+        S2: Storage<N, R2, C2>,
+        DefaultAllocator: Allocator<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let mut res = b.clone_owned();
+        self.ad_solve_lower_triangular_unchecked_mut(&mut res);
+        res
+    }
+
+    /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only
+    /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[inline]
+    pub fn ad_solve_upper_triangular_unchecked<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &Matrix<N, R2, C2, S2>,
+    ) -> MatrixMN<N, R2, C2>
+    where
+        S2: Storage<N, R2, C2>,
+        DefaultAllocator: Allocator<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let mut res = b.clone_owned();
+        self.ad_solve_upper_triangular_unchecked_mut(&mut res);
+        res
+    }
+
+    /// Solves the linear system `self.adjoint() . x = b` where `x` is the unknown and only the
+    /// lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    pub fn ad_solve_lower_triangular_unchecked_mut<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        S2: StorageMut<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        for i in 0..b.ncols() {
+            self.xx_solve_lower_triangular_vector_unchecked_mut(
+                &mut b.column_mut(i),
+                |e| e.simd_conjugate(),
+                |a, b| a.dotc(b),
+            )
+        }
+    }
+
+    /// Solves the linear system `self.adjoint() . x = b` where `x` is the unknown and only the
+    /// upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    pub fn ad_solve_upper_triangular_unchecked_mut<R2: Dim, C2: Dim, S2>(
+        &self,
+        b: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        S2: StorageMut<N, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        for i in 0..b.ncols() {
+            self.xx_solve_upper_triangular_vector_unchecked_mut(
+                &mut b.column_mut(i),
+                |e| e.simd_conjugate(),
+                |a, b| a.dotc(b),
+            )
+        }
+    }
+
+    #[inline(always)]
+    fn xx_solve_lower_triangular_vector_unchecked_mut<R2: Dim, S2>(
+        &self,
+        b: &mut Vector<N, R2, S2>,
+        conjugate: impl Fn(N) -> N,
+        dot: impl Fn(
+            &DVectorSlice<N, S::RStride, S::CStride>,
+            &DVectorSlice<N, S2::RStride, S2::CStride>,
+        ) -> N,
+    ) where
+        S2: StorageMut<N, R2, U1>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        let dim = self.nrows();
+
+        for i in (0..dim).rev() {
+            let dot = dot(&self.slice_range(i + 1.., i), &b.slice_range(i + 1.., 0));
+
+            unsafe {
+                let b_i = b.vget_unchecked_mut(i);
+                let diag = conjugate(*self.get_unchecked((i, i)));
+                *b_i = (*b_i - dot) / diag;
+            }
+        }
+    }
+
+    #[inline(always)]
+    fn xx_solve_upper_triangular_vector_unchecked_mut<R2: Dim, S2>(
+        &self,
+        b: &mut Vector<N, R2, S2>,
+        conjugate: impl Fn(N) -> N,
+        dot: impl Fn(
+            &DVectorSlice<N, S::RStride, S::CStride>,
+            &DVectorSlice<N, S2::RStride, S2::CStride>,
+        ) -> N,
+    ) where
+        S2: StorageMut<N, R2, U1>,
+        ShapeConstraint: SameNumberOfRows<R2, D>,
+    {
+        for i in 0..self.nrows() {
+            let dot = dot(&self.slice_range(..i, i), &b.slice_range(..i, 0));
+
+            unsafe {
+                let b_i = b.vget_unchecked_mut(i);
+                let diag = conjugate(*self.get_unchecked((i, i)));
+                *b_i = (*b_i - dot) / diag;
+            }
+        }
+    }
+}

From 2c2d1e4f07500d1de7ef5b456ec1346642ced971 Mon Sep 17 00:00:00 2001
From: sebcrozet <developer@crozet.re>
Date: Sun, 7 Jun 2020 09:07:25 +0200
Subject: [PATCH 18/20] Run cargo fmt.

---
 src/base/array_storage.rs      |  4 +++-
 src/base/construction_slice.rs | 12 ++++++------
 src/geometry/quaternion.rs     |  2 +-
 src/linalg/cholesky.rs         | 12 ++++++++----
 4 files changed, 18 insertions(+), 12 deletions(-)

diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs
index 8aa1a7bf..1743f06b 100644
--- a/src/base/array_storage.rs
+++ b/src/base/array_storage.rs
@@ -57,7 +57,9 @@ where
     N: Default,
 {
     fn default() -> Self {
-        ArrayStorage { data: Default::default() }
+        ArrayStorage {
+            data: Default::default(),
+        }
     }
 }
 
diff --git a/src/base/construction_slice.rs b/src/base/construction_slice.rs
index d63b374b..06b8be9f 100644
--- a/src/base/construction_slice.rs
+++ b/src/base/construction_slice.rs
@@ -220,9 +220,9 @@ macro_rules! impl_constructors(
 
 // FIXME: this is not very pretty. We could find a better call syntax.
 impl_constructors!(R, C;                         // Arguments for Matrix<N, ..., S>
-                   => R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
-                   R::name(), C::name();         // Arguments for `_generic` constructors.
-                   ); // Arguments for non-generic constructors.
+=> R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
+R::name(), C::name();         // Arguments for `_generic` constructors.
+); // Arguments for non-generic constructors.
 
 impl_constructors!(R, Dynamic;
                    => R: DimName;
@@ -279,9 +279,9 @@ macro_rules! impl_constructors_mut(
 
 // FIXME: this is not very pretty. We could find a better call syntax.
 impl_constructors_mut!(R, C;                         // Arguments for Matrix<N, ..., S>
-                       => R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
-                       R::name(), C::name();         // Arguments for `_generic` constructors.
-                       ); // Arguments for non-generic constructors.
+=> R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
+R::name(), C::name();         // Arguments for `_generic` constructors.
+); // Arguments for non-generic constructors.
 
 impl_constructors_mut!(R, Dynamic;
                        => R: DimName;
diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs
index c4f1617a..0c996ede 100755
--- a/src/geometry/quaternion.rs
+++ b/src/geometry/quaternion.rs
@@ -36,7 +36,7 @@ pub struct Quaternion<N: Scalar + SimdValue> {
 impl<N: RealField> Default for Quaternion<N> {
     fn default() -> Self {
         Quaternion {
-            coords: Vector4::zeros()
+            coords: Vector4::zeros(),
         }
     }
 }
diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs
index cf338a95..d9f33f27 100644
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@@ -25,7 +25,8 @@ use crate::storage::{Storage, StorageMut};
 )]
 #[derive(Clone, Debug)]
 pub struct Cholesky<N: SimdComplexField, D: Dim>
-where DefaultAllocator: Allocator<N, D, D>
+where
+    DefaultAllocator: Allocator<N, D, D>,
 {
     chol: MatrixN<N, D>,
 }
@@ -38,7 +39,8 @@ where
 }
 
 impl<N: SimdComplexField, D: Dim> Cholesky<N, D>
-where DefaultAllocator: Allocator<N, D, D>
+where
+    DefaultAllocator: Allocator<N, D, D>,
 {
     /// Computes the Cholesky decomposition of `matrix` without checking that the matrix is definite-positive.
     ///
@@ -139,7 +141,8 @@ where DefaultAllocator: Allocator<N, D, D>
 }
 
 impl<N: ComplexField, D: Dim> Cholesky<N, D>
-where DefaultAllocator: Allocator<N, D, D>
+where
+    DefaultAllocator: Allocator<N, D, D>,
 {
     /// Attempts to compute the Cholesky decomposition of `matrix`.
     ///
@@ -362,7 +365,8 @@ where DefaultAllocator: Allocator<N, D, D>
 }
 
 impl<N: ComplexField, D: DimSub<Dynamic>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
-where DefaultAllocator: Allocator<N, D, D>
+where
+    DefaultAllocator: Allocator<N, D, D>,
 {
     /// Attempts to compute the Cholesky decomposition of this matrix.
     ///

From 0be9a07f8bd94c21975c54fab4edcbd9ccb9b31d Mon Sep 17 00:00:00 2001
From: sebcrozet <developer@crozet.re>
Date: Sun, 7 Jun 2020 09:28:39 +0200
Subject: [PATCH 19/20] Use the #[rustfmt::skip] attribute instead of
 rustfmt_skip.

---
 nalgebra-glm/src/constructors.rs | 24 +++++++++++++++++-------
 tests/core/edition.rs            | 24 +++++++++++++++---------
 tests/core/matrix_slice.rs       | 27 +++++++++++++++++----------
 tests/core/mod.rs                |  2 +-
 tests/linalg/eigen.rs            |  3 +--
 tests/linalg/full_piv_lu.rs      | 12 +++++-------
 tests/linalg/inverse.rs          |  9 +++++++--
 tests/linalg/lu.rs               |  6 +++---
 tests/linalg/schur.rs            | 11 +++++++----
 tests/linalg/svd.rs              | 19 ++++++++++++++-----
 10 files changed, 87 insertions(+), 50 deletions(-)

diff --git a/nalgebra-glm/src/constructors.rs b/nalgebra-glm/src/constructors.rs
index 949ea9e4..175626a1 100644
--- a/nalgebra-glm/src/constructors.rs
+++ b/nalgebra-glm/src/constructors.rs
@@ -1,9 +1,8 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
-
-use na::{Scalar, RealField, U2, U3, U4};
-use crate::aliases::{TMat, Qua, TVec1, TVec2, TVec3, TVec4, TMat2, TMat2x3, TMat2x4, TMat3, TMat3x2, TMat3x4,
-              TMat4, TMat4x2, TMat4x3};
-
+use crate::aliases::{
+    Qua, TMat, TMat2, TMat2x3, TMat2x4, TMat3, TMat3x2, TMat3x4, TMat4, TMat4x2, TMat4x3, TVec1,
+    TVec2, TVec3, TVec4,
+};
+use na::{RealField, Scalar, U2, U3, U4};
 
 /// Creates a new 1D vector.
 ///
@@ -34,8 +33,8 @@ pub fn vec4<N: Scalar>(x: N, y: N, z: N, w: N) -> TVec4<N> {
     TVec4::new(x, y, z, w)
 }
 
-
 /// Create a new 2x2 matrix.
+#[rustfmt::skip]
 pub fn mat2<N: Scalar>(m11: N, m12: N,
                        m21: N, m22: N) -> TMat2<N> {
     TMat::<N, U2, U2>::new(
@@ -45,6 +44,7 @@ pub fn mat2<N: Scalar>(m11: N, m12: N,
 }
 
 /// Create a new 2x2 matrix.
+#[rustfmt::skip]
 pub fn mat2x2<N: Scalar>(m11: N, m12: N,
                          m21: N, m22: N) -> TMat2<N> {
     TMat::<N, U2, U2>::new(
@@ -54,6 +54,7 @@ pub fn mat2x2<N: Scalar>(m11: N, m12: N,
 }
 
 /// Create a new 2x3 matrix.
+#[rustfmt::skip]
 pub fn mat2x3<N: Scalar>(m11: N, m12: N, m13: N,
                          m21: N, m22: N, m23: N) -> TMat2x3<N> {
     TMat::<N, U2, U3>::new(
@@ -63,6 +64,7 @@ pub fn mat2x3<N: Scalar>(m11: N, m12: N, m13: N,
 }
 
 /// Create a new 2x4 matrix.
+#[rustfmt::skip]
 pub fn mat2x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
                          m21: N, m22: N, m23: N, m24: N) -> TMat2x4<N> {
     TMat::<N, U2, U4>::new(
@@ -72,6 +74,7 @@ pub fn mat2x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
 }
 
 /// Create a new 3x3 matrix.
+#[rustfmt::skip]
 pub fn mat3<N: Scalar>(m11: N, m12: N, m13: N,
                        m21: N, m22: N, m23: N,
                        m31: N, m32: N, m33: N) -> TMat3<N> {
@@ -83,6 +86,7 @@ pub fn mat3<N: Scalar>(m11: N, m12: N, m13: N,
 }
 
 /// Create a new 3x2 matrix.
+#[rustfmt::skip]
 pub fn mat3x2<N: Scalar>(m11: N, m12: N,
                          m21: N, m22: N,
                          m31: N, m32: N) -> TMat3x2<N> {
@@ -94,6 +98,7 @@ pub fn mat3x2<N: Scalar>(m11: N, m12: N,
 }
 
 /// Create a new 3x3 matrix.
+#[rustfmt::skip]
 pub fn mat3x3<N: Scalar>(m11: N, m12: N, m13: N,
                          m21: N, m22: N, m23: N,
                          m31: N, m32: N, m33: N) -> TMat3<N> {
@@ -105,6 +110,7 @@ pub fn mat3x3<N: Scalar>(m11: N, m12: N, m13: N,
 }
 
 /// Create a new 3x4 matrix.
+#[rustfmt::skip]
 pub fn mat3x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
                          m21: N, m22: N, m23: N, m24: N,
                          m31: N, m32: N, m33: N, m34: N) -> TMat3x4<N> {
@@ -116,6 +122,7 @@ pub fn mat3x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
 }
 
 /// Create a new 4x2 matrix.
+#[rustfmt::skip]
 pub fn mat4x2<N: Scalar>(m11: N, m12: N,
                          m21: N, m22: N,
                          m31: N, m32: N,
@@ -129,6 +136,7 @@ pub fn mat4x2<N: Scalar>(m11: N, m12: N,
 }
 
 /// Create a new 4x3 matrix.
+#[rustfmt::skip]
 pub fn mat4x3<N: Scalar>(m11: N, m12: N, m13: N,
                          m21: N, m22: N, m23: N,
                          m31: N, m32: N, m33: N,
@@ -142,6 +150,7 @@ pub fn mat4x3<N: Scalar>(m11: N, m12: N, m13: N,
 }
 
 /// Create a new 4x4 matrix.
+#[rustfmt::skip]
 pub fn mat4x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
                          m21: N, m22: N, m23: N, m24: N,
                          m31: N, m32: N, m33: N, m34: N,
@@ -155,6 +164,7 @@ pub fn mat4x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
 }
 
 /// Create a new 4x4 matrix.
+#[rustfmt::skip]
 pub fn mat4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
                        m21: N, m22: N, m23: N, m24: N,
                        m31: N, m32: N, m33: N, m34: N,
diff --git a/tests/core/edition.rs b/tests/core/edition.rs
index 7c2dee0f..dac25d45 100644
--- a/tests/core/edition.rs
+++ b/tests/core/edition.rs
@@ -1,13 +1,11 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
-
-use na::{Matrix,
-         DMatrix,
-         Matrix3, Matrix4, Matrix5,
-         Matrix4x3, Matrix3x4, Matrix5x3, Matrix3x5, Matrix4x5, Matrix5x4};
+use na::{
+    DMatrix, Matrix, Matrix3, Matrix3x4, Matrix3x5, Matrix4, Matrix4x3, Matrix4x5, Matrix5,
+    Matrix5x3, Matrix5x4,
+};
 use na::{Dynamic, U2, U3, U5};
 
-
 #[test]
+#[rustfmt::skip]
 fn upper_lower_triangular() {
     let m = Matrix4::new(
         11.0, 12.0, 13.0, 14.0,
@@ -173,6 +171,7 @@ fn upper_lower_triangular() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn swap_rows() {
     let mut m = Matrix5x3::new(
         11.0, 12.0, 13.0,
@@ -194,6 +193,7 @@ fn swap_rows() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn swap_columns() {
     let mut m = Matrix3x5::new(
         11.0, 12.0, 13.0, 14.0, 15.0,
@@ -211,6 +211,7 @@ fn swap_columns() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn remove_columns() {
     let m = Matrix3x5::new(
         11, 12, 13, 14, 15,
@@ -261,6 +262,7 @@ fn remove_columns() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn remove_columns_at() {
     let m = DMatrix::from_row_slice(5, 5, &[
         11, 12, 13, 14, 15,
@@ -317,8 +319,8 @@ fn remove_columns_at() {
     assert_eq!(m.remove_columns_at(&[0,3,4]), expected3);
 }
 
-
 #[test]
+#[rustfmt::skip]
 fn remove_rows() {
     let m = Matrix5x3::new(
         11, 12, 13,
@@ -374,6 +376,7 @@ fn remove_rows() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn remove_rows_at() {
     let m = DMatrix::from_row_slice(5, 5, &[
         11, 12, 13, 14, 15,
@@ -424,8 +427,8 @@ fn remove_rows_at() {
     assert_eq!(m.remove_rows_at(&[0,3,4]), expected3);
 }
 
-
 #[test]
+#[rustfmt::skip]
 fn insert_columns() {
     let m = Matrix5x3::new(
         11, 12, 13,
@@ -490,6 +493,7 @@ fn insert_columns() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn insert_columns_to_empty_matrix() {
     let m1 = DMatrix::repeat(0, 0, 0);
     let m2 = DMatrix::repeat(3, 0, 0);
@@ -502,6 +506,7 @@ fn insert_columns_to_empty_matrix() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn insert_rows() {
     let m = Matrix3x5::new(
         11, 12, 13, 14, 15,
@@ -573,6 +578,7 @@ fn insert_rows_to_empty_matrix() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn resize() {
     let m = Matrix3x5::new(
         11, 12, 13, 14, 15,
diff --git a/tests/core/matrix_slice.rs b/tests/core/matrix_slice.rs
index 3a6a1aa5..5fd90b0f 100644
--- a/tests/core/matrix_slice.rs
+++ b/tests/core/matrix_slice.rs
@@ -1,18 +1,15 @@
 #![allow(non_snake_case)]
-#![cfg_attr(rustfmt, rustfmt_skip)]
 
+use na::{
+    DMatrix, DMatrixSlice, DMatrixSliceMut, Matrix2, Matrix2x3, Matrix2x4, Matrix2x6, Matrix3,
+    Matrix3x2, Matrix3x4, Matrix4x2, Matrix6x2, MatrixSlice2, MatrixSlice2x3, MatrixSlice2xX,
+    MatrixSlice3, MatrixSlice3x2, MatrixSliceMut2, MatrixSliceMut2x3, MatrixSliceMut2xX,
+    MatrixSliceMut3, MatrixSliceMut3x2, MatrixSliceMutXx3, MatrixSliceXx3, RowVector4, Vector3,
+};
 use na::{U2, U3, U4};
-use na::{DMatrix,
-         RowVector4,
-         Vector3,
-         Matrix2, Matrix3,
-         Matrix2x3, Matrix3x2, Matrix3x4, Matrix4x2, Matrix2x4, Matrix6x2, Matrix2x6,
-         MatrixSlice2, MatrixSlice3, MatrixSlice2x3, MatrixSlice3x2,
-         MatrixSliceXx3, MatrixSlice2xX, DMatrixSlice,
-         MatrixSliceMut2, MatrixSliceMut3, MatrixSliceMut2x3, MatrixSliceMut3x2,
-         MatrixSliceMutXx3, MatrixSliceMut2xX, DMatrixSliceMut};
 
 #[test]
+#[rustfmt::skip]
 fn nested_fixed_slices() {
     let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0,
                            21.0, 22.0, 23.0, 24.0,
@@ -38,6 +35,7 @@ fn nested_fixed_slices() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn nested_slices() {
     let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0,
                            21.0, 22.0, 23.0, 24.0,
@@ -63,6 +61,7 @@ fn nested_slices() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn slice_mut() {
     let mut a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0,
                                21.0, 22.0, 23.0, 24.0,
@@ -82,6 +81,7 @@ fn slice_mut() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn nested_row_slices() {
     let a = Matrix6x2::new(11.0, 12.0,
                            21.0, 22.0,
@@ -105,6 +105,7 @@ fn nested_row_slices() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn row_slice_mut() {
     let mut a = Matrix6x2::new(11.0, 12.0,
                                21.0, 22.0,
@@ -129,6 +130,7 @@ fn row_slice_mut() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn nested_col_slices() {
     let a = Matrix2x6::new(11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
                            21.0, 22.0, 23.0, 24.0, 25.0, 26.0);
@@ -146,6 +148,7 @@ fn nested_col_slices() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn col_slice_mut() {
     let mut a = Matrix2x6::new(11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
                                21.0, 22.0, 23.0, 24.0, 25.0, 26.0);
@@ -163,6 +166,7 @@ fn col_slice_mut() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn rows_range_pair() {
     let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0,
                            21.0, 22.0, 23.0, 24.0,
@@ -180,6 +184,7 @@ fn rows_range_pair() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn columns_range_pair() {
     let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0,
                            21.0, 22.0, 23.0, 24.0,
@@ -198,6 +203,7 @@ fn columns_range_pair() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn new_slice() {
     let data = [ 1.0, 2.0,  3.0,  4.0,
                  5.0, 6.0,  7.0,  8.0,
@@ -228,6 +234,7 @@ fn new_slice() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn new_slice_mut() {
     let data = [ 1.0, 2.0,  3.0,  4.0,
                  5.0, 6.0,  7.0,  8.0,
diff --git a/tests/core/mod.rs b/tests/core/mod.rs
index 6e897c85..ec1c4e3e 100644
--- a/tests/core/mod.rs
+++ b/tests/core/mod.rs
@@ -3,9 +3,9 @@ mod abomonation;
 mod blas;
 mod conversion;
 mod edition;
+mod empty;
 mod matrix;
 mod matrix_slice;
-mod empty;
 #[cfg(feature = "mint")]
 mod mint;
 mod serde;
diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs
index 8b2f8ed1..1269ed45 100644
--- a/tests/linalg/eigen.rs
+++ b/tests/linalg/eigen.rs
@@ -1,5 +1,3 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
-
 use na::DMatrix;
 
 #[cfg(feature = "arbitrary")]
@@ -67,6 +65,7 @@ mod quickcheck_tests {
 
 // Test proposed on the issue #176 of rulinalg.
 #[test]
+#[rustfmt::skip]
 fn symmetric_eigen_singular_24x24() {
     let m = DMatrix::from_row_slice(
         24,
diff --git a/tests/linalg/full_piv_lu.rs b/tests/linalg/full_piv_lu.rs
index 320fe241..0bb832cd 100644
--- a/tests/linalg/full_piv_lu.rs
+++ b/tests/linalg/full_piv_lu.rs
@@ -1,8 +1,7 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
-
 use na::Matrix3;
 
 #[test]
+#[rustfmt::skip]
 fn full_piv_lu_simple() {
     let m = Matrix3::new(
         2.0, -1.0,  0.0,
@@ -22,11 +21,11 @@ fn full_piv_lu_simple() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn full_piv_lu_simple_with_pivot() {
-    let m = Matrix3::new(
-        0.0, -1.0,  2.0,
-       -1.0,  2.0, -1.0,
-        2.0, -1.0,  0.0);
+    let m = Matrix3::new(0.0, -1.0, 2.0,
+                                        -1.0, 2.0, -1.0,
+                                        2.0, -1.0, 0.0);
 
     let lu = m.full_piv_lu();
     assert_eq!(lu.determinant(), -4.0);
@@ -175,7 +174,6 @@ mod quickcheck_tests {
     gen_tests!(f64, RandScalar<f64>);
 }
 
-
 /*
 #[test]
 fn swap_rows() {
diff --git a/tests/linalg/inverse.rs b/tests/linalg/inverse.rs
index f0be4dd7..ab641a79 100644
--- a/tests/linalg/inverse.rs
+++ b/tests/linalg/inverse.rs
@@ -1,5 +1,3 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
-
 use na::{Matrix1, Matrix2, Matrix3, Matrix4, Matrix5};
 
 #[test]
@@ -11,6 +9,7 @@ fn matrix1_try_inverse() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn matrix2_try_inverse() {
     let a = Matrix2::new(  5.0,  -2.0,
                          -10.0,   1.0);
@@ -23,6 +22,7 @@ fn matrix2_try_inverse() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn matrix3_try_inverse() {
     let a = Matrix3::new(-3.0,  2.0,  0.0,
                          -6.0,  9.0, -2.0,
@@ -37,6 +37,7 @@ fn matrix3_try_inverse() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn matrix4_try_inverse_issue_214() {
     let m1 = Matrix4::new(
         -0.34727043,         0.00000005397217,   -0.000000000000003822135, -0.000000000000003821371,
@@ -58,6 +59,7 @@ fn matrix4_try_inverse_issue_214() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn matrix5_try_inverse() {
     // Dimension 5 is chosen so that the inversion happens by Gaussian elimination.
     // (at the time of writing dimensions <= 3 are implemented as analytic formulas, but we choose
@@ -90,6 +92,7 @@ fn matrix1_try_inverse_scaled_identity() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn matrix2_try_inverse_scaled_identity() {
     // A perfectly invertible matrix with
     // very small coefficients
@@ -103,6 +106,7 @@ fn matrix2_try_inverse_scaled_identity() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn matrix3_try_inverse_scaled_identity() {
     // A perfectly invertible matrix with
     // very small coefficients
@@ -118,6 +122,7 @@ fn matrix3_try_inverse_scaled_identity() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn matrix5_try_inverse_scaled_identity() {
     // A perfectly invertible matrix with
     // very small coefficients
diff --git a/tests/linalg/lu.rs b/tests/linalg/lu.rs
index 69c387d2..7fab6b01 100644
--- a/tests/linalg/lu.rs
+++ b/tests/linalg/lu.rs
@@ -1,8 +1,7 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
-
 use na::Matrix3;
 
 #[test]
+#[rustfmt::skip]
 fn lu_simple() {
     let m = Matrix3::new(
         2.0, -1.0,  0.0,
@@ -21,6 +20,7 @@ fn lu_simple() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn lu_simple_with_pivot() {
     let m = Matrix3::new(
         0.0, -1.0,  2.0,
@@ -41,7 +41,7 @@ fn lu_simple_with_pivot() {
 #[cfg(feature = "arbitrary")]
 mod quickcheck_tests {
     #[allow(unused_imports)]
-    use crate::core::helper::{RandScalar, RandComplex};
+    use crate::core::helper::{RandComplex, RandScalar};
 
     macro_rules! gen_tests(
         ($module: ident, $scalar: ty) => {
diff --git a/tests/linalg/schur.rs b/tests/linalg/schur.rs
index f9c923a2..2086ce2d 100644
--- a/tests/linalg/schur.rs
+++ b/tests/linalg/schur.rs
@@ -1,12 +1,11 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
-
 use na::{DMatrix, Matrix3, Matrix4};
 
 #[test]
+#[rustfmt::skip]
 fn schur_simpl_mat3() {
     let m = Matrix3::new(-2.0, -4.0, 2.0,
-                         -2.0,  1.0, 2.0,
-                          4.0,  2.0, 5.0);
+                                        -2.0,  1.0, 2.0,
+                                        4.0,  2.0, 5.0);
 
     let schur = m.schur();
     let (vecs, vals) = schur.unpack();
@@ -83,6 +82,7 @@ mod quickcheck_tests {
 }
 
 #[test]
+#[rustfmt::skip]
 fn schur_static_mat4_fail() {
     let m = Matrix4::new(
          33.32699857679677,  46.794945978960044, -20.792148817005838,   84.73945485997737,
@@ -95,6 +95,7 @@ fn schur_static_mat4_fail() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn schur_static_mat4_fail2() {
     let m = Matrix4::new(
         14.623586538485966, 7.646156622760756, -52.11923331576265, -97.50030223503413,
@@ -107,6 +108,7 @@ fn schur_static_mat4_fail2() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn schur_static_mat3_fail() {
     let m = Matrix3::new(
         -21.58457553143394,   -67.3881542667948, -14.619829849784338,
@@ -119,6 +121,7 @@ fn schur_static_mat3_fail() {
 
 // Test proposed on the issue #176 of rulinalg.
 #[test]
+#[rustfmt::skip]
 fn schur_singular() {
     let m = DMatrix::from_row_slice(24, 24, &[
         1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  0.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs
index ca7bab4c..cd44b61d 100644
--- a/tests/linalg/svd.rs
+++ b/tests/linalg/svd.rs
@@ -1,4 +1,3 @@
-#![cfg_attr(rustfmt, rustfmt_skip)]
 use na::{DMatrix, Matrix6};
 
 #[cfg(feature = "arbitrary")]
@@ -160,9 +159,9 @@ mod quickcheck_tests {
     gen_tests!(f64, RandScalar<f64>);
 }
 
-
 // Test proposed on the issue #176 of rulinalg.
 #[test]
+#[rustfmt::skip]
 fn svd_singular() {
     let m = DMatrix::from_row_slice(24, 24, &[
         1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  0.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
@@ -202,6 +201,7 @@ fn svd_singular() {
 
 // Same as the previous test but with one additional row.
 #[test]
+#[rustfmt::skip]
 fn svd_singular_vertical() {
     let m = DMatrix::from_row_slice(25, 24, &[
         1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  0.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
@@ -241,6 +241,7 @@ fn svd_singular_vertical() {
 
 // Same as the previous test but with one additional column.
 #[test]
+#[rustfmt::skip]
 fn svd_singular_horizontal() {
     let m = DMatrix::from_row_slice(24, 25, &[
         1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  0.0,  1.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,   0.0,
@@ -299,6 +300,7 @@ fn svd_identity() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn svd_with_delimited_subproblem() {
     let mut m = DMatrix::<f64>::from_element(10, 10, 0.0);
     m[(0,0)] = 1.0;  m[(0,1)] = 2.0;
@@ -334,6 +336,7 @@ fn svd_with_delimited_subproblem() {
 }
 
 #[test]
+#[rustfmt::skip]
 fn svd_fail() {
     let m = Matrix6::new(
         0.9299319121545955,   0.9955870335651049,   0.8824725266413644,  0.28966880207132295,  0.06102723649846409,   0.9311880746048009,
@@ -351,6 +354,12 @@ fn svd_fail() {
 fn svd_err() {
     let m = DMatrix::from_element(10, 10, 0.0);
     let svd = m.clone().svd(false, false);
-    assert_eq!(Err("SVD recomposition: U and V^t have not been computed."), svd.clone().recompose());
-    assert_eq!(Err("SVD pseudo inverse: the epsilon must be non-negative."), svd.clone().pseudo_inverse(-1.0));
-}
\ No newline at end of file
+    assert_eq!(
+        Err("SVD recomposition: U and V^t have not been computed."),
+        svd.clone().recompose()
+    );
+    assert_eq!(
+        Err("SVD pseudo inverse: the epsilon must be non-negative."),
+        svd.clone().pseudo_inverse(-1.0)
+    );
+}

From 2198b0e6b42b218c55a0cc7df26cf4432698023d Mon Sep 17 00:00:00 2001
From: sebcrozet <developer@crozet.re>
Date: Sun, 7 Jun 2020 10:29:10 +0200
Subject: [PATCH 20/20] Release v0.21.1

---
 CHANGELOG.md | 9 +++++++++
 Cargo.toml   | 2 +-
 2 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7487b236..57f18f25 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,15 @@ documented here.
 
 This project adheres to [Semantic Versioning](https://semver.org/).
 
+## [0.22.0] - WIP
+
+### Added
+ * `Cholesky::new_unchecked` which build a Cholesky decomposition without checking that its input is
+ positive-definite. It can be use with SIMD types.
+ * The `Default` trait is now implemented for matrices, and quaternions. They are all filled with zeros,
+ except for `UnitQuaternion` which is initialized with the identity.
+ * Matrix exponential `matrix.exp()`.
+ 
 
 ## [0.21.0]
 In this release, we are no longer relying on traits from the __alga__ crate for our generic code.
diff --git a/Cargo.toml b/Cargo.toml
index 21a365e3..4adf2ab1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name    = "nalgebra"
-version = "0.21.0"
+version = "0.21.1"
 authors = [ "Sébastien Crozet <developer@crozet.re>" ]
 
 description = "Linear algebra library with transformations and statically-sized or dynamically-sized matrices."