WIP use Complex instead of Real whenever possible on the linalg module.

2019-03-02 19:33:49 +01:00 · 2019-03-02 19:33:49 +01:00 · 77f048b6b9
parent 7c91f2eeb5
commit 77f048b6b9
34 changed files with 878 additions and 409 deletions
--- a/Cargo.toml
+++ b/Cargo.toml
@ -26,7 +26,7 @@ arbitrary       = [ "quickcheck" ]
 serde-serialize = [ "serde", "serde_derive", "num-complex/serde" ]
 abomonation-serialize = [ "abomonation" ]
 sparse = [ ]
-debug = [ ]
+debug = [ "approx/num-complex", "rand/std" ]
 alloc = [ ]
 io = [ "pest", "pest_derive" ]

@ -53,3 +53,6 @@ rand_xorshift = "0.1"

 [workspace]
 members = [ "nalgebra-lapack", "nalgebra-glm" ]
+
+[patch.crates-io]
+alga = { path = "../alga/alga" }
--- a/src/base/blas.rs
+++ b/src/base/blas.rs
@ -13,6 +13,39 @@ use base::dimension::{Dim, Dynamic, U1, U2, U3, U4};
 use base::storage::{Storage, StorageMut};
 use base::{DefaultAllocator, Matrix, Scalar, SquareMatrix, Vector};

+
+// FIXME: find a way to avoid code duplication just for complex number support.
+impl<N: Complex, D: Dim, S: Storage<N, D>> Vector<N, D, S> {
+    /// Computes the index of the vector component with the largest complex or real absolute value.
+    ///
+    /// # Examples:
+    ///
+    /// ```
+    /// # use num_complex::Complex;
+    /// # use nalgebra::Vector3;
+    /// let vec = Vector3::new(Complex::new(11.0, 3.0), Complex::new(-15.0, 0.0), Complex::new(13.0, 5.0));
+    /// assert_eq!(vec.icamax(), 2);
+    /// ```
+    #[inline]
+    pub fn icamax(&self) -> usize {
+        assert!(!self.is_empty(), "The input vector must not be empty.");
+
+        let mut the_max = unsafe { self.vget_unchecked(0).asum() };
+        let mut the_i = 0;
+
+        for i in 1..self.nrows() {
+            let val = unsafe { self.vget_unchecked(i).asum() };
+
+            if val > the_max {
+                the_max = val;
+                the_i = i;
+            }
+        }
+
+        the_i
+    }
+}
+
 impl<N: Scalar + PartialOrd, D: Dim, S: Storage<N, D>> Vector<N, D, S> {
    /// Computes the index and value of the vector component with the largest value.
    ///
@ -157,6 +190,41 @@ impl<N: Scalar + PartialOrd, D: Dim, S: Storage<N, D>> Vector<N, D, S> {
    }
 }

+// FIXME: find a way to avoid code duplication just for complex number support.
+impl<N: Complex, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
+    /// Computes the index of the matrix component with the largest absolute value.
+    ///
+    /// # Examples:
+    ///
+    /// ```
+    /// # use nalgebra::Matrix2x3;
+    /// let mat = Matrix2x3::new(Complex::new(11.0, 1.0), Complex::new(-12.0, 2.0), Complex::new(13.0, 3.0),
+    ///                          Complex::new(21.0, 43.0), Complex::new(22.0, 5.0), Complex::new(-23.0, 0.0);
+    /// assert_eq!(mat.iamax_full(), (1, 0));
+    /// ```
+    #[inline]
+    pub fn icamax_full(&self) -> (usize, usize) {
+        assert!(!self.is_empty(), "The input matrix must not be empty.");
+
+        let mut the_max = unsafe { self.get_unchecked((0, 0)).asum() };
+        let mut the_ij = (0, 0);
+
+        for j in 0..self.ncols() {
+            for i in 0..self.nrows() {
+                let val = unsafe { self.get_unchecked((i, j)).asum() };
+
+                if val > the_max {
+                    the_max = val;
+                    the_ij = (i, j);
+                }
+            }
+        }
+
+        the_ij
+    }
+}
+
+
 impl<N: Scalar + PartialOrd + Signed, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
    /// Computes the index of the matrix component with the largest absolute value.
    ///
@ -705,6 +773,56 @@ where
    }
 }

+// FIXME: duplicate code
+impl<N, R1: Dim, C1: Dim, S: StorageMut<N, R1, C1>> Matrix<N, R1, C1, S>
+    where N: Complex + Zero + ClosedAdd + ClosedMul
+{
+    /// Computes `self = alpha * x * y.transpose() + beta * self`.
+    ///
+    /// If `beta` is zero, `self` is never read.
+    ///
+    /// # Examples:
+    ///
+    /// ```
+    /// # use nalgebra::{Matrix2x3, Vector2, Vector3};
+    /// let mut mat = Matrix2x3::repeat(4.0);
+    /// let vec1 = Vector2::new(1.0, 2.0);
+    /// let vec2 = Vector3::new(0.1, 0.2, 0.3);
+    /// let expected = vec1 * vec2.transpose() * 10.0 + mat * 5.0;
+    ///
+    /// mat.ger(10.0, &vec1, &vec2, 5.0);
+    /// assert_eq!(mat, expected);
+    /// ```
+    #[inline]
+    pub fn gerc<D2: Dim, D3: Dim, SB, SC>(
+        &mut self,
+        alpha: N,
+        x: &Vector<N, D2, SB>,
+        y: &Vector<N, D3, SC>,
+        beta: N,
+    ) where
+        N: One,
+        SB: Storage<N, D2>,
+        SC: Storage<N, D3>,
+        ShapeConstraint: DimEq<R1, D2> + DimEq<C1, D3>,
+    {
+        let (nrows1, ncols1) = self.shape();
+        let dim2 = x.nrows();
+        let dim3 = y.nrows();
+
+        assert!(
+            nrows1 == dim2 && ncols1 == dim3,
+            "ger: dimensions mismatch."
+        );
+
+        for j in 0..ncols1 {
+            // FIXME: avoid bound checks.
+            let val = unsafe { y.vget_unchecked(j).conjugate() };
+            self.column_mut(j).axpy(alpha * val, x, beta);
+        }
+    }
+}
+
 impl<N, R1: Dim, C1: Dim, S: StorageMut<N, R1, C1>> Matrix<N, R1, C1, S>
 where N: Scalar + Zero + ClosedAdd + ClosedMul
 {
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@ -944,6 +944,47 @@ impl<N: Complex, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
            res
        }
    }
+
+    /// The conjugate of `self`.
+    #[inline]
+    pub fn conjugate(&self) -> MatrixMN<N, R, C>
+        where DefaultAllocator: Allocator<N, R, C> {
+        self.map(|e| e.conjugate())
+    }
+
+    /// Divides each component of `self` by the given real.
+    #[inline]
+    pub fn unscale(&self, real: N::Real) -> MatrixMN<N, R, C>
+        where DefaultAllocator: Allocator<N, R, C> {
+        self.map(|e| e.unscale(real))
+    }
+
+    /// Multiplies each component of `self` by the given real.
+    #[inline]
+    pub fn scale(&self, real: N::Real) -> MatrixMN<N, R, C>
+        where DefaultAllocator: Allocator<N, R, C> {
+        self.map(|e| e.scale(real))
+    }
+}
+
+impl<N: Complex, R: Dim, C: Dim, S: StorageMut<N, R, C>> Matrix<N, R, C, S> {
+    /// The conjugate of `self` computed in-place.
+    #[inline]
+    pub fn conjugate_mut(&mut self) {
+        self.apply(|e| e.conjugate())
+    }
+
+    /// Divides each component of `self` by the given real.
+    #[inline]
+    pub fn unscale_mut(&mut self, real: N::Real) {
+        self.apply(|e| e.unscale(real))
+    }
+
+    /// Multiplies each component of `self` by the given real.
+    #[inline]
+    pub fn scale_mut(&mut self, real: N::Real) {
+        self.apply(|e| e.scale(real))
+    }
 }

 impl<N: Complex, D: Dim, S: StorageMut<N, D, D>> Matrix<N, D, D, S> {
@ -975,7 +1016,7 @@ 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.
    #[inline]
    pub fn diagonal(&self) -> VectorN<N, D>
-    where DefaultAllocator: Allocator<N, D> {
+        where DefaultAllocator: Allocator<N, D> {
        assert!(
            self.is_square(),
            "Unable to get the diagonal of a non-square matrix."
@ -996,7 +1037,7 @@ impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
    /// Computes a trace of a square matrix, i.e., the sum of its diagonal elements.
    #[inline]
    pub fn trace(&self) -> N
-    where N: Ring {
+        where N: Ring {
        assert!(
            self.is_square(),
            "Cannot compute the trace of non-square matrix."
@ -1013,6 +1054,34 @@ impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
    }
 }

+impl<N: Complex, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
+    /// The symmetric part of `self`, i.e., `0.5 * (self + self.transpose())`.
+    #[inline]
+    pub fn symmetric_part(&self) -> MatrixMN<N, D, D>
+        where DefaultAllocator: Allocator<N, D, D> {
+        assert!(self.is_square(), "Cannot compute the symmetric part of a non-square matrix.");
+        let mut tr = self.transpose();
+        tr += self;
+        tr *= ::convert::<_, N>(0.5);
+        tr
+    }
+
+    /// The hermitian part of `self`, i.e., `0.5 * (self + self.conjugate_transpose())`.
+    #[inline]
+    pub fn hermitian_part(&self) -> MatrixMN<N, D, D>
+        where DefaultAllocator: Allocator<N, D, D> {
+        assert!(self.is_square(), "Cannot compute the hermitian part of a non-square matrix.");
+        let nrows = self.data.shape().0;
+
+        unsafe {
+            let mut tr = self.conjugate_transpose();
+            tr += self;
+            tr *= ::convert::<_, N>(0.5);
+            tr
+        }
+    }
+}
+
 impl<N: Scalar + One + Zero, D: DimAdd<U1> + IsNotStaticOne, S: Storage<N, D, D>> Matrix<N, D, D, S> {

    /// Yields the homogeneous matrix for this matrix, i.e., appending an additional dimension and
--- a/src/base/ops.rs
+++ b/src/base/ops.rs
@ -5,7 +5,7 @@ use std::ops::{
    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
 };

-use alga::general::{ClosedAdd, ClosedDiv, ClosedMul, ClosedNeg, ClosedSub};
+use alga::general::{Complex, ClosedAdd, ClosedDiv, ClosedMul, ClosedNeg, ClosedSub};

 use base::allocator::{Allocator, SameShapeAllocator, SameShapeC, SameShapeR};
 use base::constraint::{
@ -760,6 +760,25 @@ where
    }
 }

+// XXX: avoid code duplication.
+impl<N: Complex, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
+    /// Returns the absolute value of the component with the largest absolute value.
+    #[inline]
+    pub fn camax(&self) -> N::Real {
+        let mut max = N::Real::zero();
+
+        for e in self.iter() {
+            let ae = e.asum();
+
+            if ae > max {
+                max = ae;
+            }
+        }
+
+        max
+    }
+}
+
 impl<N: Scalar + PartialOrd + Signed, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
    /// Returns the absolute value of the component with the largest absolute value.
    #[inline]
--- a/src/base/properties.rs
+++ b/src/base/properties.rs
@ -2,7 +2,7 @@
 use approx::RelativeEq;
 use num::{One, Zero};

-use alga::general::{ClosedAdd, ClosedMul, Real};
+use alga::general::{ClosedAdd, ClosedMul, Real, Complex};

 use base::allocator::Allocator;
 use base::dimension::{Dim, DimMin};
@ -82,20 +82,23 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {

        true
    }
+}

+impl<N: Complex, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
    /// Checks that `Mᵀ × M = Id`.
    ///
    /// In this definition `Id` is approximately equal to the identity matrix with a relative error
    /// equal to `eps`.
    #[inline]
    pub fn is_orthogonal(&self, eps: N::Epsilon) -> bool
-    where
-        N: Zero + One + ClosedAdd + ClosedMul + RelativeEq,
-        S: Storage<N, R, C>,
-        N::Epsilon: Copy,
-        DefaultAllocator: Allocator<N, C, C>,
+        where
+            N: Zero + One + ClosedAdd + ClosedMul + RelativeEq,
+            S: Storage<N, R, C>,
+            N::Epsilon: Copy,
+            DefaultAllocator: Allocator<N, R, C> + Allocator<N, C, C>,
    {
-        (self.tr_mul(self)).is_identity(eps)
+        // FIXME: add a conjugate-transpose-mul
+        (self.conjugate().tr_mul(self)).is_identity(eps)
    }
 }

--- a/src/geometry/reflection.rs
+++ b/src/geometry/reflection.rs
@ -1,4 +1,4 @@
-use alga::general::Real;
+use alga::general::Complex;
 use base::allocator::Allocator;
 use base::constraint::{AreMultipliable, DimEq, SameNumberOfRows, ShapeConstraint};
 use base::{DefaultAllocator, Matrix, Scalar, Unit, Vector};
@ -13,7 +13,7 @@ pub struct Reflection<N: Scalar, D: Dim, S: Storage<N, D>> {
    bias: N,
 }

-impl<N: Real, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
+impl<N: Complex, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
    /// Creates a new reflection wrt the plane orthogonal to the given axis and bias.
    ///
    /// The bias is the position of the plane on the axis. In particular, a bias equal to zero
@ -21,7 +21,7 @@ impl<N: Real, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
    pub fn new(axis: Unit<Vector<N, D, S>>, bias: N) -> Self {
        Self {
            axis: axis.into_inner(),
-            bias: bias,
+            bias,
        }
    }

@ -35,7 +35,7 @@ impl<N: Real, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
        D: DimName,
        DefaultAllocator: Allocator<N, D>,
    {
-        let bias = pt.coords.dot(axis.as_ref());
+        let bias = axis.cdot(&pt.coords);
        Self::new(axis, bias)
    }

@ -56,7 +56,7 @@ impl<N: Real, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
            // dot product, and then mutably. Somehow, this allows significantly
            // better optimizations of the dot product from the compiler.
            let m_two: N = ::convert(-2.0f64);
-            let factor = (rhs.column(i).dot(&self.axis) - self.bias) * m_two;
+            let factor = (self.axis.cdot(&rhs.column(i)) - self.bias) * m_two;
            rhs.column_mut(i).axpy(factor, &self.axis, N::one());
        }
    }
