diff --git a/.github/workflows/codequality.yml b/.github/workflows/codequality.yml
index 23401df0e..c32ff518b 100644
--- a/.github/workflows/codequality.yml
+++ b/.github/workflows/codequality.yml
@@ -36,11 +36,14 @@ jobs:
         uses: actions-rs/cargo@v1
         with:
           command: clippy
-          args: --all-targets --all-features -- -D warnings
+          args: --all-targets -- -D warnings
 
   coverage:
     name: coverage
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-18.04
+    container:
+      image: rustmath/mkl-rust:1.43.0
+      options: --security-opt seccomp=unconfined
 
     steps:
       - name: Checkout sources
@@ -72,7 +75,7 @@ jobs:
 
       - name: Generate code coverage
         run: |
-          cargo +nightly tarpaulin --verbose --all-features --timeout 120 --out Xml
+          cargo +nightly tarpaulin --verbose --features intel-mkl-system --timeout 120 --out Xml --all --release
 
       - name: Upload to codecov.io
         uses: codecov/codecov-action@v1
diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index 643fe5bc8..b0fdf4b87 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -5,7 +5,10 @@ name: Run Tests
 jobs:
   testing:
     name: testing
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-18.04
+    container:
+      image: rustmath/mkl-rust:1.43.0
+      options: --security-opt seccomp=unconfined
     strategy:
       matrix:
         toolchain:
@@ -15,9 +18,6 @@ jobs:
       - name: Checkout sources
         uses: actions/checkout@v2
 
-      - name: Install ubuntu packages
-        run: sudo apt-get install libssl-dev gfortran libopenblas-dev liblapack-dev liblapacke-dev libatlas-base-dev
-
       - name: Install toolchain
         uses: actions-rs/toolchain@v1
         with:
@@ -29,6 +29,4 @@ jobs:
         uses: actions-rs/cargo@v1
         with:
           command: test
-          args: --all --release
-        env:
-          RUSTFLAGS: -C link-arg=-llapacke -C link-arg=-lcblas
+          args: --all --release --features intel-mkl-system
diff --git a/Cargo.toml b/Cargo.toml
index 9ccd50919..8671790c1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -18,10 +18,40 @@ categories = ["algorithms", "mathematics", "science"]
 
 exclude = [".github/"]
 
+[features]
+default = []
+
+netlib-static = ["ndarray-linalg", "netlib-src/static"]
+netlib-system = ["ndarray-linalg", "netlib-src/system"]
+
+openblas-static = ["ndarray-linalg", "openblas-src/static"]
+openblas-system = ["ndarray-linalg", "openblas-src/system"]
+
+intel-mkl-static = ["ndarray-linalg", "intel-mkl-src/mkl-static-lp64-seq", "intel-mkl-src/download"]
+intel-mkl-system = ["ndarray-linalg", "intel-mkl-src/mkl-dynamic-lp64-seq"]
+
 [dependencies]
 num-traits = "0.2"
 rand = "0.7"
-ndarray = { version = "0.13", default-features = false }
+ndarray = { version = "0.13", default-features = false, features = ["blas"] }
+ndarray-linalg = { version = "0.12", optional = true }
+
+[dependencies.intel-mkl-src]
+version = "0.6.0"
+default-features = false
+optional = true
+
+[dependencies.netlib-src]
+version = "0.8.0"
+optional = true
+features = ["cblas"]
+default-features = false
+
+[dependencies.openblas-src]
+version = "0.9.0"
+optional = true
+default-features = false
+features = ["cblas"]
 
 [dev-dependencies]
 ndarray-rand = "0.12"
