diff --git a/nalgebra-lapack/Cargo.toml b/nalgebra-lapack/Cargo.toml
index 66a613a9..b35a65aa 100644
--- a/nalgebra-lapack/Cargo.toml
+++ b/nalgebra-lapack/Cargo.toml
@@ -16,25 +16,24 @@ serde-serialize = [ "serde", "serde_derive" ]
 
 # For BLAS/LAPACK
 default    = ["openblas"]
-openblas   = ["lapack/openblas"]
-netlib     = ["lapack/netlib"]
-accelerate = ["lapack/accelerate"]
+openblas   = ["lapack-src/openblas"]
+netlib     = ["lapack-src/netlib"]
+accelerate = ["lapack-src/accelerate"]
+intel-mkl  = ["lapack-src/intel-mkl"]
 
 [dependencies]
 nalgebra      = { version = "0.14", path = ".." }
-num-traits    = "0.1"
-num-complex   = "0.1"
-alga          = "0.5"
-serde         = { version = "0.9", optional = true }
-serde_derive  = { version = "0.9", optional = true }
+num-traits    = "0.2"
+num-complex   = { version = "0.2.0-git", git = "https://github.com/rust-num/num-complex", default-features = false }
+alga          = { version = "0.5", git = "https://github.com/sebcrozet/alga.git", branch = "no_std", default-features = false }
+serde         = { version = "1.0", optional = true }
+serde_derive  = { version = "1.0", optional = true }
+lapack        = { version = "0.15", default-features = false }
+lapack-src    = { version = "0.1", default-features = false }
 # clippy = "*"
 
-[dependencies.lapack]
-version          = "0.11"
-default-features = false
-
 [dev-dependencies]
 nalgebra   = { version = "0.14", path = "..", features = [ "arbitrary" ] }
 quickcheck = "0.4"
-approx     = "0.1"
-rand       = "0.4"
+approx     = "0.2"
+rand       = "0.5"
diff --git a/nalgebra-lapack/src/cholesky.rs b/nalgebra-lapack/src/cholesky.rs
index d4ae2591..ba7830a5 100644
--- a/nalgebra-lapack/src/cholesky.rs
+++ b/nalgebra-lapack/src/cholesky.rs
@@ -4,21 +4,33 @@ use serde;
 use num::Zero;
 use num_complex::Complex;
 
-use na::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar};
+use na::allocator::Allocator;
 use na::dimension::Dim;
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar};
 
-use lapack::fortran as interface;
+use lapack;
 
 /// The cholesky decomposion of a symmetric-definite-positive matrix.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, D>,
-         MatrixN<N, D>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(deserialize = "DefaultAllocator: Allocator<N, D>,
-         MatrixN<N, D>: serde::Deserialize<'de>")))]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, D>,
+         MatrixN<N, D>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            deserialize = "DefaultAllocator: Allocator<N, D>,
+         MatrixN<N, D>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct Cholesky<N: Scalar, D: Dim>
 where
@@ -195,34 +207,24 @@ macro_rules! cholesky_scalar_impl(
         impl CholeskyScalar for $N {
             #[inline]
             fn xpotrf(uplo: u8, n: i32, a: &mut [Self], lda: i32, info: &mut i32) {
-                $xpotrf(uplo, n, a, lda, info)
+                unsafe { $xpotrf(uplo, n, a, lda, info) }
             }
 
             #[inline]
             fn xpotrs(uplo: u8, n: i32, nrhs: i32, a: &[Self], lda: i32,
                       b: &mut [Self], ldb: i32, info: &mut i32) {
-                $xpotrs(uplo, n, nrhs, a, lda, b, ldb, info)
+                unsafe { $xpotrs(uplo, n, nrhs, a, lda, b, ldb, info) }
             }
 
             #[inline]
             fn xpotri(uplo: u8, n: i32, a: &mut [Self], lda: i32, info: &mut i32) {
-                $xpotri(uplo, n, a, lda, info)
+                unsafe { $xpotri(uplo, n, a, lda, info) }
             }
         }
     )
 );
 
