diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index 1e402d29..73fe3fa4 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -1026,10 +1026,19 @@ impl<N: Complex, D: Dim, S: StorageMut<N, D, D>> Matrix<N, D, D, S> {
 }
 
 impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
-    /// Creates a square matrix with its diagonal set to `diag` and all other entries set to 0.
+    /// The diagonal of this matrix.
     #[inline]
     pub fn diagonal(&self) -> VectorN<N, D>
         where DefaultAllocator: Allocator<N, D> {
+        self.map_diagonal(|e| e)
+    }
+
+    /// Apply the given function to this matrix's diagonal and returns it.
+    ///
+    /// This is a more efficient version of `self.diagonal().map(f)` since this
+    /// allocates only once.
+    pub fn map_diagonal<N2: Scalar>(&self, mut f: impl FnMut(N) -> N2) -> VectorN<N2, D>
+        where DefaultAllocator: Allocator<N2, D> {
         assert!(
             self.is_square(),
             "Unable to get the diagonal of a non-square matrix."
@@ -1040,7 +1049,7 @@ impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
 
         for i in 0..dim.value() {
             unsafe {
-                *res.vget_unchecked_mut(i) = *self.get_unchecked((i, i));
+                *res.vget_unchecked_mut(i) = f(*self.get_unchecked((i, i)));
             }
         }
 
diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs
index 51087e96..7de0b887 100644
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@@ -272,13 +272,15 @@ where
     }
 
     /// The diagonal part of this decomposed matrix.
-    pub fn diagonal(&self) -> VectorN<N, DimMinimum<R, C>> {
-        self.diagonal.map(|e| N::from_real(e.modulus()))
+    pub fn diagonal(&self) -> VectorN<N::Real, DimMinimum<R, C>>
+        where DefaultAllocator: Allocator<N::Real, DimMinimum<R, C>> {
+        self.diagonal.map(|e| e.modulus())
     }
 
     /// The off-diagonal part of this decomposed matrix.
-    pub fn off_diagonal(&self) -> VectorN<N, DimDiff<DimMinimum<R, C>, U1>> {
-        self.off_diagonal.map(|e| N::from_real(e.modulus()))
+    pub fn off_diagonal(&self) -> VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>
+        where DefaultAllocator: Allocator<N::Real, DimDiff<DimMinimum<R, C>, U1>> {
+        self.off_diagonal.map(|e| e.modulus())
     }
 
     #[doc(hidden)]
diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs
index 72f4a427..5943e0c1 100644
--- a/src/linalg/givens.rs
+++ b/src/linalg/givens.rs
@@ -28,6 +28,16 @@ impl<N: Complex> GivensRotation<N> {
         }
     }
 
+    /// Initializes a Givens rotation from its components.
+    ///
+    /// The components are copies as-is. It is not checked whether they describe
+    /// an actually valid Givens rotation.
+    pub fn new_unchecked(c: N::Real, s: N) -> Self {
+       Self {
+           c, s
+       }
+    }
+
     /// Initializes a Givens rotation from its non-normalized cosine an sine components.
     pub fn new(c: N, s: N) -> (Self, N) {
         Self::try_new(c, s, N::Real::zero()).unwrap()
diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs
index 852878f7..8de75829 100644
--- a/src/linalg/svd.rs
+++ b/src/linalg/svd.rs
@@ -2,11 +2,11 @@
 use serde::{Deserialize, Serialize};
 
 use num_complex::Complex as NumComplex;
-use num::Zero;
+use num::{Zero, One};
 use std::ops::MulAssign;
 use approx::AbsDiffEq;
 
-use alga::general::Complex;
+use alga::general::{Real, Complex};
 use allocator::Allocator;
 use base::{DefaultAllocator, Matrix, Matrix2x3, MatrixMN, Vector2, VectorN};
 use constraint::{SameNumberOfRows, ShapeConstraint};
@@ -23,51 +23,47 @@ use linalg::givens::GivensRotation;
 #[cfg_attr(
     feature = "serde-serialize",
     serde(bound(
-        serialize = "DefaultAllocator: Allocator<N, R, C>    +
-                           Allocator<N, DimMinimum<R, C>>    +
+        serialize = "DefaultAllocator: Allocator<N::Real, DimMinimum<R, C>>    +
                            Allocator<N, DimMinimum<R, C>, C> +
                            Allocator<N, R, DimMinimum<R, C>>,
          MatrixMN<N, R, DimMinimum<R, C>>: Serialize,
          MatrixMN<N, DimMinimum<R, C>, C>: Serialize,
-         VectorN<N, DimMinimum<R, C>>: Serialize"
+         VectorN<N::Real, DimMinimum<R, C>>: Serialize"
     ))
 )]
 #[cfg_attr(
     feature = "serde-serialize",
     serde(bound(
-        deserialize = "DefaultAllocator: Allocator<N, R, C>  +
-                           Allocator<N, DimMinimum<R, C>>    +
+        deserialize = "DefaultAllocator: Allocator<N::Real, DimMinimum<R, C>>    +
                            Allocator<N, DimMinimum<R, C>, C> +
                            Allocator<N, R, DimMinimum<R, C>>,
          MatrixMN<N, R, DimMinimum<R, C>>: Deserialize<'de>,
          MatrixMN<N, DimMinimum<R, C>, C>: Deserialize<'de>,
-         VectorN<N, DimMinimum<R, C>>: Deserialize<'de>"
+         VectorN<N::Real, DimMinimum<R, C>>: Deserialize<'de>"
     ))
 )]
 #[derive(Clone, Debug)]
 pub struct SVD<N: Complex, R: DimMin<C>, C: Dim>
 where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
         + Allocator<N, R, DimMinimum<R, C>>