@ -70,8 +70,9 @@ impl<N: Real, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
        S2: StorageMut<N, R2, C2>,
        S3: StorageMut<N, R2>,
        ShapeConstraint: DimEq<C2, D> + AreMultipliable<R2, C2, D, U1>,
+        DefaultAllocator: Allocator<N, D>
    {
-        rhs.mul_to(&self.axis, work);
+        rhs.mul_to(&self.axis.conjugate(), work);

        if !self.bias.is_zero() {
            work.add_scalar_mut(-self.bias);
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::Real;
+use alga::general::Complex;
 use allocator::Allocator;
 use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Unit, VectorN};
 use constraint::{DimEq, ShapeConstraint};
@ -38,7 +38,7 @@ use linalg::householder;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct Bidiagonal<N: Real, R: DimMin<C>, C: Dim>
+pub struct Bidiagonal<N: Complex, R: DimMin<C>, C: Dim>
 where
    DimMinimum<R, C>: DimSub<U1>,
    DefaultAllocator: Allocator<N, R, C>
@ -55,7 +55,7 @@ where
    upper_diagonal: bool,
 }

-impl<N: Real, R: DimMin<C>, C: Dim> Copy for Bidiagonal<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> Copy for Bidiagonal<N, R, C>
 where
    DimMinimum<R, C>: DimSub<U1>,
    DefaultAllocator: Allocator<N, R, C>
@ -66,7 +66,7 @@ where
    VectorN<N, DimDiff<DimMinimum<R, C>, U1>>: Copy,
 {}

-impl<N: Real, R: DimMin<C>, C: Dim> Bidiagonal<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> Bidiagonal<N, R, C>
 where
    DimMinimum<R, C>: DimSub<U1>,
    DefaultAllocator: Allocator<N, R, C>
@ -273,7 +273,7 @@ where
    }
 }

-// impl<N: Real, D: DimMin<D, Output = D> + DimSub<Dynamic>> Bidiagonal<N, D, D>
+// impl<N: Complex, D: DimMin<D, Output = D> + DimSub<Dynamic>> Bidiagonal<N, D, D>
 //     where DefaultAllocator: Allocator<N, D, D> +
 //                             Allocator<N, D> {
 //     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
@ -346,7 +346,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>,
    DefaultAllocator: Allocator<N, R, C>
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@ -1,7 +1,8 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::Real;
+use num::Zero;
+use alga::general::Complex;

 use allocator::Allocator;
 use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix};
@ -26,19 +27,19 @@ use storage::{Storage, StorageMut};
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct Cholesky<N: Real, D: Dim>
+pub struct Cholesky<N: Complex, D: Dim>
 where DefaultAllocator: Allocator<N, D, D>
 {
    chol: MatrixN<N, D>,
 }

-impl<N: Real, D: Dim> Copy for Cholesky<N, D>
+impl<N: Complex, D: Dim> Copy for Cholesky<N, D>
 where
    DefaultAllocator: Allocator<N, D, D>,
    MatrixN<N, D>: Copy,
 {}

-impl<N: Real, D: DimSub<Dynamic>> Cholesky<N, D>
+impl<N: Complex, D: DimSub<Dynamic>> Cholesky<N, D>
 where DefaultAllocator: Allocator<N, D, D>
 {
    /// Attempts to compute the Cholesky decomposition of `matrix`.
@ -62,7 +63,7 @@ where DefaultAllocator: Allocator<N, D, D>
            }

            let diag = unsafe { *matrix.get_unchecked((j, j)) };
-            if diag > N::zero() {
+            if diag.real() > N::Real::zero() {
                let denom = diag.sqrt();
                unsafe {
                    *matrix.get_unchecked_mut((j, j)) = denom;
@ -144,7 +145,7 @@ where DefaultAllocator: Allocator<N, D, D>
    }
 }

-impl<N: Real, D: DimSub<Dynamic>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
+impl<N: Complex, D: DimSub<Dynamic>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
 where DefaultAllocator: Allocator<N, D, D>
 {
    /// Attempts to compute the Cholesky decomposition of this matrix.
--- a/src/linalg/determinant.rs
+++ b/src/linalg/determinant.rs
@ -1,4 +1,4 @@
-use alga::general::Real;
+use alga::general::Complex;

 use base::allocator::Allocator;
 use base::dimension::DimMin;
@ -7,7 +7,7 @@ use base::{DefaultAllocator, SquareMatrix};

 use linalg::LU;

-impl<N: Real, D: DimMin<D, Output = D>, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
+impl<N: Complex, D: DimMin<D, Output = D>, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
    /// Computes the matrix determinant.
    ///
    /// If the matrix has a dimension larger than 3, an LU decomposition is used.
--- a/src/linalg/eigen.rs
+++ b/src/linalg/eigen.rs
@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Serialize, Deserialize};

-use alga::general::Real;
+use alga::general::Complex;
 use num_complex::Complex;
 use std::cmp;
 use std::fmt::Display;
@ -17,7 +17,7 @@ use geometry::{Reflection, UnitComplex};
 use linalg::householder;
 use linalg::RealSchur;

-/// Eigendecomposition of a matrix with real eigenvalues.
+/// Eigendecomposition of a real matrix with real eigenvalues (or complex eigen values for complex matrices).
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
 #[cfg_attr(
    feature = "serde-serialize",
@ -40,7 +40,7 @@ use linalg::RealSchur;
    )
 )]
 #[derive(Clone, Debug)]
-pub struct RealEigen<N: Real, D: Dim>
+pub struct Eigen<N: Complex, D: Dim>
 where
    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
@ -48,7 +48,7 @@ where
    pub eigenvalues: VectorN<N, D>,
 }

-impl<N: Real, D: Dim> Copy for RealEigen<N, D>
+impl<N: Complex, D: Dim> Copy for Eigen<N, D>
 where
    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
    MatrixN<N, D>: Copy,
@ -56,7 +56,7 @@ where
 {
 }

-impl<N: Real, D: Dim> RealEigen<N, D>
+impl<N: Complex, D: Dim> Eigen<N, D>
 where
    D: DimSub<U1>,                                   // For Hessenberg.
    ShapeConstraint: DimEq<Dynamic, DimDiff<D, U1>>, // For Hessenberg.