-cholesky_scalar_impl!(f32, interface::spotrf, interface::spotrs, interface::spotri);
-cholesky_scalar_impl!(f64, interface::dpotrf, interface::dpotrs, interface::dpotri);
-cholesky_scalar_impl!(
-    Complex<f32>,
-    interface::cpotrf,
-    interface::cpotrs,
-    interface::cpotri
-);
-cholesky_scalar_impl!(
-    Complex<f64>,
-    interface::zpotrf,
-    interface::zpotrs,
-    interface::zpotri
-);
+cholesky_scalar_impl!(f32, lapack::spotrf, lapack::spotrs, lapack::spotri);
+cholesky_scalar_impl!(f64, lapack::dpotrf, lapack::dpotrs, lapack::dpotri);
+cholesky_scalar_impl!(Complex<f32>, lapack::cpotrf, lapack::cpotrs, lapack::cpotri);
+cholesky_scalar_impl!(Complex<f64>, lapack::zpotrf, lapack::zpotrs, lapack::zpotri);
diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs
index b368d87e..e406f570 100644
--- a/nalgebra-lapack/src/eigen.rs
+++ b/nalgebra-lapack/src/eigen.rs
@@ -6,24 +6,36 @@ use num_complex::Complex;
 
 use alga::general::Real;
 
-use ComplexHelper;
-use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use na::allocator::Allocator;
 use na::dimension::{Dim, U1};
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use ComplexHelper;
 
-use lapack::fortran as interface;
+use lapack;
 
 /// Eigendecomposition of a real square matrix with real eigenvalues.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
          VectorN<N, D>: serde::Serialize,
-         MatrixN<N, D>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(deserialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+         MatrixN<N, D>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            deserialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
          VectorN<N, D>: serde::Serialize,
-         MatrixN<N, D>: serde::Deserialize<'de>")))]
+         MatrixN<N, D>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct Eigen<N: Scalar, D: Dim>
 where
@@ -350,7 +362,7 @@ macro_rules! real_eigensystem_scalar_impl (
                      wr: &mut [Self], wi: &mut [Self],
                      vl: &mut [Self], ldvl: i32, vr: &mut [Self], ldvr: i32,
                      work: &mut [Self], lwork: i32, info: &mut i32) {
-                $xgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
+                unsafe { $xgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info) }
             }
 
 
@@ -361,16 +373,16 @@ macro_rules! real_eigensystem_scalar_impl (
                 let mut work = [ Zero::zero() ];
                 let lwork = -1 as i32;
 
-                $xgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, &mut work, lwork, info);
+                unsafe { $xgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, &mut work, lwork, info) };
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
     )
 );
 
-real_eigensystem_scalar_impl!(f32, interface::sgeev);
-real_eigensystem_scalar_impl!(f64, interface::dgeev);
+real_eigensystem_scalar_impl!(f32, lapack::sgeev);
+real_eigensystem_scalar_impl!(f64, lapack::dgeev);
 
 //// FIXME: decomposition of complex matrix and matrices with complex eigenvalues.
-// eigensystem_complex_impl!(f32, interface::cgeev);
-// eigensystem_complex_impl!(f64, interface::zgeev);
+// eigensystem_complex_impl!(f32, lapack::cgeev);
+// eigensystem_complex_impl!(f64, lapack::zgeev);
diff --git a/nalgebra-lapack/src/hessenberg.rs b/nalgebra-lapack/src/hessenberg.rs
index 2dae1e1e..6c1925f4 100644
--- a/nalgebra-lapack/src/hessenberg.rs
+++ b/nalgebra-lapack/src/hessenberg.rs
@@ -1,26 +1,38 @@
 use num::Zero;
 use num_complex::Complex;
 
-use ComplexHelper;
-use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use na::allocator::Allocator;
 use na::dimension::{DimDiff, DimSub, U1};
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use ComplexHelper;
 
-use lapack::fortran as interface;
+use lapack;
 
 /// The Hessenberg decomposition of a general matrix.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, D, D> +
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, D, D> +
                            Allocator<N, DimDiff<D, U1>>,
          MatrixN<N, D>: serde::Serialize,
