diff --git a/src/geometry/reflection.rs b/src/geometry/reflection.rs
index fca0f548..237054b6 100644
--- a/src/geometry/reflection.rs
+++ b/src/geometry/reflection.rs
@@ -61,10 +61,10 @@ impl<N: Complex, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
         }
     }
 
-    /// Applies the reflection to the rows of `rhs`.
+    /// Applies the reflection to the rows of `lhs`.
     pub fn reflect_rows<R2: Dim, C2: Dim, S2, S3>(
         &self,
-        rhs: &mut Matrix<N, R2, C2, S2>,
+        lhs: &mut Matrix<N, R2, C2, S2>,
         work: &mut Vector<N, R2, S3>,
     ) where
         S2: StorageMut<N, R2, C2>,
@@ -72,13 +72,13 @@ impl<N: Complex, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
         ShapeConstraint: DimEq<C2, D> + AreMultipliable<R2, C2, D, U1>,
         DefaultAllocator: Allocator<N, D>
     {
-        rhs.mul_to(&self.axis, work);
+        lhs.mul_to(&self.axis, work);
 
         if !self.bias.is_zero() {
             work.add_scalar_mut(-self.bias);
         }
 
         let m_two: N = ::convert(-2.0f64);
-        rhs.ger(m_two, &work, &self.axis.conjugate(), N::one());
+        lhs.ger(m_two, &work, &self.axis.conjugate(), N::one());
     }
 }
diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs
index a44be398..e487758c 100644
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@@ -49,9 +49,9 @@ where
     // contiguous. This prevents some useless copies.
     uv: MatrixMN<N, R, C>,
     /// The diagonal elements of the decomposed matrix.
-    pub diagonal: VectorN<N, DimMinimum<R, C>>,
+    diagonal: VectorN<N, DimMinimum<R, C>>,
     /// The off-diagonal elements of the decomposed matrix.
-    pub off_diagonal: VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
+    off_diagonal: VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
     upper_diagonal: bool,
 }
 
@@ -143,9 +143,9 @@ where
 
         Bidiagonal {
             uv: matrix,
-            diagonal: diagonal,
-            off_diagonal: off_diagonal,
-            upper_diagonal: upper_diagonal,
+            diagonal,
+            off_diagonal,
+            upper_diagonal,
         }
     }
 
@@ -231,7 +231,7 @@ where
         res
     }
 
-    /// Computes the orthogonal matrix `V` of this `U * D * V` decomposition.
+    /// Computes the orthogonal matrix `V_t` of this `U * D * V_t` decomposition.
     pub fn v_t(&self) -> MatrixMN<N, DimMinimum<R, C>, C>
     where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C> {
         let (nrows, ncols) = self.uv.data.shape();
@@ -258,13 +258,13 @@ where
     }
 
     /// The diagonal part of this decomposed matrix.