diff --git a/README.md b/README.md
index d61431844..9fb1ff0b6 100644
--- a/README.md
+++ b/README.md
@@ -41,6 +41,22 @@ We believe that only a significant community effort can nurture, build, and sust
 
 If this strikes a chord with you, please take a look at the [roadmap](https://github.com/rust-ml/linfa/issues/7) and get involved!
 
+## BLAS/Lapack backend
+
+At the moment you can choose between the following BLAS/LAPACK backends: `openblas`, `netblas` or `intel-mkl`
+
+|Backend  | Linux | Windows | macOS |
+|:--------|:-----:|:-------:|:-----:|
+|OpenBLAS |✔️      |-        |-      |
+|Netlib   |✔️      |-        |-      |
+|Intel MKL|✔️      |✔️        |✔️      |
+
+For example if you want to use the system IntelMKL library for the PCA example, then pass the corresponding feature:
+```
+cd linfa-reduction && cargo run --release --example pca --features linfa/intel-mkl-system
+```
+This selects the `intel-mkl` system library as BLAS/LAPACK backend. On the other hand if you want to compile the library and link it with the generated artifacts, pass `intel-mkl-static`.
+
 # License
 Dual-licensed to be compatible with the Rust project.
 
diff --git a/linfa-clustering/Cargo.toml b/linfa-clustering/Cargo.toml
index d1422a95d..5c5e68762 100644
--- a/linfa-clustering/Cargo.toml
+++ b/linfa-clustering/Cargo.toml
@@ -40,7 +40,6 @@ linfa = { version = "0.2.1", path = ".." }
 partitions = "0.2.4"
 
 [dev-dependencies]
-openblas-src = { version = "0.9", default-features = false, features = ["system"] }
 ndarray-npy = { version = "0.5", default-features = false }
 criterion = "0.3"
 serde_json = "1"
diff --git a/linfa-clustering/src/gaussian_mixture/algorithm.rs b/linfa-clustering/src/gaussian_mixture/algorithm.rs
index 6169f163f..e0d206150 100644
--- a/linfa-clustering/src/gaussian_mixture/algorithm.rs
+++ b/linfa-clustering/src/gaussian_mixture/algorithm.rs
@@ -54,7 +54,6 @@ use serde_crate::{Deserialize, Serialize};
 /// Let's do a walkthrough of a training-predict-save example.
 ///
 /// ```rust, ignore
-/// extern crate openblas_src;
 /// use linfa::Dataset;
 /// use linfa::traits::{Fit, Predict};
 /// use linfa_clustering::{GmmHyperParams, GaussianMixtureModel, generate_blobs};
@@ -478,8 +477,6 @@ impl<F: Float + Lapack + Scalar, D: Data<Elem = F>, T: Targets>
 
 #[cfg(test)]
 mod tests {
-    extern crate openblas_src;
-
     use super::*;
     use crate::generate_blobs;
     use approx::assert_abs_diff_eq;
diff --git a/linfa-clustering/src/k_means/algorithm.rs b/linfa-clustering/src/k_means/algorithm.rs
index 1e873c68f..1490c89aa 100644
--- a/linfa-clustering/src/k_means/algorithm.rs
+++ b/linfa-clustering/src/k_means/algorithm.rs
@@ -404,8 +404,6 @@ fn get_random_centroids<F: Float, D: Data<Elem = F>>(
 
 #[cfg(test)]
 mod tests {
-    extern crate openblas_src;
-
     use super::*;
     use approx::assert_abs_diff_eq;
     use ndarray::{array, stack, Array, Array1, Array2, Axis};
diff --git a/linfa-ica/Cargo.toml b/linfa-ica/Cargo.toml
index a227eafc2..dc62a6e5a 100644
--- a/linfa-ica/Cargo.toml
+++ b/linfa-ica/Cargo.toml
@@ -15,6 +15,7 @@ categories = ["algorithms", "mathematics", "science"]
 [features]
 default = []
 serde = ["serde_crate", "ndarray/serde"]                                                     
+
 [dependencies.serde_crate]
 package = "serde"
 optional = true
@@ -24,15 +25,14 @@ features = ["std", "derive"]
 
 [dependencies]
 ndarray = { version = "0.13", default-features = false }
+ndarray-linalg = "0.12"
 ndarray-rand = "0.11"
 ndarray-stats = "0.3"
-ndarray-linalg = "0.12"
 num-traits = "0.2"
 rand_isaac = "0.2.0"
 
 linfa = { version = "0.2.1", path = ".." }
 
 [dev-dependencies]
-openblas-src = { version = "0.9", default-features = false, features = ["system"] }
 ndarray-npy = { version = "0.5", default-features = false }
 paste = "1.0"
diff --git a/linfa-ica/examples/fast_ica.rs b/linfa-ica/examples/fast_ica.rs
index ce90d19fc..859d77c7f 100644
--- a/linfa-ica/examples/fast_ica.rs
+++ b/linfa-ica/examples/fast_ica.rs
@@ -1,5 +1,3 @@
-extern crate openblas_src;
-
 use linfa::{
     dataset::Dataset,
     traits::{Fit, Predict},
diff --git a/linfa-ica/src/fast_ica.rs b/linfa-ica/src/fast_ica.rs
index 25066df2e..df83e35e5 100644
--- a/linfa-ica/src/fast_ica.rs
+++ b/linfa-ica/src/fast_ica.rs
@@ -307,8 +307,6 @@ impl GFunc {
 
 #[cfg(test)]
 mod tests {
-    extern crate openblas_src;
-
     use super::*;
     use linfa::traits::{Fit, Predict};
 
diff --git a/linfa-linear/Cargo.toml b/linfa-linear/Cargo.toml
index ad3142f99..c0b0083fa 100644
--- a/linfa-linear/Cargo.toml
+++ b/linfa-linear/Cargo.toml
@@ -27,7 +27,6 @@ serde = { version = "1.0", default-features = false, features = ["derive"] }
 linfa = { version = "0.2.1", path = ".." }
 
 [dev-dependencies]
-openblas-src = { version = "0.9", default-features = false, features = ["system"] }
 csv = "1.1"
 ndarray-csv = "0.4"
 approx = "0.3.2"
diff --git a/linfa-linear/examples/diabetes.rs b/linfa-linear/examples/diabetes.rs
index c68acfb54..316336146 100644
--- a/linfa-linear/examples/diabetes.rs
+++ b/linfa-linear/examples/diabetes.rs
@@ -1,5 +1,3 @@
-extern crate openblas_src;
-
 use std::error::Error;
 use std::fs::File;
 
diff --git a/linfa-linear/examples/glm.rs b/linfa-linear/examples/glm.rs
index f55846e7d..97fd14fea 100644
--- a/linfa-linear/examples/glm.rs
+++ b/linfa-linear/examples/glm.rs
@@ -1,5 +1,3 @@
-extern crate openblas_src;
-
 use csv::ReaderBuilder;
 use flate2::read::GzDecoder;
 use linfa_linear::TweedieRegressor;
diff --git a/linfa-linear/src/ols.rs b/linfa-linear/src/ols.rs
index ae0f7d98a..c9878124e 100644
--- a/linfa-linear/src/ols.rs
+++ b/linfa-linear/src/ols.rs
@@ -255,8 +255,6 @@ impl<F: Float, D: Data<Elem = F>> Predict<&ArrayBase<D, Ix2>, Array1<F>>
 
 #[cfg(test)]
 mod tests {
-    extern crate openblas_src;
-
     use super::*;
     use approx::abs_diff_eq;
     use ndarray::array;
diff --git a/linfa-logistic/Cargo.toml b/linfa-logistic/Cargo.toml
index ba5e2feaf..e77cefa22 100644
--- a/linfa-logistic/Cargo.toml
+++ b/linfa-logistic/Cargo.toml
@@ -14,12 +14,13 @@ keywords = ["machine-learning", "linfa", "ai", "ml", "linear"]
 categories = ["algorithms", "mathematics", "science"]
 
 [dependencies]
-ndarray = {version = "0.13", features = ["blas", "approx"]}
+ndarray = {version = "0.13", features = ["approx", "blas"]}
 ndarray-linalg = "0.12"
 num-traits = "0.2"
 argmin = {version="0.3.1", features=["ndarrayl"]}
 serde = "1.0"
 
+linfa = { version = "0.2.1", path = ".." }
+
 [dev-dependencies]
-openblas-src = { version = "0.9", default-features = false, features = ["system"] }
 approx = "0.3.2"
diff --git a/linfa-logistic/src/lib.rs b/linfa-logistic/src/lib.rs
index 71c8a8cd1..f36b5ecc4 100644
--- a/linfa-logistic/src/lib.rs
+++ b/linfa-logistic/src/lib.rs
@@ -510,7 +510,7 @@ impl<'a, F: Float, A: Data<Elem = F>> ArgminOp for LogisticRegressionProblem<'a,
 
 #[cfg(test)]
 mod test {
-    extern crate openblas_src;
+    extern crate linfa;
 
     use super::*;
     use approx::AbsDiffEq;
diff --git a/linfa-reduction/Cargo.toml b/linfa-reduction/Cargo.toml
index f351df2ba..f4ddb0d0b 100644
--- a/linfa-reduction/Cargo.toml
+++ b/linfa-reduction/Cargo.toml
@@ -14,6 +14,7 @@ categories = ["algorithms", "mathematics", "science"]
 
 [features]
 default = []
+
 serde = ["serde_crate", "ndarray/serde"]
 
 [dependencies.serde_crate]
@@ -24,16 +25,15 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = "0.13"
+ndarray = { version = "0.13", default-features = false }
+ndarray-linalg = "0.12"
 ndarray-rand = "0.11"
 ndarray-stats = "0.3"
-ndarray-linalg = "0.12"
 num-traits = "0.2"
 
 linfa = { version = "0.2.1", path = ".." }
 linfa-kernel = { version = "0.2.1", path = "../linfa-kernel" }
 
 [dev-dependencies]
-openblas-src = { version = "0.9", default-features=false, features = ["system"] }
 rand_isaac = "0.2.0"
 ndarray-npy = { version = "0.5", default-features = false }
diff --git a/linfa-reduction/examples/diffusion_map.rs b/linfa-reduction/examples/diffusion_map.rs
index 5d0a5a32e..cd6fe1513 100644
--- a/linfa-reduction/examples/diffusion_map.rs
+++ b/linfa-reduction/examples/diffusion_map.rs
@@ -1,5 +1,3 @@
-extern crate openblas_src;
-
 use linfa::traits::Transformer;
 use linfa_kernel::{Kernel, KernelMethod, KernelType};
 use linfa_reduction::utils::generate_convoluted_rings2d;
diff --git a/linfa-reduction/examples/pca.rs b/linfa-reduction/examples/pca.rs
index 9b4d4350d..c03380c78 100644
--- a/linfa-reduction/examples/pca.rs
+++ b/linfa-reduction/examples/pca.rs
@@ -1,5 +1,3 @@
-extern crate openblas_src;
-
 use linfa::prelude::*;
 use linfa_reduction::{utils::generate_blobs, Pca};
 use ndarray::array;
diff --git a/linfa-reduction/src/diffusion_map/algorithms.rs b/linfa-reduction/src/diffusion_map/algorithms.rs
index 5427687ff..8ec5eaa27 100644
--- a/linfa-reduction/src/diffusion_map/algorithms.rs
+++ b/linfa-reduction/src/diffusion_map/algorithms.rs
@@ -1,11 +1,12 @@
 use ndarray::{Array1, Array2, ArrayView2};
 use ndarray_linalg::{
-    eigh::EighInto, lapack::UPLO, lobpcg, lobpcg::LobpcgResult, Lapack, Scalar, TruncatedOrder,
+    eigh::EighInto, lobpcg, lobpcg::LobpcgResult, Lapack, Scalar, TruncatedOrder, UPLO,
 };
 use ndarray_rand::{rand_distr::Uniform, RandomExt};
 use num_traits::NumCast;
 
 use linfa::{traits::Transformer, Float};
+
 use linfa_kernel::Kernel;
 
 use super::hyperparameters::{DiffusionMapHyperParams, DiffusionMapHyperParamsBuilder};
diff --git a/linfa-svm/Cargo.toml b/linfa-svm/Cargo.toml
index 1a3c85000..f66443647 100644
--- a/linfa-svm/Cargo.toml
+++ b/linfa-svm/Cargo.toml
@@ -25,6 +25,7 @@ features = ["std", "derive"]
 
 [dependencies]
 ndarray = { version = "0.13", default-features=false, features=["blas"] }
+ndarray-linalg = "0.12"
 ndarray-rand = "0.11"
 num-traits = "0.1.32"
 
@@ -32,7 +33,6 @@ linfa = { version = "0.2.1", path = ".." }
 linfa-kernel = { version = "0.2.1", path = "../linfa-kernel" }
 
 [dev-dependencies]
-openblas-src = { version = "0.9", default-features = false, features = ["system"] }
 csv = "1.1"
 ndarray-csv = "0.4"
 flate2 = "1.0"
diff --git a/linfa-svm/examples/winequality.rs b/linfa-svm/examples/winequality.rs
index 0aef42d42..452dd3062 100644
--- a/linfa-svm/examples/winequality.rs
+++ b/linfa-svm/examples/winequality.rs
@@ -1,5 +1,3 @@
-extern crate openblas_src;
-
 use std::error::Error;
 use std::fs::File;
 
diff --git a/linfa-svm/src/classification.rs b/linfa-svm/src/classification.rs
index 728872e00..6a11f2710 100644
--- a/linfa-svm/src/classification.rs
+++ b/linfa-svm/src/classification.rs
@@ -305,8 +305,6 @@ impl<'a, F: Float, T: Targets, D: Data<Elem = F>>
 
 #[cfg(test)]
 mod tests {
-    extern crate openblas_src;
-
     use super::Svm;
     use linfa::dataset::Dataset;
     use linfa::metrics::ToConfusionMatrix;
diff --git a/src/lib.rs b/src/lib.rs
index 00a74fc25..5b027194e 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -37,6 +37,18 @@ pub mod traits;
 
 pub use dataset::{Dataset, Float, Label};
 
+#[cfg(feature = "ndarray-linalg")]
+pub use ndarray_linalg as linalg;
+
+#[cfg(any(feature = "intel-mkl-system", feature = "intel-mkl-static"))]
+extern crate intel_mkl_src;
+
+#[cfg(any(feature = "openblas-system", feature = "openblas-static"))]
+extern crate openblas_src;
+
+#[cfg(any(feature = "netblas-system", feature = "netblas-static"))]
+extern crate netblas_src;
+
 /// Common metrics functions for classification and regression
 pub mod metrics {
     pub use crate::metrics_classification::{