-         VectorN<N, DimDiff<D, U1>>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(deserialize = "DefaultAllocator: Allocator<N, D, D> +
+         VectorN<N, DimDiff<D, U1>>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            deserialize = "DefaultAllocator: Allocator<N, D, D> +
                            Allocator<N, DimDiff<D, U1>>,
          MatrixN<N, D>: serde::Deserialize<'de>,
-         VectorN<N, DimDiff<D, U1>>: serde::Deserialize<'de>")))]
+         VectorN<N, DimDiff<D, U1>>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct Hessenberg<N: Scalar, D: DimSub<U1>>
 where
@@ -188,7 +200,7 @@ macro_rules! hessenberg_scalar_impl(
             #[inline]
             fn xgehrd(n: i32, ilo: i32, ihi: i32, a: &mut [Self], lda: i32,
                       tau: &mut [Self], work: &mut [Self], lwork: i32, info: &mut i32) {
-                $xgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
+                unsafe { $xgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info) }
             }
 
             #[inline]
@@ -197,7 +209,7 @@ macro_rules! hessenberg_scalar_impl(
                 let mut work = [ Zero::zero() ];
                 let lwork = -1 as i32;
 
-                $xgehrd(n, ilo, ihi, a, lda, tau, &mut work, lwork, info);
+                unsafe { $xgehrd(n, ilo, ihi, a, lda, tau, &mut work, lwork, info) };
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
@@ -210,7 +222,7 @@ macro_rules! hessenberg_real_impl(
             #[inline]
             fn xorghr(n: i32, ilo: i32, ihi: i32, a: &mut [Self], lda: i32, tau: &[Self],
                       work: &mut [Self], lwork: i32, info: &mut i32) {
-                $xorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
+                unsafe { $xorghr(n, ilo, ihi, a, lda, tau, work, lwork, info) }
             }
 
             #[inline]
@@ -219,17 +231,17 @@ macro_rules! hessenberg_real_impl(
                 let mut work = [ Zero::zero() ];
                 let lwork = -1 as i32;
 
-                $xorghr(n, ilo, ihi, a, lda, tau, &mut work, lwork, info);
+                unsafe { $xorghr(n, ilo, ihi, a, lda, tau, &mut work, lwork, info) };
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
     )
 );
 
-hessenberg_scalar_impl!(f32, interface::sgehrd);
-hessenberg_scalar_impl!(f64, interface::dgehrd);
-hessenberg_scalar_impl!(Complex<f32>, interface::cgehrd);
-hessenberg_scalar_impl!(Complex<f64>, interface::zgehrd);
+hessenberg_scalar_impl!(f32, lapack::sgehrd);
+hessenberg_scalar_impl!(f64, lapack::dgehrd);
+hessenberg_scalar_impl!(Complex<f32>, lapack::cgehrd);
+hessenberg_scalar_impl!(Complex<f64>, lapack::zgehrd);
 
-hessenberg_real_impl!(f32, interface::sorghr);
-hessenberg_real_impl!(f64, interface::dorghr);
+hessenberg_real_impl!(f32, lapack::sorghr);
+hessenberg_real_impl!(f64, lapack::dorghr);
diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs
index d0c905f5..3af37843 100644
--- a/nalgebra-lapack/src/lib.rs
+++ b/nalgebra-lapack/src/lib.rs
@@ -72,30 +72,32 @@
 
 extern crate alga;
 extern crate lapack;
+extern crate lapack_src;
 extern crate nalgebra as na;
 extern crate num_complex;
 extern crate num_traits as num;
 
 mod lapack_check;
-mod svd;
-mod eigen;
-mod symmetric_eigen;
+
 mod cholesky;
+mod eigen;
+mod hessenberg;
 mod lu;
 mod qr;
-mod hessenberg;
 mod schur;
+mod svd;
+mod symmetric_eigen;
 
 use num_complex::Complex;
 
-pub use self::svd::SVD;
 pub use self::cholesky::{Cholesky, CholeskyScalar};
-pub use self::lu::{LUScalar, LU};
 pub use self::eigen::Eigen;
-pub use self::symmetric_eigen::SymmetricEigen;
-pub use self::qr::QR;
 pub use self::hessenberg::Hessenberg;
+pub use self::lu::{LUScalar, LU};
+pub use self::qr::QR;
 pub use self::schur::RealSchur;
+pub use self::svd::SVD;
+pub use self::symmetric_eigen::SymmetricEigen;
 
 trait ComplexHelper {
     type RealPart;
diff --git a/nalgebra-lapack/src/lu.rs b/nalgebra-lapack/src/lu.rs
index 1ba9e3ba..4d3f1d94 100644
--- a/nalgebra-lapack/src/lu.rs
+++ b/nalgebra-lapack/src/lu.rs
@@ -1,13 +1,13 @@
 use num::{One, Zero};
 use num_complex::Complex;
 
-use ComplexHelper;
-use na::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, VectorN};
+use na::allocator::Allocator;
 use na::dimension::{Dim, DimMin, DimMinimum, U1};
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, VectorN};
+use ComplexHelper;
 
-use lapack::fortran as interface;
+use lapack;
 
 /// LU decomposition with partial pivoting.
 ///
@@ -18,16 +18,28 @@ use lapack::fortran as interface;
 ///
 /// Those are such that `M == P * L * U`.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, R, C> +
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, R, C> +
                            Allocator<i32, DimMinimum<R, C>>,
          MatrixMN<N, R, C>: serde::Serialize,
-         PermutationSequence<DimMinimum<R, C>>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(deserialize = "DefaultAllocator: Allocator<N, R, C> +
+         PermutationSequence<DimMinimum<R, C>>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            deserialize = "DefaultAllocator: Allocator<N, R, C> +
                            Allocator<i32, DimMinimum<R, C>>,
          MatrixMN<N, R, C>: serde::Deserialize<'de>,
-         PermutationSequence<DimMinimum<R, C>>: serde::Deserialize<'de>")))]
+         PermutationSequence<DimMinimum<R, C>>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct LU<N: Scalar, R: DimMin<C>, C: Dim>
 where
@@ -344,24 +356,24 @@ macro_rules! lup_scalar_impl(
         impl LUScalar for $N {
             #[inline]
             fn xgetrf(m: i32, n: i32, a: &mut [Self], lda: i32, ipiv: &mut [i32], info: &mut i32) {
-                $xgetrf(m, n, a, lda, ipiv, info)
+                unsafe { $xgetrf(m, n, a, lda, ipiv, info) }
             }
 
             #[inline]
             fn xlaswp(n: i32, a: &mut [Self], lda: i32, k1: i32, k2: i32, ipiv: &[i32], incx: i32) {
-                $xlaswp(n, a, lda, k1, k2, ipiv, incx)
+                unsafe { $xlaswp(n, a, lda, k1, k2, ipiv, incx) }
             }
 
             #[inline]
             fn xgetrs(trans: u8, n: i32, nrhs: i32, a: &[Self], lda: i32, ipiv: &[i32],
                       b: &mut [Self], ldb: i32, info: &mut i32) {
-                $xgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
+                unsafe { $xgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info) }
             }
 
             #[inline]
             fn xgetri(n: i32, a: &mut [Self], lda: i32, ipiv: &[i32],
                       work: &mut [Self], lwork: i32, info: &mut i32) {
-                $xgetri(n, a, lda, ipiv, work, lwork, info)
+                unsafe { $xgetri(n, a, lda, ipiv, work, lwork, info) }
             }
 
             #[inline]
@@ -369,7 +381,7 @@ macro_rules! lup_scalar_impl(
                 let mut work = [ Zero::zero() ];
                 let lwork = -1 as i32;
 
-                $xgetri(n, a, lda, ipiv, &mut work, lwork, info);
+                unsafe { $xgetri(n, a, lda, ipiv, &mut work, lwork, info); }
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
@@ -378,29 +390,29 @@ macro_rules! lup_scalar_impl(
 
 lup_scalar_impl!(
     f32,
-    interface::sgetrf,
-    interface::slaswp,
-    interface::sgetrs,
-    interface::sgetri
+    lapack::sgetrf,
+    lapack::slaswp,
+    lapack::sgetrs,
+    lapack::sgetri
 );
 lup_scalar_impl!(
     f64,
-    interface::dgetrf,
-    interface::dlaswp,
-    interface::dgetrs,
-    interface::dgetri
+    lapack::dgetrf,
+    lapack::dlaswp,
+    lapack::dgetrs,
+    lapack::dgetri
 );
 lup_scalar_impl!(
     Complex<f32>,
-    interface::cgetrf,
-    interface::claswp,
-    interface::cgetrs,
-    interface::cgetri
+    lapack::cgetrf,
+    lapack::claswp,
+    lapack::cgetrs,
+    lapack::cgetri
 );
 lup_scalar_impl!(
     Complex<f64>,
-    interface::zgetrf,
-    interface::zlaswp,
-    interface::zgetrs,
-    interface::zgetri
+    lapack::zgetrf,
+    lapack::zlaswp,
+    lapack::zgetrs,
+    lapack::zgetri
 );
diff --git a/nalgebra-lapack/src/qr.rs b/nalgebra-lapack/src/qr.rs
index df4e139e..128a2f03 100644
--- a/nalgebra-lapack/src/qr.rs
+++ b/nalgebra-lapack/src/qr.rs
@@ -1,29 +1,41 @@
 #[cfg(feature = "serde-serialize")]
 use serde;
 
-use num_complex::Complex;
 use num::Zero;
+use num_complex::Complex;
 
-use ComplexHelper;
-use na::{DefaultAllocator, Matrix, MatrixMN, Scalar, VectorN};
+use na::allocator::Allocator;
 use na::dimension::{Dim, DimMin, DimMinimum, U1};
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixMN, Scalar, VectorN};
+use ComplexHelper;
 
-use lapack::fortran as interface;
+use lapack;
 
 /// The QR decomposition of a general matrix.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, R, C> +
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, R, C> +
                            Allocator<N, DimMinimum<R, C>>,
          MatrixMN<N, R, C>: serde::Serialize,
-         VectorN<N, DimMinimum<R, C>>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(deserialize = "DefaultAllocator: Allocator<N, R, C> +
+         VectorN<N, DimMinimum<R, C>>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            deserialize = "DefaultAllocator: Allocator<N, R, C> +
                            Allocator<N, DimMinimum<R, C>>,
          MatrixMN<N, R, C>: serde::Deserialize<'de>,
-         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>")))]
+         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct QR<N: Scalar, R: DimMin<C>, C: Dim>
 where
@@ -217,7 +229,7 @@ macro_rules! qr_scalar_impl(
             #[inline]
             fn xgeqrf(m: i32, n: i32, a: &mut [Self], lda: i32, tau: &mut [Self],
                       work: &mut [Self], lwork: i32, info: &mut i32) {
-                $xgeqrf(m, n, a, lda, tau, work, lwork, info)
+                unsafe { $xgeqrf(m, n, a, lda, tau, work, lwork, info) }
             }
 
             #[inline]
@@ -226,7 +238,7 @@ macro_rules! qr_scalar_impl(
                 let mut work = [ Zero::zero() ];
                 let lwork = -1 as i32;
 
-                $xgeqrf(m, n, a, lda, tau, &mut work, lwork, info);
+                unsafe { $xgeqrf(m, n, a, lda, tau, &mut work, lwork, info); }
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
@@ -239,7 +251,7 @@ macro_rules! qr_real_impl(
             #[inline]
             fn xorgqr(m: i32, n: i32, k: i32, a: &mut [Self], lda: i32, tau: &[Self],
                       work: &mut [Self], lwork: i32, info: &mut i32) {
-                $xorgqr(m, n, k, a, lda, tau, work, lwork, info)
+                unsafe { $xorgqr(m, n, k, a, lda, tau, work, lwork, info) }
             }
 
             #[inline]
@@ -248,17 +260,17 @@ macro_rules! qr_real_impl(
                 let mut work = [ Zero::zero() ];
                 let lwork = -1 as i32;
 
-                $xorgqr(m, n, k, a, lda, tau, &mut work, lwork, info);
+                unsafe { $xorgqr(m, n, k, a, lda, tau, &mut work, lwork, info); }
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
     )
 );
 
-qr_scalar_impl!(f32, interface::sgeqrf);
-qr_scalar_impl!(f64, interface::dgeqrf);
-qr_scalar_impl!(Complex<f32>, interface::cgeqrf);
-qr_scalar_impl!(Complex<f64>, interface::zgeqrf);
+qr_scalar_impl!(f32, lapack::sgeqrf);
+qr_scalar_impl!(f64, lapack::dgeqrf);
+qr_scalar_impl!(Complex<f32>, lapack::cgeqrf);
+qr_scalar_impl!(Complex<f64>, lapack::zgeqrf);
 
-qr_real_impl!(f32, interface::sorgqr);
-qr_real_impl!(f64, interface::dorgqr);
+qr_real_impl!(f32, lapack::sorgqr);
+qr_real_impl!(f64, lapack::dorgqr);
diff --git a/nalgebra-lapack/src/schur.rs b/nalgebra-lapack/src/schur.rs
index 9cc8cf90..eccc60d6 100644
--- a/nalgebra-lapack/src/schur.rs
+++ b/nalgebra-lapack/src/schur.rs
@@ -6,24 +6,36 @@ use num_complex::Complex;
 
 use alga::general::Real;
 
-use ComplexHelper;
-use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use na::allocator::Allocator;
 use na::dimension::{Dim, U1};
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use ComplexHelper;
 
-use lapack::fortran as interface;
+use lapack;
 
 /// Eigendecomposition of a real square matrix with real eigenvalues.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
          VectorN<N, D>: serde::Serialize,
-         MatrixN<N, D>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(deserialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
+         MatrixN<N, D>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            deserialize = "DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
          VectorN<N, D>: serde::Serialize,
-         MatrixN<N, D>: serde::Deserialize<'de>")))]
+         MatrixN<N, D>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct RealSchur<N: Scalar, D: Dim>
 where
@@ -216,7 +228,7 @@ macro_rules! real_eigensystem_scalar_impl (
                      lwork:  i32,
                      bwork:  &mut [i32],
                      info:   &mut i32) {
-                $xgees(jobvs, sort, None, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info);
+                unsafe { $xgees(jobvs, sort, None, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info); }
             }
 
 
@@ -238,12 +250,12 @@ macro_rules! real_eigensystem_scalar_impl (
                 let mut work = [ Zero::zero() ];
                 let lwork    = -1 as i32;
 
-                $xgees(jobvs, sort, None, n, a, lda, sdim, wr, wi, vs, ldvs, &mut work, lwork, bwork, info);
+                unsafe { $xgees(jobvs, sort, None, n, a, lda, sdim, wr, wi, vs, ldvs, &mut work, lwork, bwork, info); }
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
     )
 );
 
-real_eigensystem_scalar_impl!(f32, interface::sgees);
-real_eigensystem_scalar_impl!(f64, interface::dgees);
+real_eigensystem_scalar_impl!(f32, lapack::sgees);
+real_eigensystem_scalar_impl!(f64, lapack::dgees);
diff --git a/nalgebra-lapack/src/svd.rs b/nalgebra-lapack/src/svd.rs
index 61ec218e..ebfe1aae 100644
--- a/nalgebra-lapack/src/svd.rs
+++ b/nalgebra-lapack/src/svd.rs
@@ -1,32 +1,44 @@
 #[cfg(feature = "serde-serialize")]
 use serde;
 
-use std::cmp;
 use num::Signed;
+use std::cmp;
 
-use na::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, VectorN};
+use na::allocator::Allocator;
 use na::dimension::{Dim, DimMin, DimMinimum, U1};
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, VectorN};
 
-use lapack::fortran as interface;
+use lapack;
 
 /// The SVD decomposition of a general matrix.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
                            Allocator<N, R, R> +
                            Allocator<N, C, C>,
          MatrixN<N, R>: serde::Serialize,
          MatrixN<N, C>: serde::Serialize,
-         VectorN<N, DimMinimum<R, C>>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
+         VectorN<N, DimMinimum<R, C>>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
                            Allocator<N, R, R> +
                            Allocator<N, C, C>,
          MatrixN<N, R>: serde::Deserialize<'de>,
          MatrixN<N, C>: serde::Deserialize<'de>,
-         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>")))]
+         VectorN<N, DimMinimum<R, C>>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct SVD<N: Scalar, R: DimMin<C>, C: Dim>
 where
@@ -107,17 +119,22 @@ macro_rules! svd_impl(
                 let mut info  = 0;
                 let mut iwork = unsafe { ::uninitialized_vec(8 * cmp::min(nrows.value(), ncols.value())) };
 
-                $lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(),
+                unsafe {
+                    $lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(),
                     lda, &mut s.as_mut_slice(), u.as_mut_slice(), ldu as i32, vt.as_mut_slice(),
                     ldvt as i32, &mut work, lwork, &mut iwork, &mut info);
+                }
                 lapack_check!(info);
 
                 lwork = work[0] as i32;
                 let mut work = unsafe { ::uninitialized_vec(lwork as usize) };
 
+                unsafe {
                 $lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(),
                     lda, &mut s.as_mut_slice(), u.as_mut_slice(), ldu as i32, vt.as_mut_slice(),
                     ldvt as i32, &mut work, lwork, &mut iwork, &mut info);
+                }
+
                 lapack_check!(info);
 
                 Some(SVD { u: u, singular_values: s, vt: vt })
@@ -274,7 +291,7 @@ macro_rules! svd_complex_impl(
 );
 */
 
-svd_impl!(f32, interface::sgesdd);
-svd_impl!(f64, interface::dgesdd);
-// svd_complex_impl!(lapack_svd_complex_f32, f32, interface::cgesvd);
-// svd_complex_impl!(lapack_svd_complex_f64, f64, interface::zgesvd);
+svd_impl!(f32, lapack::sgesdd);
+svd_impl!(f64, lapack::dgesdd);
+// svd_complex_impl!(lapack_svd_complex_f32, f32, lapack::cgesvd);
+// svd_complex_impl!(lapack_svd_complex_f64, f64, lapack::zgesvd);
diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs
index 8b21511d..bf041bbe 100644
--- a/nalgebra-lapack/src/symmetric_eigen.rs
+++ b/nalgebra-lapack/src/symmetric_eigen.rs
@@ -6,26 +6,38 @@ use std::ops::MulAssign;
 
 use alga::general::Real;
 
-use ComplexHelper;
-use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use na::allocator::Allocator;
 use na::dimension::{Dim, U1};
 use na::storage::Storage;
-use na::allocator::Allocator;
+use na::{DefaultAllocator, Matrix, MatrixN, Scalar, VectorN};
+use ComplexHelper;
 
-use lapack::fortran as interface;
+use lapack;
 
 /// Eigendecomposition of a real square symmetric matrix with real eigenvalues.
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(serialize = "DefaultAllocator: Allocator<N, D, D> +
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            serialize = "DefaultAllocator: Allocator<N, D, D> +
                            Allocator<N, D>,
          VectorN<N, D>: serde::Serialize,
-         MatrixN<N, D>: serde::Serialize")))]
-#[cfg_attr(feature = "serde-serialize",
-           serde(bound(deserialize = "DefaultAllocator: Allocator<N, D, D> +
+         MatrixN<N, D>: serde::Serialize"
+        )
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(
+            deserialize = "DefaultAllocator: Allocator<N, D, D> +
                            Allocator<N, D>,
          VectorN<N, D>: serde::Deserialize<'de>,
-         MatrixN<N, D>: serde::Deserialize<'de>")))]
+         MatrixN<N, D>: serde::Deserialize<'de>"
+        )
+    )
+)]
 #[derive(Clone, Debug)]
 pub struct SymmetricEigen<N: Scalar, D: Dim>
 where
@@ -187,7 +199,7 @@ macro_rules! real_eigensystem_scalar_impl (
             #[inline]
             fn xsyev(jobz: u8, uplo: u8, n: i32, a: &mut [Self], lda: i32, w: &mut [Self], work: &mut [Self],
                      lwork: i32, info: &mut i32) {
-                $xsyev(jobz, uplo, n, a, lda, w, work, lwork, info)
+                unsafe { $xsyev(jobz, uplo, n, a, lda, w, work, lwork, info) }
             }
 
 
@@ -197,12 +209,12 @@ macro_rules! real_eigensystem_scalar_impl (
                 let mut w    = [ Zero::zero() ];
                 let lwork    = -1 as i32;
 
-                $xsyev(jobz, uplo, n, a, lda, &mut w, &mut work, lwork, info);
+                unsafe { $xsyev(jobz, uplo, n, a, lda, &mut w, &mut work, lwork, info); }
                 ComplexHelper::real_part(work[0]) as i32
             }
         }
     )
 );
 
-real_eigensystem_scalar_impl!(f32, interface::ssyev);
-real_eigensystem_scalar_impl!(f64, interface::dsyev);
+real_eigensystem_scalar_impl!(f32, lapack::ssyev);
+real_eigensystem_scalar_impl!(f64, lapack::dsyev);