-        + Allocator<N, DimMinimum<R, C>>
+        + Allocator<N::Real, DimMinimum<R, C>>
 {
     /// The left-singular vectors `U` of this SVD.
     pub u: Option<MatrixMN<N, R, DimMinimum<R, C>>>,
     /// The right-singular vectors `V^t` of this SVD.
     pub v_t: Option<MatrixMN<N, DimMinimum<R, C>, C>>,
-    // FIXME: this should a vector of N::Real
-    // because singular values are necessarily real.
     /// The singular values of this SVD.
-    pub singular_values: VectorN<N, DimMinimum<R, C>>,
+    pub singular_values: VectorN<N::Real, DimMinimum<R, C>>,
 }
 
 impl<N: Complex, R: DimMin<C>, C: Dim> Copy for SVD<N, R, C>
 where
     DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
         + Allocator<N, R, DimMinimum<R, C>>
-        + Allocator<N, DimMinimum<R, C>>,
+        + Allocator<N::Real, DimMinimum<R, C>>,
     MatrixMN<N, R, DimMinimum<R, C>>: Copy,
     MatrixMN<N, DimMinimum<R, C>, C>: Copy,
-    VectorN<N, DimMinimum<R, C>>: Copy,
+    VectorN<N::Real, DimMinimum<R, C>>: Copy,
 {}
 
 impl<N: Complex, R: DimMin<C>, C: Dim> SVD<N, R, C>
@@ -79,7 +75,9 @@ where
         + Allocator<N, DimDiff<DimMinimum<R, C>, U1>>
         + Allocator<N, DimMinimum<R, C>, C>
         + Allocator<N, R, DimMinimum<R, C>>