@ -68,8 +68,8 @@ where
    DefaultAllocator: Allocator<usize, D, D>,
    MatrixN<N, D>: Display,
 {
-    /// Computes the eigendecomposition of a diagonalizable matrix with real eigenvalues.
-    pub fn new(m: MatrixN<N, D>) -> Option<RealEigen<N, D>> {
+    /// Computes the eigendecomposition of a diagonalizable matrix with Complex eigenvalues.
+    pub fn new(m: MatrixN<N, D>) -> Option<Eigen<N, D>> {
        assert!(
            m.is_square(),
            "Unable to compute the eigendecomposition of a non-square matrix."
@ -80,7 +80,7 @@ where

        println!("Schur eigenvalues: {}", eigenvalues);

-        // Check that the eigenvalues are all real.
+        // Check that the eigenvalues are all Complex.
        for i in 0..dim - 1 {
            if !eigenvalues[(i + 1, i)].is_zero() {
                return None;
@ -112,8 +112,8 @@ where
            let _ = eigenvectors.column_mut(i).normalize_mut();
        }

-        Some(RealEigen {
-            eigenvectors: eigenvectors,
+        Some(Eigen {
+            eigenvectors,
            eigenvalues: eigenvalues.diagonal(),
        })
    }
--- a/src/linalg/full_piv_lu.rs
+++ b/src/linalg/full_piv_lu.rs
@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::Real;
+use alga::general::Complex;
 use allocator::Allocator;
 use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN};
 use constraint::{SameNumberOfRows, ShapeConstraint};
@ -32,7 +32,7 @@ use linalg::PermutationSequence;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct FullPivLU<N: Real, R: DimMin<C>, C: Dim>
+pub struct FullPivLU<N: Complex, R: DimMin<C>, C: Dim>
 where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>
 {
    lu: MatrixMN<N, R, C>,
@ -40,14 +40,14 @@ where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimu
    q: PermutationSequence<DimMinimum<R, C>>,
 }

-impl<N: Real, R: DimMin<C>, C: Dim> Copy for FullPivLU<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> Copy for FullPivLU<N, R, C>
 where
    DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
    MatrixMN<N, R, C>: Copy,
    PermutationSequence<DimMinimum<R, C>>: Copy,
 {}

-impl<N: Real, R: DimMin<C>, C: Dim> FullPivLU<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> FullPivLU<N, R, C>
 where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>
 {
    /// Computes the LU decomposition with full pivoting of `matrix`.
@ -69,7 +69,7 @@ where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimu
        }

        for i in 0..min_nrows_ncols.value() {
-            let piv = matrix.slice_range(i.., i..).iamax_full();
+            let piv = matrix.slice_range(i.., i..).icamax_full();
            let row_piv = piv.0 + i;
            let col_piv = piv.1 + i;
            let diag = matrix[(row_piv, col_piv)];
@ -156,7 +156,7 @@ where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimu
    }
 }

-impl<N: Real, D: DimMin<D, Output = D>> FullPivLU<N, D, D>
+impl<N: Complex, D: DimMin<D, Output = D>> FullPivLU<N, D, D>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), D>
 {
    /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
@ -261,7 +261,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), D>
    }
 }

-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 DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>
 {
    /// Computes the LU decomposition with full pivoting of `matrix`.
--- a/src/linalg/givens.rs
+++ b/src/linalg/givens.rs
@ -1,36 +1,175 @@
 //! Construction of givens rotations.

-use alga::general::Real;
-use num_complex::Complex;
+use alga::general::{Complex, Real};
+use num_complex::Complex as NumComplex;

-use base::dimension::U2;
-use base::storage::Storage;
-use base::Vector;
+use base::dimension::{Dim, U2};
+use base::constraint::{ShapeConstraint, DimEq};
+use base::storage::{Storage, StorageMut};
+use base::{Vector, Matrix};

 use geometry::UnitComplex;

+/// A Givens rotation.
+pub struct GivensRotation<N: Complex> {
+    c: N,
+    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 = Complex::new(v[0], -v[1]);
+        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 = Complex::new(v[1], v[0]);
+        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> {
+    /// Initializes a Givens rotation form its non-normalized cosine an sine components.
+    pub fn new(c: N, s: N) -> Self {
+        let denom = (c.modulus_squared() + s.modulus_squared()).sqrt();
+        Self {
+            c: c.unscale(denom),
+            s: s.unscale(denom)
+        }
+    }
+
+    /// 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> {
+        let denom = (c.modulus_squared() + s.modulus_squared()).sqrt();
+
+        if denom > eps {
+            Some(Self {
+                c: c.unscale(denom),
+                s: s.unscale(denom)
+            })
+        } else {
+            None
+        }
+    }
+
+    /// 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<S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(Self, 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 r = sign0.scale(denom);
+            Some((Self { c, s }, r))
+        } else {
+            None
+        }
+    }
+
+    /// 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<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);
+            Some((Self { c, s }, r))
+        } else {
+            None
+        }
+    }
+
+    /// The cos part of this roration.
+    pub fn c(&self) -> N {
+        self.c
+    }
+
+    /// The sin part of this roration.
+    pub fn s(&self) -> N {
+        self.s
+    }
+
+    /// The inverse of this givens rotation.
+    pub fn inverse(&self) -> Self {
+        Self { c: self.c, s: -self.s.conjugate() }
+    }
+
+    /// Performs the multiplication `rhs = self * rhs` in-place.
+    pub fn rotate<R2: Dim, C2: Dim, S2: StorageMut<N, R2, C2>>(
+        &self,
+        rhs: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        ShapeConstraint: DimEq<R2, U2>,
+    {
+        assert_eq!(
+            rhs.nrows(),
+            2,
+            "Unit complex rotation: the input matrix must have exactly two rows."
+        );
+        let s = self.s;
+        let c = self.c;
+
+        for j in 0..rhs.ncols() {
+            unsafe {
+                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.conjugate() * b;
+            }
+        }
+    }
+
+    /// Performs the multiplication `lhs = lhs * self` in-place.
+    pub fn rotate_rows<R2: Dim, C2: Dim, S2: StorageMut<N, R2, C2>>(
+        &self,
+        lhs: &mut Matrix<N, R2, C2, S2>,
+    ) where
+        ShapeConstraint: DimEq<C2, U2>,
+    {
+        assert_eq!(
+            lhs.ncols(),
+            2,
+            "Unit complex rotation: the input matrix must have exactly two columns."
+        );
+        let s = self.s;
+        let c = self.c;
+
+        // FIXME: can we optimize that to iterate on one column at a time ?
+        for j in 0..lhs.nrows() {
+            unsafe {
+                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.conjugate() * b;
+            }
+        }
+    }
+}
+
--- a/src/linalg/hessenberg.rs
+++ b/src/linalg/hessenberg.rs
@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::Real;
+use alga::general::Complex;
 use allocator::Allocator;
 use base::{DefaultAllocator, MatrixMN, MatrixN, SquareMatrix, VectorN};
 use constraint::{DimEq, ShapeConstraint};
@ -31,21 +31,21 @@ use linalg::householder;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct Hessenberg<N: Real, D: DimSub<U1>>
+pub struct Hessenberg<N: Complex, D: DimSub<U1>>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
 {
    hess: MatrixN<N, D>,
    subdiag: VectorN<N, DimDiff<D, U1>>,
 }

-impl<N: Real, D: DimSub<U1>> Copy for Hessenberg<N, D>
+impl<N: Complex, D: DimSub<U1>> Copy for Hessenberg<N, D>
 where
    DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
    MatrixN<N, D>: Copy,
    VectorN<N, DimDiff<D, U1>>: Copy,
 {}

-impl<N: Real, D: DimSub<U1>> Hessenberg<N, D>
+impl<N: Complex, D: DimSub<U1>> Hessenberg<N, D>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimDiff<D, U1>>
 {
    /// Computes the Hessenberg decomposition using householder reflections.
@ -137,7 +137,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimD
    }
 }