-    pub fn diagonal(&self) -> &VectorN<N, DimMinimum<R, C>> {
-        &self.diagonal
+    pub fn diagonal(&self) -> VectorN<N, DimMinimum<R, C>> {
+        self.diagonal.map(|e| N::from_real(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
+    pub fn off_diagonal(&self) -> VectorN<N, DimDiff<DimMinimum<R, C>, U1>> {
+        self.off_diagonal.map(|e| N::from_real(e.modulus()))
     }
 
     #[doc(hidden)]
diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs
index 472f3bb3..72f4a427 100644
--- a/src/linalg/givens.rs
+++ b/src/linalg/givens.rs
@@ -1,7 +1,7 @@
 //! Construction of givens rotations.
 
 use alga::general::{Complex, Real};
-use num::Zero;
+use num::{Zero, One};
 use num_complex::Complex as NumComplex;
 
 use base::dimension::{Dim, U2};
@@ -12,60 +12,37 @@ use base::{Vector, Matrix};
 use geometry::UnitComplex;
 
 /// A Givens rotation.
-#[derive(Debug)]
+#[derive(Debug, Clone, Copy)]
 pub struct GivensRotation<N: Complex> {
-    // FIXME: c should be a `N::Real`.
-    c: N,
+    c: N::Real,
     s: N
 }
 
-// XXX: remove this
-/// Computes the rotation `R` required such that the `y` component of `R * v` is zero.
-///
-/// Returns `None` if no rotation is needed (i.e. if `v.y == 0`). Otherwise, this returns the norm
-/// of `v` and the rotation `r` such that `R * v = [ |v|, 0.0 ]^t` where `|v|` is the norm of `v`.
-pub fn cancel_y<N: Real, S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(UnitComplex<N>, N)> {
-    if !v[1].is_zero() {
-        let c = NumComplex::new(v[0], -v[1]);
-        Some(UnitComplex::from_complex_and_get(c))
-    } else {
-        None
-    }
-}
-
-// XXX: remove this
-/// Computes the rotation `R` required such that the `x` component of `R * v` is zero.
-///
-/// Returns `None` if no rotation is needed (i.e. if `v.x == 0`). Otherwise, this returns the norm
-/// of `v` and the rotation `r` such that `R * v = [ 0.0, |v| ]^t` where `|v|` is the norm of `v`.
-pub fn cancel_x<N: Real, S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(UnitComplex<N>, N)> {
-    if !v[0].is_zero() {
-        let c = NumComplex::new(v[1], v[0]);
-        Some(UnitComplex::from_complex_and_get(c))
-    } else {
-        None
-    }
-}
-
-
 // Matrix = UnitComplex * Matrix
 impl<N: Complex> GivensRotation<N> {
+    /// The Givents rotation that does nothing.
+    pub fn identity() -> Self {
+        Self {
+            c: N::Real::one(),
+            s: N::zero()
+        }
+    }
+
     /// Initializes a Givens rotation from its non-normalized cosine an sine components.
-    pub fn new(c: N, s: N) -> Self {
-        let res = Self::try_new(c, s, N::Real::zero()).unwrap();
-        println!("The rot: {:?}", res);
-        res
+    pub fn new(c: N, s: N) -> (Self, N) {
+        Self::try_new(c, s, N::Real::zero()).unwrap()
     }
 
     /// Initializes a Givens rotation form its non-normalized cosine an sine components.
-    pub fn try_new(c: N, s: N, eps: N::Real) -> Option<Self> {
+    pub fn try_new(c: N, s: N, eps: N::Real) -> Option<(Self, N)> {
         let (mod0, sign0) = c.to_exp();
         let denom = (mod0 * mod0 + s.modulus_squared()).sqrt();
 
         if denom > eps {
-            let c = N::from_real(mod0 / denom);
-            let s = s / sign0.scale(denom);
-            Some(Self { c, s })
+            let norm = sign0.scale(denom);
+            let c = mod0 / denom;
+            let s = s / norm;
+            Some((Self { c, s }, norm))
         } else {
             None
         }
@@ -79,8 +56,8 @@ impl<N: Complex> GivensRotation<N> {
         if !v[1].is_zero() {
             let (mod0, sign0) = v[0].to_exp();
             let denom = (mod0 * mod0 + v[1].modulus_squared()).sqrt();
-            let c = N::from_real(mod0 / denom);
-            let s = (sign0 * v[1].conjugate()).unscale(-denom);
+            let c = mod0 / denom;
+            let s = -v[1] / sign0.scale(denom);
             let r = sign0.scale(denom);
             Some((Self { c, s }, r))
         } else {
@@ -94,11 +71,11 @@ impl<N: Complex> GivensRotation<N> {
     /// of `v` and the rotation `r` such that `R * v = [ 0.0, |v| ]^t` where `|v|` is the norm of `v`.
     pub fn cancel_x<S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(Self, N)> {
         if !v[0].is_zero() {
-            let (mod0, sign0) = v[0].to_exp();
-            let denom = (mod0 * mod0 + v[1].modulus_squared()).sqrt();
-            let c = N::from_real(mod0 / denom);
-            let s = (sign0 * v[1].conjugate()).unscale(denom);
-            let r = sign0.scale(denom);
+            let (mod1, sign1) = v[1].to_exp();
+            let denom = (mod1 * mod1 + v[0].modulus_squared()).sqrt();
+            let c = mod1 / denom;
+            let s = (v[0].conjugate() * sign1).unscale(denom);
+            let r = sign1.scale(denom);
             Some((Self { c, s }, r))
         } else {
             None
@@ -106,7 +83,7 @@ impl<N: Complex> GivensRotation<N> {
     }
 
     /// The cos part of this roration.
-    pub fn c(&self) -> N {
+    pub fn c(&self) -> N::Real {
         self.c
     }
 
@@ -140,8 +117,8 @@ impl<N: Complex> GivensRotation<N> {
                 let a = *rhs.get_unchecked((0, j));
                 let b = *rhs.get_unchecked((1, j));
 
-                *rhs.get_unchecked_mut((0, j)) = c * a + -s.conjugate() * b;
-                *rhs.get_unchecked_mut((1, j)) = s * a + c * b;
+                *rhs.get_unchecked_mut((0, j)) = a.scale(c) - s.conjugate() * b;
+                *rhs.get_unchecked_mut((1, j)) = s * a + b.scale(c);
             }
         }
     }
@@ -167,8 +144,8 @@ impl<N: Complex> GivensRotation<N> {
                 let a = *lhs.get_unchecked((j, 0));
                 let b = *lhs.get_unchecked((j, 1));
 
-                *lhs.get_unchecked_mut((j, 0)) = c * a + s * b;
-                *lhs.get_unchecked_mut((j, 1)) = -s.conjugate() * a + c * b;
+                *lhs.get_unchecked_mut((j, 0)) = a.scale(c) + s * b;
+                *lhs.get_unchecked_mut((j, 1)) = -s.conjugate() * a + b.scale(c);
             }
         }
     }
diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs
index 0fe46499..dc97b9b3 100644
--- a/src/linalg/householder.rs
+++ b/src/linalg/householder.rs
@@ -23,20 +23,20 @@ pub fn reflection_axis_mut<N: Complex, D: Dim, S: StorageMut<N, D>>(
     let reflection_norm = reflection_sq_norm.sqrt();
 
     let factor;
-    let scaled_norm;
+    let signed_norm;
 
     unsafe {
-        let (modulus, exp) = column.vget_unchecked(0).to_exp();
-        scaled_norm = exp.scale(reflection_norm);
+        let (modulus, sign) = column.vget_unchecked(0).to_exp();
+        signed_norm = sign.scale(reflection_norm);
         factor = (reflection_sq_norm + modulus * reflection_norm) * ::convert(2.0);
-        *column.vget_unchecked_mut(0) += scaled_norm;
+        *column.vget_unchecked_mut(0) += signed_norm;
     };
 
     if !factor.is_zero() {
         column.unscale_mut(factor.sqrt());
-        (-scaled_norm, true)
+        (-signed_norm, true)
     } else {
-        (-scaled_norm, false)
+        (-signed_norm, false)
     }
 }
 
@@ -67,7 +67,7 @@ pub fn clear_column_unchecked<N: Complex, R: Dim, C: Dim>(
     }
 }
 
-/// Uses an hoseholder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
+/// Uses an householder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
 /// superdiagonal element.
 #[doc(hidden)]
 pub fn clear_row_unchecked<N: Complex, R: Dim, C: Dim>(
@@ -85,6 +85,7 @@ pub fn clear_row_unchecked<N: Complex, R: Dim, C: Dim>(
     axis.tr_copy_from(&top.columns_range(irow + shift..));
 
     let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
+    axis.conjugate_mut(); // So that reflect_rows actually cancels the first row.
     *diag_elt = reflection_norm;
 
     if not_zero {
diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs
index 4418b283..6c7f7194 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -30,6 +30,6 @@ pub use self::lu::*;
 pub use self::permutation_sequence::*;
 pub use self::qr::*;
 pub use self::schur::*;
-pub use self::svd::*;
+//pub use self::svd::*;
 pub use self::symmetric_eigen::*;
 pub use self::symmetric_tridiagonal::*;
diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs
index 2eca9843..b815131e 100644
--- a/src/linalg/schur.rs
+++ b/src/linalg/schur.rs
@@ -417,10 +417,11 @@ where
             if compute_q {
                 // XXX: we have to build the matrix manually because
                 // rot.to_rotation_matrix().unwrap() causes an ICE.
+                let c = N::from_real(rot.c());
                 q = Some(MatrixN::from_column_slice_generic(
                     dim,
                     dim,
-                    &[rot.c(), rot.s(), -rot.s().conjugate(), rot.c()],
+                    &[c, rot.s(), -rot.s().conjugate(), c],
                 ));
             }
         }
@@ -479,9 +480,9 @@ fn compute_2x2_basis<N: Complex, S: Storage<N, U2, U2>>(
         // This is necessary for numerical stability of the normalization of the complex
         // number.
         if x1.modulus() > x2.modulus() {
-            Some(GivensRotation::new(x1, h10))
+            Some(GivensRotation::new(x1, h10).0)
         } else {
-            Some(GivensRotation::new(x2, h10))
+            Some(GivensRotation::new(x2, h10).0)
         }
     } else {
         None
diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs
index 67b49604..bcc93333 100644
--- a/src/linalg/svd.rs
+++ b/src/linalg/svd.rs
@@ -1,27 +1,29 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};
 
-use num_complex::Complex;
+use num_complex::Complex as NumComplex;
+use num::Zero;
 use std::ops::MulAssign;
+use approx::AbsDiffEq;
 
-use alga::general::Real;
+use alga::general::Complex;
 use allocator::Allocator;
 use base::{DefaultAllocator, Matrix, Matrix2x3, MatrixMN, Vector2, VectorN};
 use constraint::{SameNumberOfRows, ShapeConstraint};
 use dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1, U2};
 use storage::Storage;
 
-use geometry::UnitComplex;
 use linalg::givens;
 use linalg::symmetric_eigen;
 use linalg::Bidiagonal;
+use linalg::givens::GivensRotation;
 
 /// Singular Value Decomposition of a general matrix.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
 #[cfg_attr(
     feature = "serde-serialize",
     serde(bound(
-        serialize = "DefaultAllocator: Allocator<N, R, C>                +
+        serialize = "DefaultAllocator: Allocator<N, R, C>    +
                            Allocator<N, DimMinimum<R, C>>    +
                            Allocator<N, DimMinimum<R, C>, C> +
                            Allocator<N, R, DimMinimum<R, C>>,
@@ -33,7 +35,7 @@ use linalg::Bidiagonal;
 #[cfg_attr(
     feature = "serde-serialize",
     serde(bound(
-        deserialize = "DefaultAllocator: Allocator<N, R, C>                +
+        deserialize = "DefaultAllocator: Allocator<N, R, C>  +
                            Allocator<N, DimMinimum<R, C>>    +
                            Allocator<N, DimMinimum<R, C>, C> +
                            Allocator<N, R, DimMinimum<R, C>>,
@@ -43,7 +45,7 @@ use linalg::Bidiagonal;
     ))
 )]
 #[derive(Clone, Debug)]
-pub struct SVD<N: Real, R: DimMin<C>, C: Dim>
+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>>
@@ -52,11 +54,13 @@ where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
     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>>,
 }
 
-impl<N: Real, R: DimMin<C>, C: Dim> Copy for SVD<N, 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>>
@@ -66,7 +70,7 @@ where
     VectorN<N, DimMinimum<R, C>>: Copy,
 {}
 
-impl<N: Real, R: DimMin<C>, C: Dim> SVD<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> SVD<N, R, C>
 where
     DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
     DefaultAllocator: Allocator<N, R, C>
@@ -79,7 +83,7 @@ where
 {
     /// 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 {
-        Self::try_new(matrix, compute_u, compute_v, N::default_epsilon(), 0).unwrap()
+        Self::try_new(matrix, compute_u, compute_v, N::Real::default_epsilon(), 0).unwrap()
     }
 
     /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
@@ -96,7 +100,7 @@ where
         mut matrix: MatrixMN<N, R, C>,
         compute_u: bool,
         compute_v: bool,
-        eps: N,
+        eps: N::Real,
         max_niter: usize,
     ) -> Option<Self>
     {
@@ -108,18 +112,20 @@ where
         let min_nrows_ncols = nrows.min(ncols);
         let dim = min_nrows_ncols.value();
 
-        let m_amax = matrix.amax();
+        let m_amax = matrix.camax();
 
         if !m_amax.is_zero() {
-            matrix /= m_amax;
+            matrix.unscale_mut(m_amax);
         }
 
         let mut b = Bidiagonal::new(matrix);
         let mut u = if compute_u { Some(b.u()) } else { None };
         let mut v_t = if compute_v { Some(b.v_t()) } else { None };
+        let mut diagonal = b.diagonal();
+        let mut off_diagonal = b.off_diagonal();
 
         let mut niter = 0;
-        let (mut start, mut end) = Self::delimit_subproblem(&mut b, &mut u, &mut v_t, dim - 1, eps);
+        let (mut start, mut end) = Self::delimit_subproblem(&mut diagonal, &mut off_diagonal, &mut u, &mut v_t, b.is_upper_diagonal(), dim - 1, eps);
 
         while end != start {
             let subdim = end - start + 1;
@@ -131,19 +137,19 @@ where
 
                 let mut vec;
                 {
-                    let dm = b.diagonal[m];
-                    let dn = b.diagonal[n];
-                    let fm = b.off_diagonal[m];
+                    let dm = diagonal[m];
+                    let dn = diagonal[n];
+                    let fm = off_diagonal[m];
 
-                    let tmm = dm * dm + b.off_diagonal[m - 1] * b.off_diagonal[m - 1];
+                    let tmm = dm * dm + off_diagonal[m - 1] * off_diagonal[m - 1];
                     let tmn = dm * fm;
                     let tnn = dn * dn + fm * fm;
 
                     let shift = symmetric_eigen::wilkinson_shift(tmm, tnn, tmn);
 
                     vec = Vector2::new(
-                        b.diagonal[start] * b.diagonal[start] - shift,
-                        b.diagonal[start] * b.off_diagonal[start],
+                        diagonal[start] * diagonal[start] - shift,
+                        diagonal[start] * off_diagonal[start],
                     );
                 }
 
@@ -151,31 +157,31 @@ where
                     let m12 = if k == n - 1 {
                         N::zero()
                     } else {
-                        b.off_diagonal[k + 1]
+                        off_diagonal[k + 1]
                     };
 
                     let mut subm = Matrix2x3::new(
-                        b.diagonal[k],
-                        b.off_diagonal[k],
+                        diagonal[k],
+                        off_diagonal[k],
                         N::zero(),
                         N::zero(),
-                        b.diagonal[k + 1],
+                        diagonal[k + 1],
                         m12,
                     );
 
-                    if let Some((rot1, norm1)) = givens::cancel_y(&vec) {
-                        rot1.conjugate()
+                    if let Some((rot1, norm1)) = GivensRotation::cancel_y(&vec) {
+                        rot1.inverse()
                             .rotate_rows(&mut subm.fixed_columns_mut::<U2>(0));
 
                         if k > start {
                             // This is not the first iteration.
-                            b.off_diagonal[k - 1] = norm1;
+                            off_diagonal[k - 1] = norm1;
                         }
 
                         let v = Vector2::new(subm[(0, 0)], subm[(1, 0)]);
                         // FIXME: does the case `v.y == 0` ever happen?
                         let (rot2, norm2) =
-                            givens::cancel_y(&v).unwrap_or((UnitComplex::identity(), subm[(0, 0)]));
+                            GivensRotation::cancel_y(&v).unwrap_or((GivensRotation::identity(), subm[(0, 0)]));
                         rot2.rotate(&mut subm.fixed_columns_mut::<U2>(1));
                         subm[(0, 0)] = norm2;
 
@@ -197,12 +203,12 @@ where
                             }
                         }
 
-                        b.diagonal[k + 0] = subm[(0, 0)];
-                        b.diagonal[k + 1] = subm[(1, 1)];
-                        b.off_diagonal[k + 0] = subm[(0, 1)];
+                        diagonal[k + 0] = subm[(0, 0)];
+                        diagonal[k + 1] = subm[(1, 1)];
+                        off_diagonal[k + 0] = subm[(0, 1)];
 
                         if k != n - 1 {
-                            b.off_diagonal[k + 1] = subm[(1, 2)];
+                            off_diagonal[k + 1] = subm[(1, 2)];
                         }
 
                         vec.x = subm[(0, 1)];
@@ -214,16 +220,16 @@ where
             } else if subdim == 2 {
                 // Solve the remaining 2x2 subproblem.
                 let (u2, s, v2) = Self::compute_2x2_uptrig_svd(
-                    b.diagonal[start],
-                    b.off_diagonal[start],
-                    b.diagonal[start + 1],
+                    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(),
                 );
 
-                b.diagonal[start + 0] = s[0];
-                b.diagonal[start + 1] = s[1];
-                b.off_diagonal[start] = N::zero();
+                diagonal[start + 0] = s[0];
+                diagonal[start + 1] = s[1];
+                off_diagonal[start] = N::zero();
 
                 if let Some(ref mut u) = u {
                     let rot = if b.is_upper_diagonal() {
@@ -247,7 +253,7 @@ where
             }
 
             // Re-delimit the subproblem in case some decoupling occurred.
-            let sub = Self::delimit_subproblem(&mut b, &mut u, &mut v_t, end, eps);
+            let sub = Self::delimit_subproblem(&mut diagonal, &mut off_diagonal, &mut u, &mut v_t, b.is_upper_diagonal(), end, eps);
             start = sub.0;
             end = sub.1;
 
@@ -257,24 +263,25 @@ where
             }
         }
 
-        b.diagonal *= m_amax;
+        diagonal.scale_mut(m_amax);
 
         // Ensure all singular value are non-negative.
         for i in 0..dim {
-            let sval = b.diagonal[i];
-            if sval < N::zero() {
-                b.diagonal[i] = -sval;
+            let sval = diagonal[i];
+            let (modulus, sign) = sval.to_exp();
+            if modulus != N::Real::zero() {
+                diagonal[i] = N::from_real(modulus);
 
                 if let Some(ref mut u) = u {
-                    u.column_mut(i).neg_mut();
+                    u.column_mut(i).mul_assign(sign);
                 }
             }
         }
 
         Some(Self {
-            u: u,
-            v_t: v_t,
-            singular_values: b.diagonal,
+            u,
+            v_t,
+            singular_values: diagonal,
         })
     }
 
@@ -287,10 +294,10 @@ where
         m22: N,
         compute_u: bool,
         compute_v: bool,
-    ) -> (Option<UnitComplex<N>>, Vector2<N>, Option<UnitComplex<N>>)
+    ) -> (Option<GivensRotation<N>>, Vector2<N>, Option<GivensRotation<N>>)
     {
-        let two: N = ::convert(2.0f64);
-        let half: N = ::convert(0.5f64);
+        let two: N::Real = ::convert(2.0f64);
+        let half: N::Real = ::convert(0.5f64);
 
         let denom = (m11 + m22).hypot(m12) + (m11 - m22).hypot(m12);
 
@@ -298,24 +305,29 @@ where
         // 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 v1 = two * m11 * m22 / denom;
-        let v2 = half * denom;
+        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 = Vector2::new(m11 * m12, v1 * v1 - m11 * m11).normalize();
+            let (csv, sgn_v) = GivensRotation::new(m11 * m12, v1 * v1 - m11 * m11);
+            v1 *= sgn_v;
+            v2 *= sgn_v;
 
             if compute_v {
-                v_t = Some(UnitComplex::new_unchecked(Complex::new(csv.x, csv.y)));
+                v_t = Some(csv);
             }
 
             if compute_u {
-                let cu = (m11 * csv.x + m12 * csv.y) / v1;
-                let su = (m22 * csv.y) / v1;
+                let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1;
+                let su = (m22 * csv.s()) / v1;
+                let (csu, sgn_u) = GivensRotation::new(cu, su);
 
-                u = Some(UnitComplex::new_unchecked(Complex::new(cu, su)));
+                v1 *= sgn_u;
+                v2 *= sgn_u;
+                u = Some(csu);
             }
         }
 
@@ -338,11 +350,13 @@ where
     */
 
     fn delimit_subproblem(
-        b: &mut Bidiagonal<N, R, C>,
+        diagonal: &mut VectorN<N, DimMinimum<R, C>>,
+        off_diagonal: &mut VectorN<N, 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,
         end: usize,
-        eps: N,
+        eps: N::Real,
     ) -> (usize, usize)
     {
         let mut n = end;
@@ -350,20 +364,20 @@ where
         while n > 0 {
             let m = n - 1;
 
-            if b.off_diagonal[m].is_zero()
-                || b.off_diagonal[m].abs() <= eps * (b.diagonal[n].abs() + b.diagonal[m].abs())
+            if off_diagonal[m].is_zero()
+                || off_diagonal[m].modulus() <= eps * (diagonal[n].modulus() + diagonal[m].modulus())
             {
-                b.off_diagonal[m] = N::zero();
-            } else if b.diagonal[m].abs() <= eps {
-                b.diagonal[m] = N::zero();
-                Self::cancel_horizontal_off_diagonal_elt(b, u, v_t, m, m + 1);
+                off_diagonal[m] = N::zero();
+            } else if diagonal[m].modulus() <= eps {
+                diagonal[m] = N::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(b, u, v_t, m - 1);
+                    Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m - 1);
                 }
-            } else if b.diagonal[n].abs() <= eps {
-                b.diagonal[n] = N::zero();
-                Self::cancel_vertical_off_diagonal_elt(b, u, v_t, m);
+            } else if diagonal[n].modulus() <= eps {
+                diagonal[n] = N::zero();
+                Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m);
             } else {
                 break;
             }
@@ -379,18 +393,18 @@ where
         while new_start > 0 {
             let m = new_start - 1;
 
-            if b.off_diagonal[m].abs() <= eps * (b.diagonal[new_start].abs() + b.diagonal[m].abs())
+            if off_diagonal[m].modulus() <= eps * (diagonal[new_start].modulus() + diagonal[m].modulus())
             {
-                b.off_diagonal[m] = N::zero();
+                off_diagonal[m] = N::zero();
                 break;
             }
             // FIXME: write a test that enters this case.
-            else if b.diagonal[m].abs() <= eps {
-                b.diagonal[m] = N::zero();
-                Self::cancel_horizontal_off_diagonal_elt(b, u, v_t, m, n);
+            else if diagonal[m].modulus() <= eps {
+                diagonal[m] = N::zero();
+                Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, n);
 
                 if m != 0 {
-                    Self::cancel_vertical_off_diagonal_elt(b, u, v_t, m - 1);
+                    Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m - 1);
                 }
                 break;
             }
@@ -403,21 +417,23 @@ where
 
     // Cancels the i-th off-diagonal element using givens rotations.
     fn cancel_horizontal_off_diagonal_elt(
-        b: &mut Bidiagonal<N, R, C>,
+        diagonal: &mut VectorN<N, DimMinimum<R, C>>,
+        off_diagonal: &mut VectorN<N, 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,
         i: usize,
         end: usize,
     )
     {
-        let mut v = Vector2::new(b.off_diagonal[i], b.diagonal[i + 1]);
-        b.off_diagonal[i] = N::zero();
+        let mut v = Vector2::new(off_diagonal[i], diagonal[i + 1]);
+        off_diagonal[i] = N::zero();
 
         for k in i..end {
-            if let Some((rot, norm)) = givens::cancel_x(&v) {
-                b.diagonal[k + 1] = norm;
+            if let Some((rot, norm)) = GivensRotation::cancel_x(&v) {
+                diagonal[k + 1] = norm;
 
-                if b.is_upper_diagonal() {
+                if is_upper_diagonal {
                     if let Some(ref mut u) = *u {
                         rot.inverse()
                             .rotate_rows(&mut u.fixed_columns_with_step_mut::<U2>(i, k - i));
@@ -427,9 +443,9 @@ where
                 }
 
                 if k + 1 != end {
-                    v.x = -rot.sin_angle() * b.off_diagonal[k + 1];
-                    v.y = b.diagonal[k + 2];
-                    b.off_diagonal[k + 1] *= rot.cos_angle();
+                    v.x = -rot.s() * off_diagonal[k + 1];
+                    v.y = diagonal[k + 2];
+                    off_diagonal[k + 1] = off_diagonal[k + 1].scale(rot.c());
                 }
             } else {
                 break;
@@ -439,20 +455,22 @@ where
 
     // Cancels the i-th off-diagonal element using givens rotations.
     fn cancel_vertical_off_diagonal_elt(
-        b: &mut Bidiagonal<N, R, C>,
+        diagonal: &mut VectorN<N, DimMinimum<R, C>>,
+        off_diagonal: &mut VectorN<N, 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,
         i: usize,
     )
     {
-        let mut v = Vector2::new(b.diagonal[i], b.off_diagonal[i]);
-        b.off_diagonal[i] = N::zero();
+        let mut v = Vector2::new(diagonal[i], off_diagonal[i]);
+        off_diagonal[i] = N::zero();
 
         for k in (0..i + 1).rev() {
-            if let Some((rot, norm)) = givens::cancel_y(&v) {
-                b.diagonal[k] = norm;
+            if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
+                diagonal[k] = norm;
 
-                if b.is_upper_diagonal() {
+                if is_upper_diagonal {
                     if let Some(ref mut v_t) = *v_t {
                         rot.rotate(&mut v_t.fixed_rows_with_step_mut::<U2>(k, i - k));
                     }
@@ -462,9 +480,9 @@ where
                 }
 
                 if k > 0 {
-                    v.x = b.diagonal[k - 1];
-                    v.y = rot.sin_angle() * b.off_diagonal[k - 1];
-                    b.off_diagonal[k - 1] *= rot.cos_angle();
+                    v.x = diagonal[k - 1];
+                    v.y = rot.s() * off_diagonal[k - 1];
+                    off_diagonal[k - 1] = off_diagonal[k - 1].scale(rot.c());
                 }
             } else {
                 break;
@@ -474,12 +492,12 @@ where
 
     /// Computes the rank of the decomposed matrix, i.e., the number of singular values greater
     /// than `eps`.
-    pub fn rank(&self, eps: N) -> usize {
+    pub fn rank(&self, eps: N::Real) -> usize {
         assert!(
-            eps >= N::zero(),
+            eps >= N::Real::zero(),
             "SVD rank: the epsilon must be non-negative."
         );
-        self.singular_values.iter().filter(|e| **e > eps).count()
+        self.singular_values.iter().filter(|e| e.asum() > eps).count()
     }
 
     /// Rebuild the original matrix.
@@ -507,18 +525,18 @@ where
     /// Any singular value smaller than `eps` is assumed to be zero.
     /// Returns `Err` if the right- and left- singular vectors have not
     /// been computed at construction-time.
-    pub fn pseudo_inverse(mut self, eps: N) -> Result<MatrixMN<N, C, R>, &'static str>
+    pub fn pseudo_inverse(mut self, eps: N::Real) -> Result<MatrixMN<N, C, R>, &'static str>
     where
         DefaultAllocator: Allocator<N, C, R>,
     {
-        if eps < N::zero() {
+        if eps < N::Real::zero() {
             Err("SVD pseudo inverse: the epsilon must be non-negative.")
         }
         else {
             for i in 0..self.singular_values.len() {
                 let val = self.singular_values[i];
 
-                if val > eps {
+                if val.asum() > eps {
                     self.singular_values[i] = N::one() / val;
                 } else {
                     self.singular_values[i] = N::zero();
@@ -537,14 +555,14 @@ where
     pub fn solve<R2: Dim, C2: Dim, S2>(
         &self,
         b: &Matrix<N, R2, C2, S2>,
-        eps: N,
+        eps: N::Real,
     ) -> Result<MatrixMN<N, C, C2>, &'static str>
     where
         S2: Storage<N, R2, C2>,
         DefaultAllocator: Allocator<N, C, C2> + Allocator<N, DimMinimum<R, C>, C2>,
         ShapeConstraint: SameNumberOfRows<R, R2>,
     {
-        if eps < N::zero() {
+        if eps < N::Real::zero() {
             Err("SVD solve: the epsilon must be non-negative.")
         }
         else {
@@ -557,7 +575,7 @@ where
 
                         for i in 0..self.singular_values.len() {
                             let val = self.singular_values[i];
-                            if val > eps {
+                            if val.asum() > eps {
                                 col[i] /= val;
                             } else {
                                 col[i] = N::zero();
@@ -575,7 +593,7 @@ where
     }
 }
 
-impl<N: Real, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
+impl<N: Complex, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
 where
     DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
     DefaultAllocator: Allocator<N, R, C>
@@ -605,7 +623,7 @@ where
         self,
         compute_u: bool,
         compute_v: bool,
-        eps: N,
+        eps: N::Real,
         max_niter: usize,
     ) -> Option<SVD<N, R, C>>
     {
@@ -620,7 +638,7 @@ where
     /// Computes the rank of this matrix.
     ///
     /// All singular values below `eps` are considered equal to 0.
-    pub fn rank(&self, eps: N) -> usize {
+    pub fn rank(&self, eps: N::Real) -> usize {
         let svd = SVD::new(self.clone_owned(), false, false);
         svd.rank(eps)
     }
@@ -628,7 +646,7 @@ where
     /// Computes the pseudo-inverse of this matrix.
     ///
     /// All singular values below `eps` are considered equal to 0.
-    pub fn pseudo_inverse(self, eps: N) -> Result<MatrixMN<N, C, R>, &'static str>
+    pub fn pseudo_inverse(self, eps: N::Real) -> Result<MatrixMN<N, C, R>, &'static str>
     where
         DefaultAllocator: Allocator<N, C, R>,
     {
diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs
index 6084c1b4..b5437d83 100644
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};
 
-use num::Zero;
+use num::{Zero, One};
 use num_complex::Complex as NumComplex;
 use approx::AbsDiffEq;
 use std::ops::MulAssign;
@@ -119,6 +119,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
             q = Some(res.0);
             diag = res.1;
             off_diag = res.2;
+
+            println!("Tridiagonalization q: {:.5?}", q);
         } else {
             let res = SymmetricTridiagonal::new(m).unpack_tridiagonal();
             q = None;
@@ -150,6 +152,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                     let j = i + 1;
 
                     if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
+                        println!("Canceling: {:.5?} with norm: {:.5?}", rot, norm);
                         if i > start {
                             // Not the first iteration.
                             off_diag[i - 1] = norm;
@@ -160,19 +163,33 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                         let mij = off_diag[i];
 
                         let cc = rot.c() * rot.c();
-                        let ss = rot.s() * rot.s().conjugate();
-                        let cs = rot.c() * rot.s();
+                        let ss = rot.s().modulus_squared(); // rot.s() * rot.s().conjugate()
+                        let cs = rot.s().scale(rot.c());
 
-                        let b = cs * mij.conjugate() + cs.conjugate() * mij;
+                        // b = cs * mij.conjugate() + cs.conjugate() * mij
+                        let b = N::from_real((cs * mij.conjugate()).real() * ::convert(2.0));
 
-                        diag[i] = (cc * mii + ss * mjj) - b;
-                        diag[j] = (ss * mii + cc * mjj) + b;
-                        off_diag[i] = cs * (mii - mjj) + mij * cc - mij.conjugate() * rot.s() * rot.s();
+                        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]);
 
                         if i != n - 1 {
                             v.x = off_diag[i];
                             v.y = -rot.s() * off_diag[i + 1];
-                            off_diag[i + 1] *= rot.c();
+                            off_diag[i + 1] = off_diag[i + 1].scale(rot.c());
                         }
 
                         if let Some(ref mut q) = q {
@@ -197,12 +214,12 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                 diag[start + 0] = eigvals[0];
                 diag[start + 1] = eigvals[1];
 
-                println!("Eigvals: {:?}", eigvals);
-                println!("m: {}", m);
-                println!("Curr q: {:?}", q);
+                println!("Eigvals: {:.5?}", eigvals);
+                println!("m: {:.5}", m);
+                println!("Curr q: {:.5?}", q);
 
                 if let Some(ref mut q) = q {
-                    if let Some(rot) = GivensRotation::try_new(basis.x, basis.y, eps) {
+                    if let Some((rot, _)) = GivensRotation::try_new(basis.x, basis.y, eps) {
                         rot.rotate_rows(&mut q.fixed_columns_mut::<U2>(start));
                     }
                 }
diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs
index 9b54c596..cd4f7c3d 100644
--- a/src/linalg/symmetric_tridiagonal.rs
+++ b/src/linalg/symmetric_tridiagonal.rs
@@ -83,6 +83,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
                 m.ger_symm(-N::one(), &p, &axis.conjugate(), N::one());
                 m.ger_symm(-N::one(), &axis, &p.conjugate(), N::one());
                 m.ger_symm(dot * ::convert(2.0), &axis, &axis.conjugate(), N::one());
+                println!("The m: {}", m);
             }
         }
 
diff --git a/tests/linalg/bidiagonal.rs b/tests/linalg/bidiagonal.rs
index 28b1e3a9..dbb0c4fb 100644
--- a/tests/linalg/bidiagonal.rs
+++ b/tests/linalg/bidiagonal.rs
@@ -41,6 +41,7 @@ quickcheck! {
         relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
     }
 
+
     fn bidiagonal_static_square(m: Matrix4<RandComplex<f64>>) -> bool {
         let m = m.map(|e| e.0);
         let bidiagonal = m.bidiagonalize();
diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs
index 28a2b0d1..024384ce 100644
--- a/tests/linalg/eigen.rs
+++ b/tests/linalg/eigen.rs
@@ -9,10 +9,9 @@ mod quickcheck_tests {
     use std::cmp;
 
     quickcheck! {
-    /*
         fn symmetric_eigen(n: usize) -> bool {
             let n      = cmp::max(1, cmp::min(n, 10));
-            let m      = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);
+            let m      = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
             let eig    = m.clone().symmetric_eigen();
             let recomp = eig.recompose();
 
@@ -23,7 +22,7 @@ mod quickcheck_tests {
 
         fn symmetric_eigen_singular(n: usize) -> bool {
             let n      = cmp::max(1, cmp::min(n, 10));
-            let mut m  = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);
+            let mut m  = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
             m.row_mut(n / 2).fill(na::zero());
             m.column_mut(n / 2).fill(na::zero());
             let eig    = m.clone().symmetric_eigen();
@@ -35,7 +34,7 @@ mod quickcheck_tests {
         }
 
         fn symmetric_eigen_static_square_4x4(m: Matrix4<RandComplex<f64>>) -> bool {
-            let m      = m.map(|e| e.0);
+            let m      = m.map(|e| e.0).hermitian_part();
             let eig    = m.symmetric_eigen();
             let recomp = eig.recompose();
 
@@ -43,7 +42,6 @@ mod quickcheck_tests {
 
             relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
         }
-        */
 
         fn symmetric_eigen_static_square_3x3(m: Matrix3<RandComplex<f64>>) -> bool {
             let m      = m.map(|e| e.0).hermitian_part();
@@ -57,7 +55,6 @@ mod quickcheck_tests {
             relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
         }
 
-/*
         fn symmetric_eigen_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
             let m      = m.map(|e| e.0).hermitian_part();
             let eig    = m.symmetric_eigen();
@@ -68,11 +65,9 @@ mod quickcheck_tests {
 
             relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
         }
-        */
     }
 }
 
-/*
 // Test proposed on the issue #176 of rulinalg.
 #[test]
 fn symmetric_eigen_singular_24x24() {
@@ -115,7 +110,7 @@ fn symmetric_eigen_singular_24x24() {
         recomp.lower_triangle(),
         epsilon = 1.0e-5
     ));
-}*/
+}
 
 //  #[cfg(feature = "arbitrary")]
 //  quickcheck! {
diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs
index e84108ed..91f2002d 100644
--- a/tests/linalg/svd.rs
+++ b/tests/linalg/svd.rs
@@ -7,8 +7,11 @@ mod quickcheck_tests {
         DMatrix, DVector, Matrix2, Matrix2x5, Matrix3, Matrix3x5, Matrix4, Matrix5x2, Matrix5x3,
     };
     use std::cmp;
+    use core::helper::{RandScalar, RandComplex};
+
 
     quickcheck! {
+    /*
         fn svd(m: DMatrix<f64>) -> bool {
             if m.len() > 0 {
                 let svd = m.clone().svd(true, true);
@@ -68,6 +71,7 @@ mod quickcheck_tests {
             relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
         }
 
+
         fn svd_static_square(m: Matrix4<f64>) -> bool {
             let svd = m.svd(true, true);
             let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
@@ -78,18 +82,27 @@ mod quickcheck_tests {
             u.is_orthogonal(1.0e-5) &&
             v_t.is_orthogonal(1.0e-5)
         }
+                        */
 
-        fn svd_static_square_2x2(m: Matrix2<f64>) -> bool {
+
+        fn svd_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
+            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);
 
-            s.iter().all(|e| *e >= 0.0) &&
+            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) &&
             relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
             u.is_orthogonal(1.0e-5) &&
             v_t.is_orthogonal(1.0e-5)
         }
 
+/*
         fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
             if m.len() > 0 {
                 let svd = m.clone().svd(true, true);
@@ -140,9 +153,11 @@ mod quickcheck_tests {
 
             true
         }
+        */
     }
 }
 
+/*
 // Test proposed on the issue #176 of rulinalg.
 #[test]
 fn svd_singular() {
@@ -342,3 +357,5 @@ fn svd_err() {
     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
diff --git a/tests/linalg/tridiagonal.rs b/tests/linalg/tridiagonal.rs
index 5ecf32ff..fee0176a 100644
--- a/tests/linalg/tridiagonal.rs
+++ b/tests/linalg/tridiagonal.rs
@@ -15,6 +15,20 @@ quickcheck! {
         relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
     }
 
+    fn symm_tridiagonal_singular(n: usize) -> bool {
+        let n = cmp::max(1, cmp::min(n, 4));
+        let mut m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
+        m.row_mut(n / 2).fill(na::zero());
+        m.column_mut(n / 2).fill(na::zero());
+        let tri = m.clone().symmetric_tridiagonalize();
+        println!("Tri: {:?}", tri);
+        let recomp = tri.recompose();
+        println!("Recomp: {:?}", recomp);
+
+
+        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
+    }
+
     fn symm_tridiagonal_static_square(m: Matrix4<RandComplex<f64>>) -> bool {
         let m = m.map(|e| e.0).hermitian_part();
         let tri = m.symmetric_tridiagonalize();