-        + Allocator<N, DimMinimum<R, C>>,
+        + Allocator<N, DimMinimum<R, C>>
+        + Allocator<N::Real, DimMinimum<R, C>>
+        + Allocator<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
 {
     /// Computes the Singular Value Decomposition of `matrix` using implicit shift.
     pub fn new(matrix: MatrixMN<N, R, C>, compute_u: bool, compute_v: bool) -> Self {
@@ -155,7 +153,7 @@ where
 
                 for k in start..n {
                     let m12 = if k == n - 1 {
-                        N::zero()
+                        N::Real::zero()
                     } else {
                         off_diagonal[k + 1]
                     };
@@ -163,15 +161,15 @@ where
                     let mut subm = Matrix2x3::new(
                         diagonal[k],
                         off_diagonal[k],
-                        N::zero(),
-                        N::zero(),
+                        N::Real::zero(),
+                        N::Real::zero(),
                         diagonal[k + 1],
                         m12,
                     );
 
                     if let Some((rot1, norm1)) = GivensRotation::cancel_y(&vec) {
-                        rot1.inverse()
-                            .rotate_rows(&mut subm.fixed_columns_mut::<U2>(0));
+                        rot1.inverse().rotate_rows(&mut subm.fixed_columns_mut::<U2>(0));
+                        let rot1 = GivensRotation::new_unchecked(rot1.c(), N::from_real(rot1.s()));
 
                         if k > start {
                             // This is not the first iteration.
@@ -182,7 +180,10 @@ where
                         // FIXME: does the case `v.y == 0` ever happen?
                         let (rot2, norm2) =
                             GivensRotation::cancel_y(&v).unwrap_or((GivensRotation::identity(), subm[(0, 0)]));
+
                         rot2.rotate(&mut subm.fixed_columns_mut::<U2>(1));
+                        let rot2 = GivensRotation::new_unchecked(rot2.c(), N::from_real(rot2.s()));
+
                         subm[(0, 0)] = norm2;
 
                         if let Some(ref mut v_t) = v_t {
@@ -219,17 +220,19 @@ where
                 }
             } else if subdim == 2 {
                 // Solve the remaining 2x2 subproblem.
-                let (u2, s, v2) = Self::compute_2x2_uptrig_svd(
+                let (u2, s, v2) = compute_2x2_uptrig_svd(
                     diagonal[start],
                     off_diagonal[start],
                     diagonal[start + 1],
                     compute_u && b.is_upper_diagonal() || compute_v && !b.is_upper_diagonal(),
                     compute_v && b.is_upper_diagonal() || compute_u && !b.is_upper_diagonal(),
                 );
+                let u2 = u2.map(|u2| GivensRotation::new_unchecked(u2.c(), N::from_real(u2.s())));
+                let v2 = v2.map(|v2| GivensRotation::new_unchecked(v2.c(), N::from_real(v2.s())));
 
                 diagonal[start + 0] = s[0];
                 diagonal[start + 1] = s[1];
-                off_diagonal[start] = N::zero();
+                off_diagonal[start] = N::Real::zero();
 
                 if let Some(ref mut u) = u {
                     let rot = if b.is_upper_diagonal() {
@@ -263,17 +266,17 @@ where
             }
         }
 
-        diagonal.scale_mut(m_amax);
+        diagonal *= m_amax;
 
         // Ensure all singular value are non-negative.
         for i in 0..dim {
             let sval = diagonal[i];
-            let (modulus, sign) = sval.to_exp();
-            if modulus != N::Real::zero() {
-                diagonal[i] = N::from_real(modulus);
+
+            if sval < N::Real::zero() {
+                diagonal[i] = -sval;
 
                 if let Some(ref mut u) = u {
-                    u.column_mut(i).mul_assign(sign);
+                    u.column_mut(i).neg_mut();
                 }
             }
         }
@@ -285,55 +288,6 @@ where
         })
     }
 
-    // Explicit formulaes inspired from the paper "Computing the Singular Values of 2-by-2 Complex
-    // Matrices", Sanzheng Qiao and Xiaohong Wang.
-    // http://www.cas.mcmaster.ca/sqrl/papers/sqrl5.pdf
-    fn compute_2x2_uptrig_svd(
-        m11: N,
-        m12: N,
-        m22: N,
-        compute_u: bool,
-        compute_v: bool,
-    ) -> (Option<GivensRotation<N>>, Vector2<N>, Option<GivensRotation<N>>)
-    {
-        let two: N::Real = ::convert(2.0f64);
-        let half: N::Real = ::convert(0.5f64);
-
-        let denom = (m11 + m22).hypot(m12) + (m11 - m22).hypot(m12);
-
-        // NOTE: v1 is the singular value that is the closest to m22.
-        // This prevents cancellation issues when constructing the vector `csv` below. If we chose
-        // otherwise, we would have v1 ~= m11 when m12 is small. This would cause catastrophic
-        // cancellation on `v1 * v1 - m11 * m11` below.
-        let mut v1 = (m11 * m22).scale(two / denom);
-        let mut v2 = N::from_real(half * denom);
-
-        let mut u = None;
-        let mut v_t = None;
-
-        if compute_u || compute_v {
-            let (csv, sgn_v) = GivensRotation::new(m11 * m12, v1 * v1 - m11 * m11);
-            v1 *= sgn_v;
-            v2 *= sgn_v;
-
-            if compute_v {
-                v_t = Some(csv);
-            }
-
-            if compute_u {
-                let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1;
-                let su = (m22 * csv.s()) / v1;
-                let (csu, sgn_u) = GivensRotation::new(cu, su);
-
-                v1 *= sgn_u;
-                v2 *= sgn_u;
-                u = Some(csu);
-            }
-        }
-
-        (u, Vector2::new(v1, v2), v_t)
-    }
-
     /*
     fn display_bidiag(b: &Bidiagonal<N, R, C>, begin: usize, end: usize) {
         for i in begin .. end {
@@ -350,8 +304,8 @@ where
     */
 
     fn delimit_subproblem(
-        diagonal: &mut VectorN<N, DimMinimum<R, C>>,
-        off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
+        diagonal: &mut VectorN<N::Real, DimMinimum<R, C>>,
+        off_diagonal: &mut VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
         u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
         v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
         is_upper_diagonal: bool,
@@ -367,16 +321,16 @@ where
             if off_diagonal[m].is_zero()
                 || off_diagonal[m].modulus() <= eps * (diagonal[n].modulus() + diagonal[m].modulus())
             {
-                off_diagonal[m] = N::zero();
+                off_diagonal[m] = N::Real::zero();
             } else if diagonal[m].modulus() <= eps {
-                diagonal[m] = N::zero();
+                diagonal[m] = N::Real::zero();
                 Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, m + 1);
 
                 if m != 0 {
                     Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m - 1);
                 }
             } else if diagonal[n].modulus() <= eps {
-                diagonal[n] = N::zero();
+                diagonal[n] = N::Real::zero();
                 Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m);
             } else {
                 break;
@@ -395,12 +349,12 @@ where
 
             if off_diagonal[m].modulus() <= eps * (diagonal[new_start].modulus() + diagonal[m].modulus())
             {
-                off_diagonal[m] = N::zero();
+                off_diagonal[m] = N::Real::zero();
                 break;
             }
             // FIXME: write a test that enters this case.
             else if diagonal[m].modulus() <= eps {
-                diagonal[m] = N::zero();
+                diagonal[m] = N::Real::zero();
                 Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, n);
 
                 if m != 0 {
@@ -417,8 +371,8 @@ where
 
     // Cancels the i-th off-diagonal element using givens rotations.
     fn cancel_horizontal_off_diagonal_elt(
-        diagonal: &mut VectorN<N, DimMinimum<R, C>>,
-        off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
+        diagonal: &mut VectorN<N::Real, DimMinimum<R, C>>,
+        off_diagonal: &mut VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
         u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
         v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
         is_upper_diagonal: bool,
@@ -427,10 +381,11 @@ where
     )
     {
         let mut v = Vector2::new(off_diagonal[i], diagonal[i + 1]);
-        off_diagonal[i] = N::zero();
+        off_diagonal[i] = N::Real::zero();
 
         for k in i..end {
             if let Some((rot, norm)) = GivensRotation::cancel_x(&v) {
+                let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
                 diagonal[k + 1] = norm;
 
                 if is_upper_diagonal {
@@ -443,9 +398,9 @@ where
                 }
 
                 if k + 1 != end {
-                    v.x = -rot.s() * off_diagonal[k + 1];
+                    v.x = -rot.s().real() * off_diagonal[k + 1];
                     v.y = diagonal[k + 2];
-                    off_diagonal[k + 1] = off_diagonal[k + 1].scale(rot.c());
+                    off_diagonal[k + 1] *= rot.c();
                 }
             } else {
                 break;
@@ -455,8 +410,8 @@ where
 
     // Cancels the i-th off-diagonal element using givens rotations.
     fn cancel_vertical_off_diagonal_elt(
-        diagonal: &mut VectorN<N, DimMinimum<R, C>>,
-        off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
+        diagonal: &mut VectorN<N::Real, DimMinimum<R, C>>,
+        off_diagonal: &mut VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
         u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
         v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
         is_upper_diagonal: bool,
@@ -464,10 +419,11 @@ where
     )
     {
         let mut v = Vector2::new(diagonal[i], off_diagonal[i]);
-        off_diagonal[i] = N::zero();
+        off_diagonal[i] = N::Real::zero();
 
         for k in (0..i + 1).rev() {
             if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
+                let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
                 diagonal[k] = norm;
 
                 if is_upper_diagonal {
@@ -481,8 +437,8 @@ where
 
                 if k > 0 {
                     v.x = diagonal[k - 1];
-                    v.y = rot.s() * off_diagonal[k - 1];
-                    off_diagonal[k - 1] = off_diagonal[k - 1].scale(rot.c());
+                    v.y = rot.s().real() * off_diagonal[k - 1];
+                    off_diagonal[k - 1] *= rot.c();
                 }
             } else {
                 break;
@@ -497,7 +453,7 @@ where
             eps >= N::Real::zero(),
             "SVD rank: the epsilon must be non-negative."
         );
-        self.singular_values.iter().filter(|e| e.asum() > eps).count()
+        self.singular_values.iter().filter(|e| **e > eps).count()
     }
 
     /// Rebuild the original matrix.
@@ -510,7 +466,7 @@ where
             (Some(mut u), Some(v_t)) => {
                 for i in 0..self.singular_values.len() {
                     let val = self.singular_values[i];
-                    u.column_mut(i).mul_assign(val);
+                    u.column_mut(i).scale_mut(val);
                 }
                 Ok(u * v_t)
             }
@@ -536,10 +492,10 @@ where
             for i in 0..self.singular_values.len() {
                 let val = self.singular_values[i];
 
-                if val.asum() > eps {
-                    self.singular_values[i] = N::one() / val;
+                if val > eps {
+                    self.singular_values[i] = N::Real::one() / val;
                 } else {
-                    self.singular_values[i] = N::zero();
+                    self.singular_values[i] = N::Real::zero();
                 }
             }
 
@@ -575,8 +531,8 @@ where
 
                         for i in 0..self.singular_values.len() {
                             let val = self.singular_values[i];
-                            if val.asum() > eps {
-                                col[i] /= val;
+                            if val > eps {
+                                col[i] = col[i].unscale(val);
                             } else {
                                 col[i] = N::zero();
                             }
@@ -602,7 +558,9 @@ where
         + Allocator<N, DimDiff<DimMinimum<R, C>, U1>>
         + Allocator<N, DimMinimum<R, C>, C>
         + Allocator<N, R, DimMinimum<R, C>>
-        + Allocator<N, DimMinimum<R, C>>,
+        + Allocator<N, DimMinimum<R, C>>
+        + Allocator<N::Real, DimMinimum<R, C>>
+        + Allocator<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
 {
     /// Computes the Singular Value Decomposition using implicit shift.
     pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<N, R, C> {
@@ -631,7 +589,7 @@ where
     }
 
     /// Computes the singular values of this matrix.
-    pub fn singular_values(&self) -> VectorN<N, DimMinimum<R, C>> {
+    pub fn singular_values(&self) -> VectorN<N::Real, DimMinimum<R, C>> {
         SVD::new(self.clone_owned(), false, false).singular_values
     }
 
@@ -653,3 +611,53 @@ where
         SVD::new(self.clone_owned(), true, true).pseudo_inverse(eps)
     }
 }
+
+
+// Explicit formulae inspired from the paper "Computing the Singular Values of 2-by-2 Complex
+// Matrices", Sanzheng Qiao and Xiaohong Wang.
+// http://www.cas.mcmaster.ca/sqrl/papers/sqrl5.pdf
+fn compute_2x2_uptrig_svd<N: Real>(
+    m11: N,
+    m12: N,
+    m22: N,
+    compute_u: bool,
+    compute_v: bool,
+) -> (Option<GivensRotation<N>>, Vector2<N>, Option<GivensRotation<N>>)
+{
+    let two: N::Real = ::convert(2.0f64);
+    let half: N::Real = ::convert(0.5f64);
+
+    let denom = (m11 + m22).hypot(m12) + (m11 - m22).hypot(m12);
+
+    // NOTE: v1 is the singular value that is the closest to m22.
+    // This prevents cancellation issues when constructing the vector `csv` below. If we chose
+    // otherwise, we would have v1 ~= m11 when m12 is small. This would cause catastrophic
+    // cancellation on `v1 * v1 - m11 * m11` below.
+    let mut v1 = m11 * m22 * two / denom;
+    let mut v2 = half * denom;
+
+    let mut u = None;
+    let mut v_t = None;
+
+    if compute_u || compute_v {
+        let (csv, sgn_v) = GivensRotation::new(m11 * m12, v1 * v1 - m11 * m11);
+        v1 *= sgn_v;
+        v2 *= sgn_v;
+
+        if compute_v {
+            v_t = Some(csv);
+        }
+
+        if compute_u {
+            let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1;
+            let su = (m22 * csv.s()) / v1;
+            let (csu, sgn_u) = GivensRotation::new(cu, su);
+
+            v1 *= sgn_u;
+            v2 *= sgn_u;
+            u = Some(csu);
+        }
+    }
+
+    (u, Vector2::new(v1, v2), v_t)
+}
\ No newline at end of file
diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs
index b5437d83..a0f614fe 100644
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@@ -22,8 +22,8 @@ use linalg::SymmetricTridiagonal;
     feature = "serde-serialize",
     serde(bound(
         serialize = "DefaultAllocator: Allocator<N, D, D> +
-                           Allocator<N, D>,
-         VectorN<N, D>: Serialize,
+                           Allocator<N::Real, D>,
+         VectorN<N::Real, D>: Serialize,
          MatrixN<N, D>: Serialize"
     ))
 )]
@@ -31,32 +31,31 @@ use linalg::SymmetricTridiagonal;
     feature = "serde-serialize",
     serde(bound(
         deserialize = "DefaultAllocator: Allocator<N, D, D> +
-                           Allocator<N, D>,
-         VectorN<N, D>: Deserialize<'de>,
+                           Allocator<N::Real, D>,
+         VectorN<N::Real, D>: Deserialize<'de>,
          MatrixN<N, D>: Deserialize<'de>"
     ))
 )]
 #[derive(Clone, Debug)]
 pub struct SymmetricEigen<N: Complex, D: Dim>
-where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
+where DefaultAllocator: Allocator<N, D, D> + Allocator<N::Real, D>
 {
     /// The eigenvectors of the decomposed matrix.
     pub eigenvectors: MatrixN<N, D>,
 
-    // FIXME: this should be a VectorN<N::Real, D>
     /// The unsorted eigenvalues of the decomposed matrix.
-    pub eigenvalues: VectorN<N, D>,
+    pub eigenvalues: VectorN<N::Real, D>,
 }
 
 impl<N: Complex, D: Dim> Copy for SymmetricEigen<N, D>
 where
-    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<N::Real, D>,
     MatrixN<N, D>: Copy,
-    VectorN<N, D>: Copy,
+    VectorN<N::Real, D>: Copy,
 {}
 
 impl<N: Complex, D: Dim> SymmetricEigen<N, D>
-where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
+where DefaultAllocator: Allocator<N, D, D> + Allocator<N::Real, D>
 {
     /// Computes the eigendecomposition of the given symmetric matrix.
     ///
@@ -64,7 +63,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
     pub fn new(m: MatrixN<N, D>) -> Self
     where
         D: DimSub<U1>,
-        DefaultAllocator: Allocator<N, DimDiff<D, U1>>,
+        DefaultAllocator: Allocator<N, DimDiff<D, U1>> + // For tridiagonalization
+        Allocator<N::Real, DimDiff<D, U1>>,
     {
         Self::try_new(m, N::Real::default_epsilon(), 0).unwrap()
     }
@@ -83,7 +83,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
     pub fn try_new(m: MatrixN<N, D>, eps: N::Real, max_niter: usize) -> Option<Self>
     where
         D: DimSub<U1>,
-        DefaultAllocator: Allocator<N, DimDiff<D, U1>>,
+        DefaultAllocator: Allocator<N, DimDiff<D, U1>> + // For tridiagonalization
+                          Allocator<N::Real, DimDiff<D, U1>>,
     {
         Self::do_decompose(m, true, eps, max_niter).map(|(vals, vecs)| SymmetricEigen {
             eigenvectors: vecs.unwrap(),
@@ -96,10 +97,11 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
         eigenvectors: bool,
         eps: N::Real,
         max_niter: usize,
-    ) -> Option<(VectorN<N, D>, Option<MatrixN<N, D>>)>
+    ) -> Option<(VectorN<N::Real, D>, Option<MatrixN<N, D>>)>
     where
         D: DimSub<U1>,
-        DefaultAllocator: Allocator<N, DimDiff<D, U1>>,
+        DefaultAllocator: Allocator<N, DimDiff<D, U1>> + // For tridiagonalization
+                          Allocator<N::Real, DimDiff<D, U1>>,
     {
         assert!(
             m.is_square(),
@@ -158,41 +160,29 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                             off_diag[i - 1] = norm;
                         }
 
+
                         let mii = diag[i];
                         let mjj = diag[j];
                         let mij = off_diag[i];
 
                         let cc = rot.c() * rot.c();
-                        let ss = rot.s().modulus_squared(); // rot.s() * rot.s().conjugate()
-                        let cs = rot.s().scale(rot.c());
+                        let ss = rot.s() * rot.s();
+                        let cs = rot.c() * rot.s();
 
-                        // b = cs * mij.conjugate() + cs.conjugate() * mij
-                        let b = N::from_real((cs * mij.conjugate()).real() * ::convert(2.0));
+                        let b = cs * ::convert(2.0) * mij;
 
-                        diag[i] = (mii.scale(cc) + mjj.scale(ss)) - b;
-                        diag[j] = (mii.scale(ss) + mjj.scale(cc)) + b;
-                        off_diag[i] = cs * (mii - mjj) + mij.scale(cc) - mij.conjugate() * rot.s() * rot.s();
-
-                        let mut mat = Matrix2::new(
-                            mii, mij.conjugate(),
-                            mij, mjj);
-                        println!("The mat before rotate: {:.5}", mat);
-                        println!("The v before rotate: {:.5?}", v);
-                        rot.rotate(&mut mat);
-                        rot.inverse().rotate_rows(&mut mat);
-                        let mut v2 = v.clone();
-                        rot.rotate(&mut v2);
-                        println!("The v: {:.5?}", v2);
-                        println!("The mat: {:.5}", mat);
-                        println!("Its components: {:.5}, {:.5}, {:.5}", diag[i], diag[j], off_diag[i]);
+                        diag[i] = (cc * mii + ss * mjj) - b;
+                        diag[j] = (ss * mii + cc * mjj) + b;
+                        off_diag[i] = cs * (mii - mjj) + mij * (cc - ss);
 
                         if i != n - 1 {
                             v.x = off_diag[i];
                             v.y = -rot.s() * off_diag[i + 1];
-                            off_diag[i + 1] = off_diag[i + 1].scale(rot.c());
+                            off_diag[i + 1] *= rot.c();
                         }
 
                         if let Some(ref mut q) = q {
+                            let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
                             rot.inverse().rotate_rows(&mut q.fixed_columns_mut::<U2>(i));
                         }
                     } else {
@@ -220,6 +210,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
 
                 if let Some(ref mut q) = q {
                     if let Some((rot, _)) = GivensRotation::try_new(basis.x, basis.y, eps) {
+                        let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
                         rot.rotate_rows(&mut q.fixed_columns_mut::<U2>(start));
                     }
                 }
@@ -245,14 +236,14 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
     }
 
     fn delimit_subproblem(
-        diag: &VectorN<N, D>,
-        off_diag: &mut VectorN<N, DimDiff<D, U1>>,
+        diag: &VectorN<N::Real, D>,
+        off_diag: &mut VectorN<N::Real, DimDiff<D, U1>>,
         end: usize,
         eps: N::Real,
     ) -> (usize, usize)
     where
         D: DimSub<U1>,
-        DefaultAllocator: Allocator<N, DimDiff<D, U1>>,
+        DefaultAllocator: Allocator<N::Real, DimDiff<D, U1>>,
     {
         let mut n = end;
 
@@ -277,7 +268,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
             if off_diag[m].is_zero()
                 || off_diag[m].modulus() <= eps * (diag[new_start].modulus() + diag[m].modulus())
             {
-                off_diag[m] = N::zero();
+                off_diag[m] = N::Real::zero();
                 break;
             }
 
@@ -294,7 +285,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
         let mut u_t = self.eigenvectors.clone();
         for i in 0..self.eigenvalues.len() {
             let val = self.eigenvalues[i];
-            u_t.column_mut(i).mul_assign(val);
+            u_t.column_mut(i).scale_mut(val);
         }
         u_t.conjugate_transpose_mut();
         &self.eigenvectors * u_t
@@ -324,7 +315,8 @@ pub fn wilkinson_shift<N: Complex>(tmm: N, tnn: N, tmn: N) -> N {
  *
  */
 impl<N: Complex, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
-where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimDiff<D, U1>>
+where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>> +
+                        Allocator<N::Real, D> + Allocator<N::Real, DimDiff<D, U1>>
 {
     /// Computes the eigendecomposition of this symmetric matrix.
     ///
@@ -351,7 +343,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimD
     /// Computes the eigenvalues of this symmetric matrix.
     ///
     /// Only the lower-triangular part of the matrix is read.
-    pub fn symmetric_eigenvalues(&self) -> VectorN<N, D> {
+    pub fn symmetric_eigenvalues(&self) -> VectorN<N::Real, D> {
         SymmetricEigen::do_decompose(self.clone_owned(), false, N::Real::default_epsilon(), 0)
             .unwrap()
             .0
diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs
index 8fd36003..0b8f5675 100644
--- a/src/linalg/symmetric_tridiagonal.rs
+++ b/src/linalg/symmetric_tridiagonal.rs
@@ -102,30 +102,30 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
 
     /// Retrieve the orthogonal transformation, diagonal, and off diagonal elements of this
     /// decomposition.
-    pub fn unpack(self) -> (MatrixN<N, D>, VectorN<N, D>, VectorN<N, DimDiff<D, U1>>)
-    where DefaultAllocator: Allocator<N, D> {
+    pub fn unpack(self) -> (MatrixN<N, D>, VectorN<N::Real, D>, VectorN<N::Real, DimDiff<D, U1>>)
+    where DefaultAllocator: Allocator<N::Real, D>
+                          + Allocator<N::Real, DimDiff<D, U1>> {
         let diag = self.diagonal();
         let q = self.q();
 
-        (q, diag, self.off_diagonal.apply_into(|e| N::from_real(e.modulus())))
+        (q, diag, self.off_diagonal.map(N::modulus))
     }
 
     /// Retrieve the diagonal, and off diagonal elements of this decomposition.
-    pub fn unpack_tridiagonal(mut self) -> (VectorN<N, D>, VectorN<N, DimDiff<D, U1>>)
-    where DefaultAllocator: Allocator<N, D> {
-        let diag = self.diagonal();
-
-        (diag, self.off_diagonal.apply_into(|e| N::from_real(e.modulus())))
+    pub fn unpack_tridiagonal(mut self) -> (VectorN<N::Real, D>, VectorN<N::Real, DimDiff<D, U1>>)
+        where DefaultAllocator: Allocator<N::Real, D>
+                              + Allocator<N::Real, DimDiff<D, U1>> {
+        (self.diagonal(), self.off_diagonal.map(N::modulus))
     }
 
     /// The diagonal components of this decomposition.
-    pub fn diagonal(&self) -> VectorN<N, D>
-    where DefaultAllocator: Allocator<N, D> { self.tri.diagonal() }
+    pub fn diagonal(&self) -> VectorN<N::Real, D>
+    where DefaultAllocator: Allocator<N::Real, D> { self.tri.map_diagonal(|e| e.real()) }
 
     /// The off-diagonal components of this decomposition.
-    pub fn off_diagonal(&self) -> VectorN<N, DimDiff<D, U1>>
-    where DefaultAllocator: Allocator<N, D> {
-        self.off_diagonal.map(|e| N::from_real(e.modulus()))
+    pub fn off_diagonal(&self) -> VectorN<N::Real, DimDiff<D, U1>>
+    where DefaultAllocator: Allocator<N::Real, DimDiff<D, U1>> {
+        self.off_diagonal.map(N::modulus)
     }
 
     /// Computes the orthogonal matrix `Q` of this decomposition.
diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs
index 860c7ce2..8e98b3f7 100644
--- a/tests/linalg/svd.rs
+++ b/tests/linalg/svd.rs
@@ -18,11 +18,11 @@ mod quickcheck_tests {
                 let svd = m.clone().svd(true, true);
                 let recomp_m = svd.clone().recompose().unwrap();
                 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-                let ds = DMatrix::from_diagonal(&s);
+                let ds = DMatrix::from_diagonal(&s.map(|e| Complex::from_real(e)));
 
                 println!("{}{}", &m, &u * &ds * &v_t);
 
-                s.iter().all(|e| e.real() >= 0.0) &&
+                s.iter().all(|e| *e >= 0.0) &&
                 relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) &&
                 relative_eq!(m, recomp_m, epsilon = 1.0e-5)
             }
@@ -35,9 +35,9 @@ mod quickcheck_tests {
             let m = m.map(|e| e.0);
             let svd = m.svd(true, true);
             let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-            let ds = Matrix3::from_diagonal(&s);
+            let ds = Matrix3::from_diagonal(&s.map(|e| Complex::from_real(e)));
 
-            s.iter().all(|e| e.real() >= 0.0) &&
+            s.iter().all(|e| *e >= 0.0) &&
             relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
             u.is_orthogonal(1.0e-5) &&
             v_t.is_orthogonal(1.0e-5)
@@ -47,9 +47,9 @@ mod quickcheck_tests {
             let m = m.map(|e| e.0);
             let svd = m.svd(true, true);
             let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-            let ds = Matrix2::from_diagonal(&s);
+            let ds = Matrix2::from_diagonal(&s.map(|e| Complex::from_real(e)));
 
-            s.iter().all(|e| e.real() >= 0.0) &&
+            s.iter().all(|e| *e >= 0.0) &&
             relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
             u.is_orthogonal(1.0e-5) &&
             v_t.is_orthogonal(1.0e-5)
@@ -60,9 +60,9 @@ mod quickcheck_tests {
             let svd = m.svd(true, true);
             let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
 
-            let ds = Matrix3::from_diagonal(&s);
+            let ds = Matrix3::from_diagonal(&s.map(|e| Complex::from_real(e)));
 
-            s.iter().all(|e| e.real() >= 0.0) &&
+            s.iter().all(|e| *e >= 0.0) &&
             relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
         }
 
@@ -70,9 +70,9 @@ mod quickcheck_tests {
             let m = m.map(|e| e.0);
             let svd = m.svd(true, true);
             let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-            let ds = Matrix2::from_diagonal(&s);
+            let ds = Matrix2::from_diagonal(&s.map(|e| Complex::from_real(e)));
 
-            s.iter().all(|e| e.real() >= 0.0) &&
+            s.iter().all(|e| *e >= 0.0) &&
             relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
         }
 
@@ -80,9 +80,9 @@ mod quickcheck_tests {
             let m = m.map(|e| e.0);
             let svd = m.svd(true, true);
             let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-            let ds = Matrix4::from_diagonal(&s);
+            let ds = Matrix4::from_diagonal(&s.map(|e| Complex::from_real(e)));
 
-            s.iter().all(|e| e.real() >= 0.0) &&
+            s.iter().all(|e| *e >= 0.0) &&
             relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
             u.is_orthogonal(1.0e-5) &&
             v_t.is_orthogonal(1.0e-5)
@@ -92,14 +92,14 @@ mod quickcheck_tests {
             let m = m.map(|e| e.0);
             let svd = m.svd(true, true);
             let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-            let ds = Matrix2::from_diagonal(&s);
+            let ds = Matrix2::from_diagonal(&s.map(|e| Complex::from_real(e)));
 
             println!("u, s, v_t: {}{}{}", u, s, v_t);
             println!("m: {}", m);
             println!("recomp: {}", u * ds * v_t);
             println!("uu_t, vv_t: {}{}", u * u.conjugate_transpose(), v_t.conjugate_transpose() * v_t);
 
-            s.iter().all(|e| e.re >= 0.0) &&
+            s.iter().all(|e| *e >= 0.0) &&
             relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
             u.is_orthogonal(1.0e-5) &&
             v_t.is_orthogonal(1.0e-5)