-impl<N: Real, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
+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>>
 {
    /// Computes the Hessenberg decomposition of this matrix using householder reflections.
--- a/src/linalg/householder.rs
+++ b/src/linalg/householder.rs
@ -1,6 +1,7 @@
 //! Construction of householder elementary reflections.

-use alga::general::Real;
+use num::Zero;
+use alga::general::Complex;
 use allocator::Allocator;
 use base::{DefaultAllocator, MatrixMN, MatrixN, Unit, Vector, VectorN};
 use dimension::Dim;
@ -15,35 +16,34 @@ use geometry::Reflection;
 /// `column` after reflection and `false` if no reflection was necessary.
 #[doc(hidden)]
 #[inline(always)]
-pub fn reflection_axis_mut<N: Real, D: Dim, S: StorageMut<N, D>>(
+pub fn reflection_axis_mut<N: Complex, D: Dim, S: StorageMut<N, D>>(
    column: &mut Vector<N, D, S>,
 ) -> (N, bool) {
    let reflection_sq_norm = column.norm_squared();
-    let mut reflection_norm = reflection_sq_norm.sqrt();
+    let reflection_norm = reflection_sq_norm.sqrt();

    let factor;
-    unsafe {
-        if *column.vget_unchecked(0) > N::zero() {
-            reflection_norm = -reflection_norm;
-        }
+    let scaled_norm;

-        factor =
-            (reflection_sq_norm - *column.vget_unchecked(0) * reflection_norm) * ::convert(2.0);
-        *column.vget_unchecked_mut(0) -= reflection_norm;
-    }
+    unsafe {
+        let (modulus, exp) = column.vget_unchecked(0).to_exp();
+        scaled_norm = exp.scale(reflection_norm);
+        factor = (reflection_sq_norm + modulus * reflection_norm) * ::convert(2.0);
+        *column.vget_unchecked_mut(0) += scaled_norm;
+    };

    if !factor.is_zero() {
-        *column /= factor.sqrt();
-        (reflection_norm, true)
+        column.unscale_mut(factor.sqrt());
+        (-scaled_norm, true)
    } else {
-        (reflection_norm, false)
+        (-scaled_norm, false)
    }
 }

 /// Uses an householder reflection to zero out the `icol`-th column, starting with the `shift + 1`-th
 /// subdiagonal element.
 #[doc(hidden)]
-pub fn clear_column_unchecked<N: Real, R: Dim, C: Dim>(
+pub fn clear_column_unchecked<N: Complex, R: Dim, C: Dim>(
    matrix: &mut MatrixMN<N, R, C>,
    diag_elt: &mut N,
    icol: usize,
@ -70,7 +70,7 @@ pub fn clear_column_unchecked<N: Real, R: Dim, C: Dim>(
 /// Uses an hoseholder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
 /// superdiagonal element.
 #[doc(hidden)]
-pub fn clear_row_unchecked<N: Real, R: Dim, C: Dim>(
+pub fn clear_row_unchecked<N: Complex, R: Dim, C: Dim>(
    matrix: &mut MatrixMN<N, R, C>,
    diag_elt: &mut N,
    axis_packed: &mut VectorN<N, C>,
@ -94,7 +94,7 @@ pub fn clear_row_unchecked<N: Real, R: Dim, C: Dim>(
            &mut work.rows_range_mut(irow + 1..),
        );
        top.columns_range_mut(irow + shift..)
-            .tr_copy_from(refl.axis());
+            .tr_copy_from(&refl.axis());
    } else {
        top.columns_range_mut(irow + shift..).tr_copy_from(&axis);
    }
@ -104,7 +104,7 @@ pub fn clear_row_unchecked<N: Real, R: Dim, C: Dim>(
 /// the lower-diagonal element of the given matrix.
 /// matrices.
 #[doc(hidden)]
-pub fn assemble_q<N: Real, D: Dim>(m: &MatrixN<N, D>) -> MatrixN<N, D>
+pub fn assemble_q<N: Complex, D: Dim>(m: &MatrixN<N, D>) -> MatrixN<N, D>
 where DefaultAllocator: Allocator<N, D, D> {
    assert!(m.is_square());
    let dim = m.data.shape().0;
--- a/src/linalg/inverse.rs
+++ b/src/linalg/inverse.rs
@ -1,4 +1,4 @@
-use alga::general::Real;
+use alga::general::Complex;

 use base::allocator::Allocator;
 use base::dimension::Dim;
@ -7,7 +7,7 @@ use base::{DefaultAllocator, MatrixN, SquareMatrix};

 use linalg::lu;

-impl<N: Real, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
+impl<N: Complex, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
    /// Attempts to invert this matrix.
    #[inline]
    pub fn try_inverse(self) -> Option<MatrixN<N, D>>
@ -21,7 +21,7 @@ impl<N: Real, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
    }
 }

-impl<N: Real, D: Dim, S: StorageMut<N, D, D>> SquareMatrix<N, D, S> {
+impl<N: Complex, D: Dim, S: StorageMut<N, D, D>> SquareMatrix<N, D, S> {
    /// Attempts to invert this matrix in-place. Returns `false` and leaves `self` untouched if
    /// inversion fails.
    #[inline]
@ -115,7 +115,7 @@ impl<N: Real, D: Dim, S: StorageMut<N, D, D>> SquareMatrix<N, D, S> {
 }

 // NOTE: this is an extremely efficient, loop-unrolled matrix inverse from MESA (MIT licensed).
-fn do_inverse4<N: Real, D: Dim, S: StorageMut<N, D, D>>(
+fn do_inverse4<N: Complex, D: Dim, S: StorageMut<N, D, D>>(
    m: &MatrixN<N, D>,
    out: &mut SquareMatrix<N, D, S>,
 ) -> bool
--- a/src/linalg/lu.rs
+++ b/src/linalg/lu.rs
@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::{Field, Real};
+use alga::general::{Field, Complex};
 use allocator::{Allocator, Reallocator};
 use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar};
 use constraint::{SameNumberOfRows, ShapeConstraint};
@ -32,14 +32,14 @@ use linalg::PermutationSequence;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct LU<N: Real, R: DimMin<C>, C: Dim>
+pub struct LU<N: Complex, R: DimMin<C>, C: Dim>
 where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>
 {
    lu: MatrixMN<N, R, C>,
    p: PermutationSequence<DimMinimum<R, C>>,
 }

-impl<N: Real, R: DimMin<C>, C: Dim> Copy for LU<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> Copy for LU<N, R, C>
 where
    DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
    MatrixMN<N, R, C>: Copy,
@ -49,7 +49,7 @@ where
 /// Performs a LU decomposition to overwrite `out` with the inverse of `matrix`.
 ///
 /// If `matrix` is not invertible, `false` is returned and `out` may contain invalid data.
-pub fn try_invert_to<N: Real, D: Dim, S>(
+pub fn try_invert_to<N: Complex, D: Dim, S>(
    mut matrix: MatrixN<N, D>,
    out: &mut Matrix<N, D, D, S>,
 ) -> bool
@ -66,7 +66,7 @@ where
    out.fill_with_identity();

    for i in 0..dim {
-        let piv = matrix.slice_range(i.., i).iamax() + i;
+        let piv = matrix.slice_range(i.., i).icamax() + i;
        let diag = matrix[(piv, i)];

        if diag.is_zero() {
@ -86,7 +86,7 @@ where
    matrix.solve_upper_triangular_mut(out)
 }

-impl<N: Real, R: DimMin<C>, C: Dim> LU<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> LU<N, R, C>
 where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>
 {
    /// Computes the LU decomposition with partial (row) pivoting of `matrix`.
@ -101,7 +101,7 @@ where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimu
        }

        for i in 0..min_nrows_ncols.value() {
-            let piv = matrix.slice_range(i.., i).iamax() + i;
+            let piv = matrix.slice_range(i.., i).icamax() + i;
            let diag = matrix[(piv, i)];

            if diag.is_zero() {
@ -197,7 +197,7 @@ where DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimu
    }
 }

-impl<N: Real, D: DimMin<D, Output = D>> LU<N, D, D>
+impl<N: Complex, D: DimMin<D, Output = D>> LU<N, D, D>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), D>
 {
    /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
@ -368,7 +368,7 @@ pub fn gauss_step_swap<N, R: Dim, C: Dim, S>(
    }
 }

-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 DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>
 {
    /// Computes the LU decomposition with partial (row) pivoting of `matrix`.
--- a/src/linalg/qr.rs
+++ b/src/linalg/qr.rs
@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::Real;
+use alga::general::Complex;
 use allocator::{Allocator, Reallocator};
 use base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Unit, VectorN};
 use constraint::{SameNumberOfRows, ShapeConstraint};
@ -32,21 +32,21 @@ use linalg::householder;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct QR<N: Real, R: DimMin<C>, C: Dim>
+pub struct QR<N: Complex, R: DimMin<C>, C: Dim>
 where DefaultAllocator: Allocator<N, R, C> + Allocator<N, DimMinimum<R, C>>
 {
    qr: MatrixMN<N, R, C>,
    diag: VectorN<N, DimMinimum<R, C>>,
 }

-impl<N: Real, R: DimMin<C>, C: Dim> Copy for QR<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> Copy for QR<N, R, C>
 where
    DefaultAllocator: Allocator<N, R, C> + Allocator<N, DimMinimum<R, C>>,
    MatrixMN<N, R, C>: Copy,
    VectorN<N, DimMinimum<R, C>>: Copy,
 {}

-impl<N: Real, R: DimMin<C>, C: Dim> QR<N, R, C>
+impl<N: Complex, R: DimMin<C>, C: Dim> QR<N, R, C>
 where DefaultAllocator: Allocator<N, R, C> + Allocator<N, R> + Allocator<N, DimMinimum<R, C>>
 {
    /// Computes the QR decomposition using householder reflections.
@ -162,7 +162,7 @@ where DefaultAllocator: Allocator<N, R, C> + Allocator<N, R> + Allocator<N, DimM
    }
 }

-impl<N: Real, D: DimMin<D, Output = D>> QR<N, D, D>
+impl<N: Complex, D: DimMin<D, Output = D>> QR<N, D, D>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
 {
    /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
@ -294,7 +294,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
    // }
 }

-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 DefaultAllocator: Allocator<N, R, C> + Allocator<N, R> + Allocator<N, DimMinimum<R, C>>
 {
    /// Computes the QR decomposition of this matrix.
--- a/src/linalg/schur.rs
+++ b/src/linalg/schur.rs
@ -1,8 +1,9 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::Real;
-use num_complex::Complex;
+use approx::AbsDiffEq;
+use alga::general::{Complex, Real};
+use num_complex::Complex as NumComplex;
 use std::cmp;

 use allocator::Allocator;
@ -14,6 +15,7 @@ use constraint::{DimEq, ShapeConstraint};
 use geometry::{Reflection, UnitComplex};
 use linalg::householder;
 use linalg::Hessenberg;
+use linalg::givens::GivensRotation;

 /// Real Schur decomposition of a square matrix.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
@ -32,20 +34,20 @@ use linalg::Hessenberg;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct RealSchur<N: Real, D: Dim>
+pub struct RealSchur<N: Complex, D: Dim>
 where DefaultAllocator: Allocator<N, D, D>
 {
    q: MatrixN<N, D>,
    t: MatrixN<N, D>,
 }

-impl<N: Real, D: Dim> Copy for RealSchur<N, D>
+impl<N: Complex, D: Dim> Copy for RealSchur<N, D>
 where
    DefaultAllocator: Allocator<N, D, D>,
    MatrixN<N, D>: Copy,
 {}

-impl<N: Real, D: Dim> RealSchur<N, D>
+impl<N: Complex, D: Dim> RealSchur<N, D>
 where
    D: DimSub<U1>,                                   // For Hessenberg.
    ShapeConstraint: DimEq<Dynamic, DimDiff<D, U1>>, // For Hessenberg.
@ -56,7 +58,7 @@ where
 {
    /// Computes the Schur decomposition of a square matrix.
    pub fn new(m: MatrixN<N, D>) -> Self {
-        Self::try_new(m, N::default_epsilon(), 0).unwrap()
+        Self::try_new(m, N::Real::default_epsilon(), 0).unwrap()
    }

    /// Attempts to compute the Schur decomposition of a square matrix.
@ -70,7 +72,7 @@ where
    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
    /// continues indefinitely until convergence.
-    pub fn try_new(m: MatrixN<N, D>, eps: N, max_niter: usize) -> Option<Self> {
+    pub fn try_new(m: MatrixN<N, D>, eps: N::Real, max_niter: usize) -> Option<Self> {
        let mut work = unsafe { VectorN::new_uninitialized_generic(m.data.shape().0, U1) };

        Self::do_decompose(m, &mut work, eps, max_niter, true).map(|(q, t)| RealSchur {
@ -82,7 +84,7 @@ where
    fn do_decompose(
        mut m: MatrixN<N, D>,
        work: &mut VectorN<N, D>,
-        eps: N,
+        eps: N::Real,
        max_niter: usize,
        compute_q: bool,
    ) -> Option<(Option<MatrixN<N, D>>, MatrixN<N, D>)>
@ -111,8 +113,8 @@ where
            return decompose_2x2(m, compute_q);
        }

-        let amax_m = m.amax();
-        m /= amax_m;
+        let amax_m = m.camax();
+        m.unscale_mut(amax_m);

        let hess = Hessenberg::new_with_workspace(m, work);
        let mut q;
@ -259,7 +261,7 @@ where
            }
        }

-        t *= amax_m;
+        t.scale_mut(amax_m);

        Some((q, t))
    }
@ -289,8 +291,9 @@ where
    }

    /// Computes the complex eigenvalues of the decomposed matrix.
-    fn do_complex_eigenvalues(t: &MatrixN<N, D>, out: &mut VectorN<Complex<N>, D>)
-    where DefaultAllocator: Allocator<Complex<N>, D> {
+    fn do_complex_eigenvalues(t: &MatrixN<N, D>, out: &mut VectorN<NumComplex<N>, D>)
+    where N: Real,
+          DefaultAllocator: Allocator<NumComplex<N>, D> {
        let dim = t.nrows();
        let mut m = 0;

@ -298,7 +301,7 @@ where
            let n = m + 1;

            if t[(n, m)].is_zero() {
-                out[m] = Complex::new(t[(m, m)], N::zero());
+                out[m] = NumComplex::new(t[(m, m)], N::zero());
                m += 1;
            } else {
                // Solve the 2x2 eigenvalue subproblem.
@ -313,21 +316,21 @@ where

                // All 2x2 blocks have negative discriminant because we already decoupled those
                // with positive eigenvalues..
-                let sqrt_discr = Complex::new(N::zero(), (-discr).sqrt());
+                let sqrt_discr = NumComplex::new(N::zero(), (-discr).sqrt());

-                out[m] = Complex::new(tra * ::convert(0.5), N::zero()) + sqrt_discr;
-                out[m + 1] = Complex::new(tra * ::convert(0.5), N::zero()) - sqrt_discr;
+                out[m] = NumComplex::new(tra * ::convert(0.5), N::zero()) + sqrt_discr;
+                out[m + 1] = NumComplex::new(tra * ::convert(0.5), N::zero()) - sqrt_discr;

                m += 2;
            }
        }

        if m == dim - 1 {
-            out[m] = Complex::new(t[(m, m)], N::zero());
+            out[m] = NumComplex::new(t[(m, m)], N::zero());
        }
    }

-    fn delimit_subproblem(t: &mut MatrixN<N, D>, eps: N, end: usize) -> (usize, usize)
+    fn delimit_subproblem(t: &mut MatrixN<N, D>, eps: N::Real, end: usize) -> (usize, usize)
    where
        D: DimSub<U1>,
        DefaultAllocator: Allocator<N, DimDiff<D, U1>>,
@ -337,7 +340,7 @@ where
        while n > 0 {
            let m = n - 1;

-            if t[(n, m)].abs() <= eps * (t[(n, n)].abs() + t[(m, m)].abs()) {
+            if t[(n, m)].modulus() <= eps * (t[(n, n)].modulus() + t[(m, m)].modulus()) {
                t[(n, m)] = N::zero();
            } else {
                break;
@ -356,7 +359,7 @@ where

            let off_diag = t[(new_start, m)];
            if off_diag.is_zero()
-                || off_diag.abs() <= eps * (t[(new_start, new_start)].abs() + t[(m, m)].abs())
+                || off_diag.modulus() <= eps * (t[(new_start, new_start)].modulus() + t[(m, m)].modulus())
            {
                t[(new_start, m)] = N::zero();
                break;
@ -387,15 +390,16 @@ where
    }

    /// Computes the complex eigenvalues of the decomposed matrix.
-    pub fn complex_eigenvalues(&self) -> VectorN<Complex<N>, D>
-    where DefaultAllocator: Allocator<Complex<N>, D> {
+    pub fn complex_eigenvalues(&self) -> VectorN<NumComplex<N>, D>
+    where N: Real,
+          DefaultAllocator: Allocator<NumComplex<N>, D> {
        let mut out = unsafe { VectorN::new_uninitialized_generic(self.t.data.shape().0, U1) };
        Self::do_complex_eigenvalues(&self.t, &mut out);
        out
    }
 }

-fn decompose_2x2<N: Real, D: Dim>(
+fn decompose_2x2<N: Complex, D: Dim>(
    mut m: MatrixN<N, D>,
    compute_q: bool,
 ) -> Option<(Option<MatrixN<N, D>>, MatrixN<N, D>)>
@ -412,13 +416,12 @@ where
            rot.rotate_rows(&mut m);

            if compute_q {
-                let c = rot.into_inner();
                // XXX: we have to build the matrix manually because
                // rot.to_rotation_matrix().unwrap() causes an ICE.
                q = Some(MatrixN::from_column_slice_generic(
                    dim,
                    dim,
-                    &[c.re, c.im, -c.im, c.re],
+                    &[rot.c(), rot.s(), -rot.s().conjugate(), rot.c().conjugate()],
                ));
            }
        }
@ -432,7 +435,7 @@ where
    Some((q, m))
 }

-fn compute_2x2_eigvals<N: Real, S: Storage<N, U2, U2>>(
+fn compute_2x2_eigvals<N: Complex, S: Storage<N, U2, U2>>(
    m: &SquareMatrix<N, U2, S>,
 ) -> Option<(N, N)> {
    // Solve the 2x2 eigenvalue subproblem.
@ -447,13 +450,10 @@ fn compute_2x2_eigvals<N: Real, S: Storage<N, U2, U2>>(
    let val = (h00 - h11) * ::convert(0.5);
    let discr = h10 * h01 + val * val;

-    if discr >= N::zero() {
-        let sqrt_discr = discr.sqrt();
+    discr.try_sqrt().map(|sqrt_discr| {
        let half_tra = (h00 + h11) * ::convert(0.5);
-        Some((half_tra + sqrt_discr, half_tra - sqrt_discr))
-    } else {
-        None
-    }
+        (half_tra + sqrt_discr, half_tra - sqrt_discr)
+    })
 }

 // Computes the 2x2 transformation that upper-triangulates a 2x2 matrix with real eigenvalues.
@ -461,9 +461,9 @@ fn compute_2x2_eigvals<N: Real, S: Storage<N, U2, U2>>(
 ///
 /// Returns `None` if the matrix has complex eigenvalues, or is upper-triangular. In both case,
 /// the basis is the identity.
-fn compute_2x2_basis<N: Real, S: Storage<N, U2, U2>>(
+fn compute_2x2_basis<N: Complex, S: Storage<N, U2, U2>>(
    m: &SquareMatrix<N, U2, S>,
-) -> Option<UnitComplex<N>> {
+) -> Option<GivensRotation<N>> {
    let h10 = m[(1, 0)];

    if h10.is_zero() {
@ -477,19 +477,17 @@ fn compute_2x2_basis<N: Real, S: Storage<N, U2, U2>>(
        // NOTE: Choose the one that yields a larger x component.
        // This is necessary for numerical stability of the normalization of the complex
        // number.
-        let basis = if x1.abs() > x2.abs() {
-            Complex::new(x1, -h10)
+        if x1.modulus() > x2.modulus() {
+            Some(GivensRotation::new(x1, -h10))
        } else {
-            Complex::new(x2, -h10)
-        };
-
-        Some(UnitComplex::from_complex(basis))
+            Some(GivensRotation::new(x2, -h10))
+        }
    } else {
        None
    }
 }

-impl<N: Real, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S>
+impl<N: Complex, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S>
 where
    D: DimSub<U1>,                                   // For Hessenberg.
    ShapeConstraint: DimEq<Dynamic, DimDiff<D, U1>>, // For Hessenberg.
@ -514,7 +512,7 @@ where
    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
    /// continues indefinitely until convergence.
-    pub fn try_real_schur(self, eps: N, max_niter: usize) -> Option<RealSchur<N, D>> {
+    pub fn try_real_schur(self, eps: N::Real, max_niter: usize) -> Option<RealSchur<N, D>> {
        RealSchur::try_new(self.into_owned(), eps, max_niter)
    }

@ -546,7 +544,7 @@ where
        let schur = RealSchur::do_decompose(
            self.clone_owned(),
            &mut work,
-            N::default_epsilon(),
+            N::Real::default_epsilon(),
            0,
            false,
        )
@ -559,9 +557,10 @@ where
    }

    /// Computes the eigenvalues of this matrix.
-    pub fn complex_eigenvalues(&self) -> VectorN<Complex<N>, D>
+    pub fn complex_eigenvalues(&self) -> VectorN<NumComplex<N>, D>
    // FIXME: add balancing?
-    where DefaultAllocator: Allocator<Complex<N>, D> {
+    where N: Real,
+          DefaultAllocator: Allocator<NumComplex<N>, D> {
        let dim = self.data.shape().0;
        let mut work = unsafe { VectorN::new_uninitialized_generic(dim, U1) };

--- a/src/linalg/solve.rs
+++ b/src/linalg/solve.rs
@ -1,4 +1,4 @@
-use alga::general::Real;
+use alga::general::Complex;

 use base::allocator::Allocator;
 use base::constraint::{SameNumberOfRows, ShapeConstraint};
@ -6,7 +6,7 @@ use base::dimension::{Dim, U1};
 use base::storage::{Storage, StorageMut};
 use base::{DefaultAllocator, Matrix, MatrixMN, SquareMatrix, Vector};

-impl<N: Real, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
+impl<N: Complex, 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]
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@ -1,17 +1,19 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use num_complex::Complex;
+use num::Zero;
+use num_complex::Complex as NumComplex;
+use approx::AbsDiffEq;
 use std::ops::MulAssign;

-use alga::general::Real;
+use alga::general::Complex;
 use allocator::Allocator;
 use base::{DefaultAllocator, Matrix2, MatrixN, SquareMatrix, Vector2, VectorN};
 use dimension::{Dim, DimDiff, DimSub, U1, U2};
 use storage::Storage;

 use geometry::UnitComplex;
-use linalg::givens;
+use linalg::givens::GivensRotation;
 use linalg::SymmetricTridiagonal;

 /// Eigendecomposition of a symmetric matrix.
@ -35,7 +37,7 @@ use linalg::SymmetricTridiagonal;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct SymmetricEigen<N: Real, D: Dim>
+pub struct SymmetricEigen<N: Complex, D: Dim>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
 {
    /// The eigenvectors of the decomposed matrix.
@ -45,14 +47,14 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
    pub eigenvalues: VectorN<N, D>,
 }

-impl<N: Real, D: Dim> Copy for SymmetricEigen<N, D>
+impl<N: Complex, D: Dim> Copy for SymmetricEigen<N, D>
 where
    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
    MatrixN<N, D>: Copy,
    VectorN<N, D>: Copy,
 {}

-impl<N: Real, D: Dim> SymmetricEigen<N, D>
+impl<N: Complex, D: Dim> SymmetricEigen<N, D>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
 {
    /// Computes the eigendecomposition of the given symmetric matrix.
@ -63,7 +65,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
        D: DimSub<U1>,
        DefaultAllocator: Allocator<N, DimDiff<D, U1>>,
    {
-        Self::try_new(m, N::default_epsilon(), 0).unwrap()
+        Self::try_new(m, N::Real::default_epsilon(), 0).unwrap()
    }

    /// Computes the eigendecomposition of the given symmetric matrix with user-specified
@ -77,7 +79,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
    /// continues indefinitely until convergence.
-    pub fn try_new(m: MatrixN<N, D>, eps: N, max_niter: usize) -> Option<Self>
+    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>>,
@ -91,7 +93,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
    fn do_decompose(
        mut m: MatrixN<N, D>,
        eigenvectors: bool,
-        eps: N,
+        eps: N::Real,
        max_niter: usize,
    ) -> Option<(VectorN<N, D>, Option<MatrixN<N, D>>)>
    where
@ -103,11 +105,10 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
            "Unable to compute the eigendecomposition of a non-square matrix."
        );
        let dim = m.nrows();
-
-        let m_amax = m.amax();
+        let m_amax = m.camax();

        if !m_amax.is_zero() {
-            m /= m_amax;
+            m.unscale_mut(m_amax);
        }

        let (mut q, mut diag, mut off_diag);
@ -125,7 +126,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
        }

        if dim == 1 {
-            diag *= m_amax;
+            diag.scale_mut(m_amax);
            return Some((diag, q));
        }

@ -147,7 +148,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                for i in start..n {
                    let j = i + 1;

-                    if let Some((rot, norm)) = givens::cancel_y(&v) {
+                    if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
                        if i > start {
                            // Not the first iteration.
                            off_diag[i - 1] = norm;
@ -157,9 +158,9 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                        let mjj = diag[j];
                        let mij = off_diag[i];

-                        let cc = rot.cos_angle() * rot.cos_angle();
-                        let ss = rot.sin_angle() * rot.sin_angle();
-                        let cs = rot.cos_angle() * rot.sin_angle();
+                        let cc = rot.c() * rot.c();
+                        let ss = rot.s() * rot.s();
+                        let cs = rot.c() * rot.s();

                        let b = cs * ::convert(2.0) * mij;

@ -169,8 +170,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>

                        if i != n - 1 {
                            v.x = off_diag[i];
-                            v.y = -rot.sin_angle() * off_diag[i + 1];
-                            off_diag[i + 1] *= rot.cos_angle();
+                            v.y = -rot.s() * off_diag[i + 1];
+                            off_diag[i + 1] *= rot.c();
                        }

                        if let Some(ref mut q) = q {
@ -181,7 +182,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                    }
                }

-                if off_diag[m].abs() <= eps * (diag[m].abs() + diag[n].abs()) {
+                if off_diag[m].modulus() <= eps * (diag[m].modulus() + diag[n].modulus()) {
                    end -= 1;
                }
            } else if subdim == 2 {
@ -198,8 +199,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
                diag[start + 1] = eigvals[1];

                if let Some(ref mut q) = q {
-                    if let Some(basis) = basis.try_normalize(eps) {
-                        let rot = UnitComplex::new_unchecked(Complex::new(basis.x, basis.y));
+                    if let Some(rot) = GivensRotation::try_new(basis.x, basis.y, eps) {
                        rot.rotate_rows(&mut q.fixed_columns_mut::<U2>(start));
                    }
                }
@ -219,7 +219,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
            }
        }

-        diag *= m_amax;
+        diag.scale_mut(m_amax);

        Some((diag, q))
    }
@ -228,7 +228,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
        diag: &VectorN<N, D>,
        off_diag: &mut VectorN<N, DimDiff<D, U1>>,
        end: usize,
-        eps: N,
+        eps: N::Real,
    ) -> (usize, usize)
    where
        D: DimSub<U1>,
@ -239,7 +239,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
        while n > 0 {
            let m = n - 1;

-            if off_diag[m].abs() > eps * (diag[n].abs() + diag[m].abs()) {
+            if off_diag[m].modulus() > eps * (diag[n].modulus() + diag[m].modulus()) {
                break;
            }

@ -255,7 +255,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
            let m = new_start - 1;

            if off_diag[m].is_zero()
-                || off_diag[m].abs() <= eps * (diag[new_start].abs() + diag[m].abs())
+                || off_diag[m].modulus() <= eps * (diag[new_start].modulus() + diag[m].modulus())
            {
                off_diag[m] = N::zero();
                break;
@ -276,7 +276,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
            let val = self.eigenvalues[i];
            u_t.column_mut(i).mul_assign(val);
        }
-        u_t.transpose_mut();
+        u_t.conjugate_transpose_mut();
        &self.eigenvectors * u_t
    }
 }
@ -287,7 +287,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
 /// The inputs are interpreted as the 2x2 matrix:
 ///     tmm  tmn
 ///     tmn  tnn
-pub fn wilkinson_shift<N: Real>(tmm: N, tnn: N, tmn: N) -> N {
+pub fn wilkinson_shift<N: Complex>(tmm: N, tnn: N, tmn: N) -> N {
    let sq_tmn = tmn * tmn;
    if !sq_tmn.is_zero() {
        // We have the guarantee that the denominator won't be zero.
@ -303,7 +303,7 @@ pub fn wilkinson_shift<N: Real>(tmm: N, tnn: N, tmn: N) -> N {
 * Computations of eigenvalues for symmetric matrices.
 *
 */
-impl<N: Real, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
+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>>
 {
    /// Computes the eigendecomposition of this symmetric matrix.
@ -324,7 +324,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimD
    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
    /// continues indefinitely until convergence.
-    pub fn try_symmetric_eigen(self, eps: N, max_niter: usize) -> Option<SymmetricEigen<N, D>> {
+    pub fn try_symmetric_eigen(self, eps: N::Real, max_niter: usize) -> Option<SymmetricEigen<N, D>> {
        SymmetricEigen::try_new(self.into_owned(), eps, max_niter)
    }

@ -332,7 +332,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimD
    ///
    /// Only the lower-triangular part of the matrix is read.
    pub fn symmetric_eigenvalues(&self) -> VectorN<N, D> {
-        SymmetricEigen::do_decompose(self.clone_owned(), false, N::default_epsilon(), 0)
+        SymmetricEigen::do_decompose(self.clone_owned(), false, N::Real::default_epsilon(), 0)
            .unwrap()
            .0
    }
--- a/src/linalg/symmetric_tridiagonal.rs
+++ b/src/linalg/symmetric_tridiagonal.rs
@ -1,7 +1,7 @@
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};

-use alga::general::Real;
+use alga::general::Complex;
 use allocator::Allocator;
 use base::{DefaultAllocator, MatrixMN, MatrixN, SquareMatrix, VectorN};
 use dimension::{DimDiff, DimSub, U1};
@ -30,21 +30,21 @@ use linalg::householder;
    ))
 )]
 #[derive(Clone, Debug)]
-pub struct SymmetricTridiagonal<N: Real, D: DimSub<U1>>
+pub struct SymmetricTridiagonal<N: Complex, D: DimSub<U1>>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
 {
    tri: MatrixN<N, D>,
    off_diagonal: VectorN<N, DimDiff<D, U1>>,
 }

-impl<N: Real, D: DimSub<U1>> Copy for SymmetricTridiagonal<N, D>
+impl<N: Complex, D: DimSub<U1>> Copy for SymmetricTridiagonal<N, D>
 where
    DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
    MatrixN<N, D>: Copy,
    VectorN<N, DimDiff<D, U1>>: Copy,
 {}

-impl<N: Real, D: DimSub<U1>> SymmetricTridiagonal<N, D>
+impl<N: Complex, D: DimSub<U1>> SymmetricTridiagonal<N, D>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
 {
    /// Computes the tridiagonalization of the symmetric matrix `m`.
@ -75,17 +75,18 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
            if not_zero {
                let mut p = p.rows_range_mut(i..);

-                p.gemv_symm(::convert(2.0), &m, &axis, N::zero());
+                p.gemv_symm(::convert(2.0), &m, &axis.conjugate(), N::zero());
                let dot = axis.dot(&p);
-                p.axpy(-dot, &axis, N::one());
-                m.ger_symm(-N::one(), &p, &axis, N::one());
+//                p.axpy(-dot, &axis.conjugate(), N::one());
+                m.ger_symm(-N::one(), &p, &axis.conjugate(), N::one());
                m.ger_symm(-N::one(), &axis, &p, N::one());
+                m.ger_symm(dot * ::convert(2.0), &axis, &axis.conjugate(), N::one());
            }
        }

        Self {
            tri: m,
-            off_diagonal: off_diagonal,
+            off_diagonal,
        }
    }

@ -138,14 +139,14 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>

        for i in 0..self.off_diagonal.len() {
            self.tri[(i + 1, i)] = self.off_diagonal[i];
-            self.tri[(i, i + 1)] = self.off_diagonal[i];
+            self.tri[(i, i + 1)] = self.off_diagonal[i].conjugate();
        }

-        &q * self.tri * q.transpose()
+        &q * self.tri * q.conjugate_transpose()
    }
 }

-impl<N: Real, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
+impl<N: Complex, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
 where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
 {
    /// Computes the tridiagonalization of this symmetric matrix.
--- a/tests/core/helper.rs
+++ b/tests/core/helper.rs
@ -0,0 +1,54 @@
+// This module implement several methods to fill some
+// missing features of num-complex when it comes to randomness.
+
+use quickcheck::{Arbitrary, Gen};
+use rand::distributions::{Standard, Distribution};
+use rand::Rng;
+use num_complex::Complex;
+use na::Real;
+
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+pub struct RandComplex<N>(pub Complex<N>);
+
+impl<N: Arbitrary + Real> Arbitrary for RandComplex<N> {
+    #[inline]
+    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+        let im = Arbitrary::arbitrary(rng);
+        let re = Arbitrary::arbitrary(rng);
+        RandComplex(Complex::new(re, im))
+    }
+}
+
+impl<N: Real> Distribution<RandComplex<N>> for Standard
+    where
+        Standard: Distribution<N>,
+{
+    #[inline]
+    fn sample<'a, G: Rng + ?Sized>(&self, rng: &'a mut G) -> RandComplex<N> {
+        let re = rng.gen();
+        let im = rng.gen();
+        RandComplex(Complex::new(re, im))
+    }
+}
+
+// This is a wrapper similar to RandComplex, but for non-complex.
+// This exists only to make generic tests easier to write.
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+pub struct RandScalar<N>(pub N);
+
+impl<N: Arbitrary> Arbitrary for RandScalar<N> {
+    #[inline]
+    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+        RandScalar(Arbitrary::arbitrary(rng))
+    }
+}
+
+impl<N: Real> Distribution<RandScalar<N>> for Standard
+    where
+        Standard: Distribution<N>,
+{
+    #[inline]
+    fn sample<'a, G: Rng + ?Sized>(&self, rng: &'a mut G) -> RandScalar<N> {
+        RandScalar(self.sample(rng))
+    }
+}
--- a/tests/core/matrix.rs
+++ b/tests/core/matrix.rs
@ -1022,7 +1022,7 @@ mod finite_dim_inner_space_tests {
     *
     */
    #[cfg(feature = "arbitrary")]
-    fn is_subspace_basis<T: FiniteDimInnerSpace<Real = f64> + Display>(vs: &[T]) -> bool {
+    fn is_subspace_basis<T: FiniteDimInnerSpace<Real = f64, Complex = f64> + Display>(vs: &[T]) -> bool {
        for i in 0..vs.len() {
            // Basis elements must be normalized.
            if !relative_eq!(vs[i].norm(), 1.0, epsilon = 1.0e-7) {
--- a/tests/core/mod.rs
+++ b/tests/core/mod.rs
@ -8,3 +8,7 @@ mod matrix_slice;
 #[cfg(feature = "mint")]
 mod mint;
 mod serde;
+
+
+#[cfg(feature = "arbitrary")]
+pub mod helper;
--- a/tests/geometry/isometry.rs
+++ b/tests/geometry/isometry.rs
@ -82,7 +82,7 @@ quickcheck!(
        r: Rotation2<f64>,
        t: Translation2<f64>,
        v: Vector2<f64>,
-        p: Point2<f64>,
+        p: Point2<f64>
    ) -> bool
    {
        // (rotation × translation) * point = rotation × (translation * point)
@ -120,7 +120,7 @@ quickcheck!(
        r: Rotation3<f64>,
        t: Translation3<f64>,
        v: Vector3<f64>,
-        p: Point3<f64>,
+        p: Point3<f64>
    ) -> bool
    {
        // (rotation × translation) * point = rotation × (translation * point)
@ -158,7 +158,7 @@ quickcheck!(
        t: Translation3<f64>,
        v: Vector3<f64>,
        p: Point3<f64>,
-        r: Rotation3<f64>,
+        r: Rotation3<f64>
    ) -> bool
    {
        let iMi = i * i;
--- a/tests/geometry/similarity.rs
+++ b/tests/geometry/similarity.rs
@ -19,7 +19,7 @@ quickcheck!(
    fn inverse_is_parts_inversion(
        t: Translation3<f64>,
        r: UnitQuaternion<f64>,
-        scaling: f64,
+        scaling: f64
    ) -> bool
    {
        if relative_eq!(scaling, 0.0) {
@ -33,7 +33,7 @@ quickcheck!(
    fn multiply_equals_alga_transform(
        s: Similarity3<f64>,
        v: Vector3<f64>,
-        p: Point3<f64>,
+        p: Point3<f64>
    ) -> bool
    {
        s * v == s.transform_vector(&v)
@ -56,7 +56,7 @@ quickcheck!(
        t: Translation3<f64>,
        v: Vector3<f64>,
        p: Point3<f64>,
-        scaling: f64,
+        scaling: f64
    ) -> bool
    {
        if relative_eq!(scaling, 0.0) {
@ -152,7 +152,7 @@ quickcheck!(
        uq: UnitQuaternion<f64>,
        t: Translation3<f64>,
        v: Vector3<f64>,
-        p: Point3<f64>,
+        p: Point3<f64>
    ) -> bool
    {
        let sMs = s * s;
--- a/tests/geometry/unit_complex.rs
+++ b/tests/geometry/unit_complex.rs
@ -72,7 +72,7 @@ quickcheck!(
        uc: UnitComplex<f64>,
        v: Vector2<f64>,
        p: Point2<f64>,
-        r: Rotation2<f64>,
+        r: Rotation2<f64>
    ) -> bool
    {
        let uv = Unit::new_normalize(v);
--- a/tests/lib.rs
+++ b/tests/lib.rs
@ -12,9 +12,10 @@ extern crate num_traits as num;
 extern crate quickcheck;
 extern crate rand;
 extern crate serde_json;
+extern crate num_complex;

 mod core;
 mod geometry;
 mod linalg;
-#[cfg(feature = "sparse")]
-mod sparse;
+//#[cfg(feature = "sparse")]
+//mod sparse;
--- a/tests/linalg/bidiagonal.rs
+++ b/tests/linalg/bidiagonal.rs
@ -1,9 +1,11 @@
 #![cfg(feature = "arbitrary")]

 use na::{DMatrix, Matrix2, Matrix3x5, Matrix4, Matrix5x3};
+use core::helper::{RandScalar, RandComplex};

 quickcheck! {
-    fn bidiagonal(m: DMatrix<f64>) -> bool {
+    fn bidiagonal(m: DMatrix<RandComplex<f64>>) -> bool {
+        let m = m.map(|e| e.0);
        if m.len() == 0  {
            return true;
        }
@ -17,7 +19,8 @@ quickcheck! {
        relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
    }

-    fn bidiagonal_static_5_3(m: Matrix5x3<f64>) -> bool {
+    fn bidiagonal_static_5_3(m: Matrix5x3<RandComplex<f64>>) -> bool {
+        let m = m.map(|e| e.0);
        let bidiagonal = m.bidiagonalize();
        let (u, d, v_t) = bidiagonal.unpack();

@ -27,7 +30,8 @@ quickcheck! {
        relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
    }

-    fn bidiagonal_static_3_5(m: Matrix3x5<f64>) -> bool {
+    fn bidiagonal_static_3_5(m: Matrix3x5<RandComplex<f64>>) -> bool {
+        let m = m.map(|e| e.0);
        let bidiagonal = m.bidiagonalize();
        let (u, d, v_t) = bidiagonal.unpack();

@ -37,7 +41,8 @@ quickcheck! {
        relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
    }

-    fn bidiagonal_static_square(m: Matrix4<f64>) -> bool {
+    fn bidiagonal_static_square(m: Matrix4<RandComplex<f64>>) -> bool {
+        let m = m.map(|e| e.0);
        let bidiagonal = m.bidiagonalize();
        let (u, d, v_t) = bidiagonal.unpack();

@ -47,7 +52,8 @@ quickcheck! {
        relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
    }

-    fn bidiagonal_static_square_2x2(m: Matrix2<f64>) -> bool {
+    fn bidiagonal_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
+        let m = m.map(|e| e.0);
        let bidiagonal = m.bidiagonalize();
        let (u, d, v_t) = bidiagonal.unpack();

--- a/tests/linalg/eigen.rs
+++ b/tests/linalg/eigen.rs
@ -5,12 +5,13 @@ use na::DMatrix;
 #[cfg(feature = "arbitrary")]
 mod quickcheck_tests {
    use na::{DMatrix, Matrix2, Matrix3, Matrix4};
+    use core::helper::{RandScalar, RandComplex};
    use std::cmp;

    quickcheck! {
        fn symmetric_eigen(n: usize) -> bool {
            let n      = cmp::max(1, cmp::min(n, 10));
-            let m      = DMatrix::<f64>::new_random(n, n);
+            let m      = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);
            let eig    = m.clone().symmetric_eigen();
            let recomp = eig.recompose();

@ -21,9 +22,9 @@ mod quickcheck_tests {

        fn symmetric_eigen_singular(n: usize) -> bool {
            let n      = cmp::max(1, cmp::min(n, 10));
-            let mut m  = DMatrix::<f64>::new_random(n, n);
-            m.row_mut(n / 2).fill(0.0);
-            m.column_mut(n / 2).fill(0.0);
+            let mut m  = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);
+            m.row_mut(n / 2).fill(na::zero());
+            m.column_mut(n / 2).fill(na::zero());
            let eig    = m.clone().symmetric_eigen();
            let recomp = eig.recompose();

@ -32,7 +33,8 @@ mod quickcheck_tests {
            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
        }

-        fn symmetric_eigen_static_square_4x4(m: Matrix4<f64>) -> bool {
+        fn symmetric_eigen_static_square_4x4(m: Matrix4<RandComplex<f64>>) -> bool {
+            let m      = m.map(|e| e.0);
            let eig    = m.symmetric_eigen();
            let recomp = eig.recompose();

@ -41,7 +43,8 @@ mod quickcheck_tests {
            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
        }

-        fn symmetric_eigen_static_square_3x3(m: Matrix3<f64>) -> bool {
+        fn symmetric_eigen_static_square_3x3(m: Matrix3<RandComplex<f64>>) -> bool {
+            let m      = m.map(|e| e.0);
            let eig    = m.symmetric_eigen();
            let recomp = eig.recompose();

@ -50,7 +53,8 @@ mod quickcheck_tests {
            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
        }

-        fn symmetric_eigen_static_square_2x2(m: Matrix2<f64>) -> bool {
+        fn symmetric_eigen_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
+            let m      = m.map(|e| e.0);
            let eig    = m.symmetric_eigen();
            let recomp = eig.recompose();

--- a/tests/linalg/hessenberg.rs
+++ b/tests/linalg/hessenberg.rs
@ -1,6 +1,7 @@
 #![cfg(feature = "arbitrary")]

 use na::{DMatrix, Matrix2, Matrix4};
+use core::helper::{RandScalar, RandComplex};
 use std::cmp;

 #[test]
@ -14,20 +15,22 @@ fn hessenberg_simple() {
 quickcheck! {
    fn hessenberg(n: usize) -> bool {
        let n = cmp::max(1, cmp::min(n, 50));
-        let m = DMatrix::<f64>::new_random(n, n);
+        let m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);

        let hess = m.clone().hessenberg();
        let (p, h) = hess.unpack();
        relative_eq!(m, &p * h * p.transpose(), epsilon = 1.0e-7)
    }

-    fn hessenberg_static_mat2(m: Matrix2<f64>) -> bool {
+    fn hessenberg_static_mat2(m: Matrix2<RandComplex<f64>>) -> bool {
+        let m = m.map(|e| e.0);
        let hess = m.hessenberg();
        let (p, h) = hess.unpack();
        relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7)
    }

-    fn hessenberg_static(m: Matrix4<f64>) -> bool {
+    fn hessenberg_static(m: Matrix4<RandComplex<f64>>) -> bool {
+        let m = m.map(|e| e.0);
        let hess = m.hessenberg();
        let (p, h) = hess.unpack();
        relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7)
--- a/tests/linalg/lu.rs
+++ b/tests/linalg/lu.rs
@ -40,114 +40,134 @@ fn lu_simple_with_pivot() {

 #[cfg(feature = "arbitrary")]
 mod quickcheck_tests {
-    use std::cmp;
-    use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
+    use core::helper::{RandScalar, RandComplex};

-    quickcheck! {
-        fn lu(m: DMatrix<f64>) -> bool {
-            let mut m = m;
-            if m.len() == 0 {
-                m = DMatrix::new_random(1, 1);
-            }
+    macro_rules! gen_tests(
+        ($module: ident, $scalar: ty) => {
+            mod $module {
+                use std::cmp;
+                use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
+                #[allow(unused_imports)]
+                use core::helper::{RandScalar, RandComplex};

-            let lu = m.clone().lu();
-            let (p, l, u) = lu.unpack();
-            let mut lu = l * u;
-            p.inv_permute_rows(&mut lu);
+                quickcheck! {
+                    fn lu(m: DMatrix<$scalar>) -> bool {
+                        let mut m = m;
+                        if m.len() == 0 {
+                             m = DMatrix::<$scalar>::new_random(1, 1);
+                        }

-            relative_eq!(m, lu, epsilon = 1.0e-7)
-        }
+                        let m = m.map(|e| e.0);

-        fn lu_static_3_5(m: Matrix3x5<f64>) -> bool {
-            let lu = m.lu();
-            let (p, l, u) = lu.unpack();
-            let mut lu = l * u;
-            p.inv_permute_rows(&mut lu);
+                        let lu = m.clone().lu();
+                        let (p, l, u) = lu.unpack();
+                        let mut lu = l * u;
+                        p.inv_permute_rows(&mut lu);

-            relative_eq!(m, lu, epsilon = 1.0e-7)
-        }
+                        relative_eq!(m, lu, epsilon = 1.0e-7)
+                    }

-        fn lu_static_5_3(m: Matrix5x3<f64>) -> bool {
-            let lu = m.lu();
-            let (p, l, u) = lu.unpack();
-            let mut lu = l * u;
-            p.inv_permute_rows(&mut lu);
+                    fn lu_static_3_5(m: Matrix3x5<$scalar>) -> bool {
+                        let m = m.map(|e| e.0);
+                        let lu = m.lu();
+                        let (p, l, u) = lu.unpack();
+                        let mut lu = l * u;
+                        p.inv_permute_rows(&mut lu);

-            relative_eq!(m, lu, epsilon = 1.0e-7)
-        }
+                        relative_eq!(m, lu, epsilon = 1.0e-7)
+                    }

-        fn lu_static_square(m: Matrix4<f64>) -> bool {
-            let lu = m.lu();
-            let (p, l, u) = lu.unpack();
-            let mut lu = l * u;
-            p.inv_permute_rows(&mut lu);
+                    fn lu_static_5_3(m: Matrix5x3<$scalar>) -> bool {
+                        let m = m.map(|e| e.0);
+                        let lu = m.lu();
+                        let (p, l, u) = lu.unpack();
+                        let mut lu = l * u;
+                        p.inv_permute_rows(&mut lu);

-            relative_eq!(m, lu, epsilon = 1.0e-7)
-        }
+                        relative_eq!(m, lu, epsilon = 1.0e-7)
+                    }

-        fn lu_solve(n: usize, nb: usize) -> bool {
-            if n != 0 && nb != 0 {
-                let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
-                let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
-                let m  = DMatrix::<f64>::new_random(n, n);
+                    fn lu_static_square(m: Matrix4<$scalar>) -> bool {
+                        let m = m.map(|e| e.0);
+                        let lu = m.lu();
+                        let (p, l, u) = lu.unpack();
+                        let mut lu = l * u;
+                        p.inv_permute_rows(&mut lu);

-                let lu = m.clone().lu();
-                let b1 = DVector::new_random(n);
-                let b2 = DMatrix::new_random(n, nb);
+                        relative_eq!(m, lu, epsilon = 1.0e-7)
+                    }

-                let sol1 = lu.solve(&b1);
-                let sol2 = lu.solve(&b2);
+                    fn lu_solve(n: usize, nb: usize) -> bool {
+                        if n != 0 && nb != 0 {
+                            let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
+                            let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
+                            let m  = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);

-                return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
-                       (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
-            }
+                            let lu = m.clone().lu();
+                            let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
+                            let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);

-            return true;
-        }
+                            let sol1 = lu.solve(&b1);
+                            let sol2 = lu.solve(&b2);

-        fn lu_solve_static(m: Matrix4<f64>) -> bool {
-             let lu = m.lu();
-             let b1 = Vector4::new_random();
-             let b2 = Matrix4x3::new_random();
+                            return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
+                                   (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+                        }

-             let sol1 = lu.solve(&b1);
-             let sol2 = lu.solve(&b2);
+                        return true;
+                    }

-             return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
-                    (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
-        }
+                    fn lu_solve_static(m: Matrix4<$scalar>) -> bool {
+                         let m = m.map(|e| e.0);
+                         let lu = m.lu();
+                         let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
+                         let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);

-        fn lu_inverse(n: usize) -> bool {
-            let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
-            let m = DMatrix::<f64>::new_random(n, n);
+                         let sol1 = lu.solve(&b1);
+                         let sol2 = lu.solve(&b2);

-            let mut l = m.lower_triangle();
-            let mut u = m.upper_triangle();
+                         return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
+                                (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+                    }

-            // Ensure the matrix is well conditioned for inversion.
-            l.fill_diagonal(1.0);
-            u.fill_diagonal(1.0);
-            let m = l * u;
+                    fn lu_inverse(n: usize) -> bool {
+                        let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
+                        let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);

-            let m1  = m.clone().lu().try_inverse().unwrap();
-            let id1 = &m  * &m1;
-            let id2 = &m1 * &m;
+                        let mut l = m.lower_triangle();
+                        let mut u = m.upper_triangle();

-            return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
-        }
+                        // Ensure the matrix is well conditioned for inversion.
+                        l.fill_diagonal(na::one());
+                        u.fill_diagonal(na::one());
+                        let m = l * u;

-        fn lu_inverse_static(m: Matrix4<f64>) -> bool {
-            let lu  = m.lu();
+                        let m1  = m.clone().lu().try_inverse().unwrap();
+                        let id1 = &m  * &m1;
+                        let id2 = &m1 * &m;

-            if let Some(m1) = lu.try_inverse() {
-                let id1 = &m  * &m1;
-                let id2 = &m1 * &m;
+                        return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
+                    }

-                id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
-            }
-            else {
-                true
+                    fn lu_inverse_static(m: Matrix4<$scalar>) -> bool {
+                        let m = m.map(|e| e.0);
+                        let lu  = m.lu();
+
+                        if let Some(m1) = lu.try_inverse() {
+                            let id1 = &m  * &m1;
+                            let id2 = &m1 * &m;
+
+                            id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
+                        }
+                        else {
+                            true
+                        }
+                    }
+                }
            }
        }
-    }
+    );
+
+    gen_tests!(complex, RandComplex<f64>);
+    gen_tests!(f64, RandScalar<f64>);
 }
--- a/tests/linalg/qr.rs
+++ b/tests/linalg/qr.rs
@ -1,112 +1,133 @@
 #![cfg(feature = "arbitrary")]

-use na::{DMatrix, DVector, Matrix3x5, Matrix4, Matrix4x3, Matrix5x3, Vector4};
-use std::cmp;
+use core::helper::{RandScalar, RandComplex};

-quickcheck! {
-    fn qr(m: DMatrix<f64>) -> bool {
-        let qr = m.clone().qr();
-        let q  = qr.q();
-        let r  = qr.r();
+macro_rules! gen_tests(
+    ($module: ident, $scalar: ty) => {
+        mod $module {
+            use na::{DMatrix, DVector, Matrix3x5, Matrix4, Matrix4x3, Matrix5x3, Vector4};
+            use std::cmp;
+            use core::helper::{RandScalar, RandComplex};

-        relative_eq!(m, &q * r, epsilon = 1.0e-7) &&
-        q.is_orthogonal(1.0e-7)
-    }
+            quickcheck! {
+                fn qr(m: DMatrix<$scalar>) -> bool {
+                    let m = m.map(|e| e.0);
+                    let qr = m.clone().qr();
+                    let q  = qr.q();
+                    let r  = qr.r();

-    fn qr_static_5_3(m: Matrix5x3<f64>) -> bool {
-        let qr = m.qr();
-        let q  = qr.q();
-        let r  = qr.r();
+                    println!("m: {}", m);
+                    println!("qr: {}", &q * &r);

-        relative_eq!(m, q * r, epsilon = 1.0e-7) &&
-        q.is_orthogonal(1.0e-7)
-    }
+                    relative_eq!(m, &q * r, epsilon = 1.0e-7) &&
+                    q.is_orthogonal(1.0e-7)
+                }

-    fn qr_static_3_5(m: Matrix3x5<f64>) -> bool {
-        let qr = m.qr();
-        let q  = qr.q();
-        let r  = qr.r();
+                fn qr_static_5_3(m: Matrix5x3<$scalar>) -> bool {
+                    let m = m.map(|e| e.0);
+                    let qr = m.qr();
+                    let q  = qr.q();
+                    let r  = qr.r();

-        relative_eq!(m, q * r, epsilon = 1.0e-7) &&
-        q.is_orthogonal(1.0e-7)
-    }
+                    relative_eq!(m, q * r, epsilon = 1.0e-7) &&
+                    q.is_orthogonal(1.0e-7)
+                }

-    fn qr_static_square(m: Matrix4<f64>) -> bool {
-        let qr = m.qr();
-        let q  = qr.q();
-        let r  = qr.r();
+                fn qr_static_3_5(m: Matrix3x5<$scalar>) -> bool {
+                    let m = m.map(|e| e.0);
+                    let qr = m.qr();
+                    let q  = qr.q();
+                    let r  = qr.r();

-        println!("{}{}{}{}", q, r, q * r, m);
+                    relative_eq!(m, q * r, epsilon = 1.0e-7) &&
+                    q.is_orthogonal(1.0e-7)
+                }

-        relative_eq!(m, q * r, epsilon = 1.0e-7) &&
-        q.is_orthogonal(1.0e-7)
-    }
+                fn qr_static_square(m: Matrix4<$scalar>) -> bool {
+                    let m = m.map(|e| e.0);
+                    let qr = m.qr();
+                    let q  = qr.q();
+                    let r  = qr.r();

-    fn qr_solve(n: usize, nb: usize) -> bool {
-        if n != 0 && nb != 0 {
-            let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
-            let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
-            let m  = DMatrix::<f64>::new_random(n, n);
+                    println!("{}{}{}{}", q, r, q * r, m);

-            let qr = m.clone().qr();
-            let b1 = DVector::new_random(n);
-            let b2 = DMatrix::new_random(n, nb);
+                    relative_eq!(m, q * r, epsilon = 1.0e-7) &&
+                    q.is_orthogonal(1.0e-7)
+                }

-            if qr.is_invertible() {
-                let sol1 = qr.solve(&b1).unwrap();
-                let sol2 = qr.solve(&b2).unwrap();
+                fn qr_solve(n: usize, nb: usize) -> bool {
+                    if n != 0 && nb != 0 {
+                        let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
+                        let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
+                        let m  = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);

-                return relative_eq!(&m * sol1, b1, epsilon = 1.0e-6) &&
-                    relative_eq!(&m * sol2, b2, epsilon = 1.0e-6)
+                        let qr = m.clone().qr();
+                        let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
+                        let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
+
+                        if qr.is_invertible() {
+                            let sol1 = qr.solve(&b1).unwrap();
+                            let sol2 = qr.solve(&b2).unwrap();
+
+                            return relative_eq!(&m * sol1, b1, epsilon = 1.0e-6) &&
+                                relative_eq!(&m * sol2, b2, epsilon = 1.0e-6)
+                        }
+                    }
+
+                    return true;
+                }
+
+                fn qr_solve_static(m: Matrix4<$scalar>) -> bool {
+                     let m = m.map(|e| e.0);
+                     let qr = m.qr();
+                     let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
+                     let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
+
+                     if qr.is_invertible() {
+                         let sol1 = qr.solve(&b1).unwrap();
+                         let sol2 = qr.solve(&b2).unwrap();
+
+                         relative_eq!(m * sol1, b1, epsilon = 1.0e-6) &&
+                         relative_eq!(m * sol2, b2, epsilon = 1.0e-6)
+                     }
+                     else {
+                         false
+                     }
+                }
+
+                fn qr_inverse(n: usize) -> bool {
+                    let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
+                    let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
+
+                    if let Some(m1) = m.clone().qr().try_inverse() {
+                        let id1 = &m  * &m1;
+                        let id2 = &m1 * &m;
+
+                        id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
+                    }
+                    else {
+                        true
+                    }
+                }
+
+                fn qr_inverse_static(m: Matrix4<$scalar>) -> bool {
+                    let m  = m.map(|e| e.0);
+                    let qr = m.qr();
+
+                    if let Some(m1) = qr.try_inverse() {
+                        let id1 = &m  * &m1;
+                        let id2 = &m1 * &m;
+
+                        id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
+                    }
+                    else {
+                        true
+                    }
+                }
            }
        }
-
-        return true;
    }
+);

-    fn qr_solve_static(m: Matrix4<f64>) -> bool {
-         let qr = m.qr();
-         let b1 = Vector4::new_random();
-         let b2 = Matrix4x3::new_random();
-
-         if qr.is_invertible() {
-             let sol1 = qr.solve(&b1).unwrap();
-             let sol2 = qr.solve(&b2).unwrap();
-
-             relative_eq!(m * sol1, b1, epsilon = 1.0e-6) &&
-             relative_eq!(m * sol2, b2, epsilon = 1.0e-6)
-         }
-         else {
-             false
-         }
-    }
-
-    fn qr_inverse(n: usize) -> bool {
-        let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
-        let m = DMatrix::<f64>::new_random(n, n);
-
-        if let Some(m1) = m.clone().qr().try_inverse() {
-            let id1 = &m  * &m1;
-            let id2 = &m1 * &m;
-
-            id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
-        }
-        else {
-            true
-        }
-    }
-
-    fn qr_inverse_static(m: Matrix4<f64>) -> bool {
-        let qr  = m.qr();
-
-        if let Some(m1) = qr.try_inverse() {
-            let id1 = &m  * &m1;
-            let id2 = &m1 * &m;
-
-            id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
-        }
-        else {
-            true
-        }
-    }
-}
+gen_tests!(complex, RandComplex<f64>);
+gen_tests!(f64, RandScalar<f64>);
--- a/tests/linalg/tridiagonal.rs
+++ b/tests/linalg/tridiagonal.rs
@ -3,12 +3,24 @@
 use std::cmp;

 use na::{DMatrix, Matrix2, Matrix4};
+use core::helper::{RandScalar, RandComplex};

 quickcheck! {
-    fn symm_tridiagonal(n: usize) -> bool {
-        let n = cmp::max(1, cmp::min(n, 50));
-        let m = DMatrix::<f64>::new_random(n, n);
-        let tri = m.clone().symmetric_tridiagonalize();
+//    fn symm_tridiagonal(n: usize) -> bool {
+//        let n = cmp::max(1, cmp::min(n, 50));
+//        let m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
+//        let tri = m.clone().symmetric_tridiagonalize();
+//        let recomp = tri.recompose();
+//
+//        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
+//
+//        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();
+        println!("Internal tri: {}{}", tri.internal_tri(), tri.off_diagonal());
        let recomp = tri.recompose();

        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
@ -16,20 +28,11 @@ quickcheck! {
        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
    }

-    fn symm_tridiagonal_static_square(m: Matrix4<f64>) -> bool {
-        let tri = m.symmetric_tridiagonalize();
-        println!("{}{}", tri.internal_tri(), tri.off_diagonal());
-        let recomp = tri.recompose();
-
-        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
-
-        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
-    }
-
-    fn symm_tridiagonal_static_square_2x2(m: Matrix2<f64>) -> bool {
-        let tri = m.symmetric_tridiagonalize();
-        let recomp = tri.recompose();
-
-        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
-    }
+//    fn symm_tridiagonal_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
+//        let m = m.map(|e| e.0).hermitian_part();
+//        let tri = m.symmetric_tridiagonalize();
+//        let recomp = tri.recompose();
+//
+//        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
+//    }
 }