diff --git a/src/NativeProviders/MKL/lapack.cpp b/src/NativeProviders/MKL/lapack.cpp
index 69b7e965..fb2601f1 100644
--- a/src/NativeProviders/MKL/lapack.cpp
+++ b/src/NativeProviders/MKL/lapack.cpp
@@ -1,5 +1,6 @@
 #include <algorithm>
 #include <complex>
+
 #define MKL_Complex8 std::complex<float>
 #define MKL_Complex16 std::complex<double>
 
@@ -11,21 +12,17 @@
 #include "mkl.h"
 #include "mkl_trans.h"
 
-template<typename T>
-inline MKL_INT lu_factor(MKL_INT m, T a[], MKL_INT ipiv[],
-                         void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*))
+template<typename T, typename GETRF>
+inline MKL_INT lu_factor(MKL_INT m, T a[], MKL_INT ipiv[], GETRF getrf)
 {
-    std::complex<double> x = 5;
 	MKL_INT info = 0;
     getrf(&m, &m, a, &m, ipiv, &info);
     shift_ipiv_down(m, ipiv);
     return info;
 };
 
-template<typename T>
-inline MKL_INT lu_inverse(MKL_INT n, T a[], T work[], MKL_INT lwork,
-                          void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*),
-                          void (*getri)(const MKL_INT*, T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GETRF, typename GETRI>
+inline MKL_INT lu_inverse(MKL_INT n, T a[], T work[], MKL_INT lwork, GETRF getrf, GETRI getri)
 {
     MKL_INT* ipiv = new MKL_INT[n];
     MKL_INT info = 0;
@@ -42,9 +39,8 @@ inline MKL_INT lu_inverse(MKL_INT n, T a[], T work[], MKL_INT lwork,
     return info;
 };
 
-template<typename T>
-inline MKL_INT lu_inverse_factored(MKL_INT n, T a[], MKL_INT ipiv[], T work[], MKL_INT lwork,
-                                   void (*getri)(const MKL_INT*, T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GETRI>
+inline MKL_INT lu_inverse_factored(MKL_INT n, T a[], MKL_INT ipiv[], T work[], MKL_INT lwork, GETRI getri)
 {
     shift_ipiv_up(n, ipiv);
     MKL_INT info = 0;
@@ -53,9 +49,8 @@ inline MKL_INT lu_inverse_factored(MKL_INT n, T a[], MKL_INT ipiv[], T work[], M
     return info;
 }
 
-template<typename T>
-inline MKL_INT lu_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], MKL_INT ipiv[], T b[],
-                                 void (*getrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GETRS>
+inline MKL_INT lu_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], MKL_INT ipiv[], T b[], GETRS getrs)
 {
     shift_ipiv_up(n, ipiv);
     MKL_INT info = 0;
@@ -65,10 +60,8 @@ inline MKL_INT lu_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], MKL_INT ipiv[],
     return info;
 }
 
-template<typename T>
-inline MKL_INT lu_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
-                        void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*),
-                        void (*getrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GETRF, typename GETRS>
+inline MKL_INT lu_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[], GETRF getrf, GETRS getrs)
 {
     T* clone = Clone(n, n, a);
     MKL_INT* ipiv = new MKL_INT[n];
@@ -90,9 +83,8 @@ inline MKL_INT lu_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
 }
 
 
-template<typename T>
-inline MKL_INT cholesky_factor(MKL_INT n, T* a,
-                               void (*potrf)(const char*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename POTRF>
+inline MKL_INT cholesky_factor(MKL_INT n, T* a, POTRF potrf)
 {
     char uplo = 'L';
     MKL_INT info = 0;
@@ -112,10 +104,8 @@ inline MKL_INT cholesky_factor(MKL_INT n, T* a,
     return info;
 }
 
-template<typename T>
-inline MKL_INT cholesky_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
-                              void (*potrf)(const char*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*),
-                              void (*potrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename POTRF, typename POTRS>
+inline MKL_INT cholesky_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[], POTRF potrf, POTRS potrs)
 {
     T* clone = Clone(n, n, a);
     char uplo = 'L';
@@ -133,9 +123,8 @@ inline MKL_INT cholesky_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
     return info;
 }
 
-template<typename T>
-inline MKL_INT cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], T b[],
-                                       void (*potrs)(const char*, const MKL_INT*, const MKL_INT*, const T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename POTRS>
+inline MKL_INT cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], T b[], POTRS potrs)
 {
     char uplo = 'L';
     MKL_INT info = 0;
@@ -143,10 +132,8 @@ inline MKL_INT cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, T a[], T b[],
     return info;
 }
 
-template<typename T>
-inline MKL_INT qr_factor(MKL_INT m, MKL_INT n, T r[], T tau[], T q[], T work[], MKL_INT len,
-                         void (*geqrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, T*, T*, const MKL_INT*, MKL_INT*),
-                         void (*orgqr)(const MKL_INT*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, const T*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GEQRF, typename ORGQR>
+inline MKL_INT qr_factor(MKL_INT m, MKL_INT n, T r[], T tau[], T q[], T work[], MKL_INT len, GEQRF geqrf, ORGQR orgqr)
 {
     MKL_INT info = 0;
     geqrf(&m, &n, r, &m, tau, work, &len, &info);
@@ -175,10 +162,8 @@ inline MKL_INT qr_factor(MKL_INT m, MKL_INT n, T r[], T tau[], T q[], T work[],
     return info;
 }
 
-template<typename T>
-inline MKL_INT qr_thin_factor(MKL_INT m, MKL_INT n, T q[], T tau[], T r[], T work[], MKL_INT len,
-                              void (*geqrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, T*, T*, const MKL_INT*, MKL_INT*),
-                              void (*orgqr)(const MKL_INT*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, const T*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GEQRF, typename ORGQR>
+inline MKL_INT qr_thin_factor(MKL_INT m, MKL_INT n, T q[], T tau[], T r[], T work[], MKL_INT len, GEQRF geqrf, ORGQR orgqr)
 {
     MKL_INT info = 0;
     geqrf(&m, &n, q, &m, tau, work, &len, &info);
@@ -198,10 +183,8 @@ inline MKL_INT qr_thin_factor(MKL_INT m, MKL_INT n, T q[], T tau[], T r[], T wor
     return info;
 }
 
-template<typename T>
-inline MKL_INT qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, T a[], T b[], T x[], T work[], MKL_INT len,
-                        void (*gels)(const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*, T*,
-                                     const MKL_INT*, T* b, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GELS>
+inline MKL_INT qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, T a[], T b[], T x[], T work[], MKL_INT len, GELS gels)
 {
     T* clone_a = Clone(m, n, a);
     T* clone_b = Clone(m, bn, b);
@@ -214,12 +197,8 @@ inline MKL_INT qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, T a[], T b[], T x[], T
     return info;
 }
 
-template<typename T>
-inline MKL_INT qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len,
-                                 void (*ormqr)(const char*, const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*,
-                                         const T*, const MKL_INT*, const T*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT* info),
-                                 void (*trsm)(const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG,
-                                         const MKL_INT, const MKL_INT, const T, const T*, const MKL_INT, T*, const MKL_INT))
+template<typename T, typename ORMQR, typename TRSM>
+inline MKL_INT qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len, ORMQR ormqr, TRSM trsm)
 {
     T* clone_b = Clone(m, bn, b);
     char side ='L';
@@ -232,12 +211,8 @@ inline MKL_INT qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[],
     return info;
 }
 
-template<typename T>
-inline MKL_INT complex_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len,
-        void (*unmqr)(const char*, const char*, const MKL_INT*, const MKL_INT*, const MKL_INT*,
-                      const T*, const MKL_INT*, const T*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT* info),
-        void (*trsm)(const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG,
-                     const MKL_INT, const MKL_INT, const void*, const void*, const MKL_INT, void*, const MKL_INT ldb))
+template<typename T, typename UNMQR, typename TRSM>
+inline MKL_INT complex_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[], T tau[], T x[], T work[], MKL_INT len, UNMQR unmqr, TRSM trsm)
 {
     T* clone_b = Clone(m, bn, b);
     char side ='L';
@@ -251,10 +226,8 @@ inline MKL_INT complex_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[]
     return info;
 }
 
-template<typename T>
-inline MKL_INT svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len,
-                          void (*gesvd)(const char*, const char*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*,
-                                        T*, T*, const MKL_INT*, T*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*))
+template<typename T, typename GESVD>
+inline MKL_INT svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len, GESVD gesvd)
 {
     MKL_INT info = 0;
     char job = compute_vectors ? 'A' : 'N';
@@ -262,10 +235,8 @@ inline MKL_INT svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s
     return info;
 }
 
-template<typename T, typename R>
-inline MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len,
-                                  void (*gesvd)(const char*, const char*, const MKL_INT*, const MKL_INT*, T*, const MKL_INT*,
-                                          R*, T*, const MKL_INT*, T*, const MKL_INT*, T*, const MKL_INT*, R*, MKL_INT*))
+template<typename T, typename R, typename GESVD>
+inline MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T a[], T s[], T u[], T v[], T work[], MKL_INT len, GESVD gesvd)
 {
     MKL_INT info = 0;
     MKL_INT dim_s = std::min(m,n);
@@ -284,12 +255,8 @@ inline MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T
     return info;
 }
 
-template<typename T>
-inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[],
-                    MKL_INT(*gees)(MKL_INT, char, char, int(*)(const T*, const T*), MKL_INT, T* a,
-                                   MKL_INT, MKL_INT*, T*, T*, T*, MKL_INT),
-                    MKL_INT(*trevc)(MKL_INT, char, char, lapack_logical*, MKL_INT, const T*,
-                                    MKL_INT, T*, MKL_INT, T*, MKL_INT, MKL_INT, MKL_INT*))
+template<typename T, typename R, typename GEES, typename TREVC>
+inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc)
 {
     T* clone_a = Clone(n, n, a);
     T* wr = new T[n];
@@ -317,7 +284,7 @@ inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[
 
     for (MKL_INT index = 0; index < n; ++index)
     {
-        values[index] = MKL_Complex16(wr[index], wi[index]);
+        values[index] = R(wr[index], wi[index]);
     }
 
     for (MKL_INT  i = 0; i < n; ++i)
@@ -341,12 +308,8 @@ inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[
     return info;
 }
 
-template<typename T>
-inline MKL_INT eigen_complex_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[],
-							MKL_INT(*gees)(MKL_INT, char, char, int(*)(const T*), MKL_INT, T* a,
-							MKL_INT, MKL_INT*, T*, T*, MKL_INT),
-							MKL_INT(*trevc)(MKL_INT, char, char, const lapack_logical*, MKL_INT, T*,
-							MKL_INT, T*, MKL_INT, T*, MKL_INT, MKL_INT, MKL_INT*))
+template<typename T, typename GEES, typename TREVC>
+inline MKL_INT eigen_complex_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], GEES gees, TREVC trevc)
 {
 	T* clone_a = Clone(n, n, a);
 	T* w = new T[n];
@@ -380,12 +343,11 @@ inline MKL_INT eigen_complex_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16
 	return info;
 }
 
-template<typename T, typename R>
-inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[],
-                    MKL_INT(*syev)(int, char, char, int, T*, int, R*))
+template<typename R, typename T, typename SYEV>
+inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], SYEV syev)
 {
     T* clone_a = Clone(n, n, a);
-    R* w = new R[n];
+	R* w = new R[n];
 
 	MKL_INT info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w);
 	if (info != 0)
@@ -394,7 +356,7 @@ inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 val
 		delete[] w;
 		return info;
 	}
-
+	
 	memcpy(vectors, clone_a, n*n*sizeof(T));
 
 	for (MKL_INT index = 0; index < n; ++index)
@@ -402,7 +364,7 @@ inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 val
 		values[index] = MKL_Complex16(w[index]);
 	}
     
-	for (MKL_INT j = 0; j < n; j++)
+	for (MKL_INT j = 0; j < n; ++j)
     {
         MKL_INT jn = j*n;
 
@@ -444,252 +406,252 @@ extern "C" {
 
     DLLEXPORT MKL_INT s_lu_factor(MKL_INT m, float a[], MKL_INT ipiv[])
     {
-        return lu_factor<float>(m, a, ipiv, sgetrf);
+        return lu_factor(m, a, ipiv, sgetrf);
     }
 
     DLLEXPORT MKL_INT d_lu_factor(MKL_INT m, double a[], MKL_INT ipiv[])
     {
-        return lu_factor<double>(m, a, ipiv, dgetrf);
+        return lu_factor(m, a, ipiv, dgetrf);
     }
 
     DLLEXPORT MKL_INT c_lu_factor(MKL_INT m, MKL_Complex8 a[], MKL_INT ipiv[])
     {
-        return lu_factor<MKL_Complex8>(m, a, ipiv, cgetrf);
+        return lu_factor(m, a, ipiv, cgetrf);
     }
 
     DLLEXPORT MKL_INT z_lu_factor(MKL_INT m, MKL_Complex16 a[], MKL_INT ipiv[])
     {
-        return lu_factor<MKL_Complex16>(m, a, ipiv, zgetrf);
+        return lu_factor(m, a, ipiv, zgetrf);
     }
 
     DLLEXPORT MKL_INT s_lu_inverse(MKL_INT n, float a[], float work[], MKL_INT lwork)
     {
-        return lu_inverse<float>(n, a, work, lwork, sgetrf, sgetri);
+        return lu_inverse(n, a, work, lwork, sgetrf, sgetri);
     }
 
     DLLEXPORT MKL_INT d_lu_inverse(MKL_INT n, double a[], double work[], MKL_INT lwork)
     {
-        return lu_inverse<double>(n, a, work, lwork, dgetrf, dgetri);
+        return lu_inverse(n, a, work, lwork, dgetrf, dgetri);
     }
 
     DLLEXPORT MKL_INT c_lu_inverse(MKL_INT n, MKL_Complex8 a[], MKL_Complex8 work[], MKL_INT lwork)
     {
-        return lu_inverse<MKL_Complex8>(n, a, work, lwork, cgetrf, cgetri);
+        return lu_inverse(n, a, work, lwork, cgetrf, cgetri);
     }
 
     DLLEXPORT MKL_INT z_lu_inverse(MKL_INT n, MKL_Complex16 a[], MKL_Complex16 work[], MKL_INT lwork)
     {
-        return lu_inverse<MKL_Complex16>(n, a, work, lwork, zgetrf, zgetri);
+        return lu_inverse(n, a, work, lwork, zgetrf, zgetri);
     }
 
     DLLEXPORT MKL_INT s_lu_inverse_factored(MKL_INT n, float a[], MKL_INT ipiv[], float work[], MKL_INT lwork)
     {
-        return lu_inverse_factored<float>(n, a, ipiv, work, lwork, sgetri);
+        return lu_inverse_factored(n, a, ipiv, work, lwork, sgetri);
     }
 
     DLLEXPORT MKL_INT d_lu_inverse_factored(MKL_INT n, double a[], MKL_INT ipiv[], double work[], MKL_INT lwork)
     {
-        return lu_inverse_factored<double>(n, a, ipiv, work, lwork, dgetri);
+        return lu_inverse_factored(n, a, ipiv, work, lwork, dgetri);
     }
 
     DLLEXPORT MKL_INT c_lu_inverse_factored(MKL_INT n, MKL_Complex8 a[], MKL_INT ipiv[], MKL_Complex8 work[], MKL_INT lwork)
     {
-        return lu_inverse_factored<MKL_Complex8>(n, a, ipiv, work, lwork, cgetri);
+        return lu_inverse_factored(n, a, ipiv, work, lwork, cgetri);
     }
 
     DLLEXPORT MKL_INT z_lu_inverse_factored(MKL_INT n, MKL_Complex16 a[], MKL_INT ipiv[], MKL_Complex16 work[], MKL_INT lwork)
     {
-        return lu_inverse_factored<MKL_Complex16>(n, a, ipiv, work, lwork, zgetri);
+        return lu_inverse_factored(n, a, ipiv, work, lwork, zgetri);
     }
 
     DLLEXPORT MKL_INT s_lu_solve_factored(MKL_INT n, MKL_INT nrhs, float a[], MKL_INT ipiv[], float b[])
     {
-        return lu_solve_factored<float>(n, nrhs, a, ipiv, b, sgetrs);
+        return lu_solve_factored(n, nrhs, a, ipiv, b, sgetrs);
     }
 
     DLLEXPORT MKL_INT  d_lu_solve_factored(MKL_INT n, MKL_INT nrhs, double a[], MKL_INT ipiv[], double b[])
     {
-        return lu_solve_factored<double>(n, nrhs, a, ipiv, b, dgetrs);
+        return lu_solve_factored(n, nrhs, a, ipiv, b, dgetrs);
     }
 
     DLLEXPORT MKL_INT c_lu_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_INT ipiv[], MKL_Complex8 b[])
     {
-        return lu_solve_factored<MKL_Complex8>(n, nrhs, a, ipiv, b, cgetrs);
+        return lu_solve_factored(n, nrhs, a, ipiv, b, cgetrs);
     }
 
     DLLEXPORT MKL_INT z_lu_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[], MKL_INT ipiv[], MKL_Complex16 b[])
     {
-        return lu_solve_factored<MKL_Complex16>(n, nrhs, a, ipiv, b, zgetrs);
+        return lu_solve_factored(n, nrhs, a, ipiv, b, zgetrs);
     }
 
     DLLEXPORT MKL_INT s_lu_solve(MKL_INT n, MKL_INT nrhs, float a[], float b[])
     {
-        return lu_solve<float>(n, nrhs, a, b, sgetrf, sgetrs);
+        return lu_solve(n, nrhs, a, b, sgetrf, sgetrs);
     }
 
     DLLEXPORT MKL_INT d_lu_solve(MKL_INT n, MKL_INT nrhs, double a[], double b[])
     {
-        return lu_solve<double>(n, nrhs, a, b, dgetrf, dgetrs);
+        return lu_solve(n, nrhs, a, b, dgetrf, dgetrs);
     }
 
     DLLEXPORT MKL_INT c_lu_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_Complex8 b[])
     {
-        return lu_solve<MKL_Complex8>(n, nrhs, a, b, cgetrf, cgetrs);
+        return lu_solve(n, nrhs, a, b, cgetrf, cgetrs);
     }
 
     DLLEXPORT MKL_INT z_lu_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[],  MKL_Complex16 b[])
     {
-        return lu_solve<MKL_Complex16>(n, nrhs, a, b, zgetrf, zgetrs);
+        return lu_solve(n, nrhs, a, b, zgetrf, zgetrs);
     }
 
     DLLEXPORT MKL_INT s_cholesky_factor(MKL_INT n, float a[])
     {
-        return cholesky_factor<float>(n, a, spotrf);
+        return cholesky_factor(n, a, spotrf);
     }
 
     DLLEXPORT MKL_INT d_cholesky_factor(MKL_INT n, double* a)
     {
-        return cholesky_factor<double>(n, a, dpotrf);
+        return cholesky_factor(n, a, dpotrf);
     }
 
     DLLEXPORT MKL_INT c_cholesky_factor(MKL_INT n, MKL_Complex8 a[])
     {
-        return cholesky_factor<MKL_Complex8>(n, a, cpotrf);
+        return cholesky_factor(n, a, cpotrf);
     }
 
     DLLEXPORT MKL_INT z_cholesky_factor(MKL_INT n, MKL_Complex16 a[])
     {
-        return cholesky_factor<MKL_Complex16>(n, a, zpotrf);
+        return cholesky_factor(n, a, zpotrf);
     }
 
     DLLEXPORT MKL_INT s_cholesky_solve(MKL_INT n, MKL_INT nrhs, float a[], float b[])
     {
-        return cholesky_solve<float>(n, nrhs, a, b, spotrf, spotrs);
+        return cholesky_solve(n, nrhs, a, b, spotrf, spotrs);
     }
 
     DLLEXPORT MKL_INT d_cholesky_solve(MKL_INT n, MKL_INT nrhs, double a[], double b[])
     {
-        return cholesky_solve<double>(n, nrhs, a, b, dpotrf, dpotrs);
+        return cholesky_solve(n, nrhs, a, b, dpotrf, dpotrs);
     }
 
     DLLEXPORT MKL_INT c_cholesky_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_Complex8 b[])
     {
-        return cholesky_solve<MKL_Complex8>(n, nrhs, a, b, cpotrf, cpotrs);
+        return cholesky_solve(n, nrhs, a, b, cpotrf, cpotrs);
     }
 
     DLLEXPORT MKL_INT z_cholesky_solve(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[], MKL_Complex16 b[])
     {
-        return cholesky_solve<MKL_Complex16>(n, nrhs, a, b, zpotrf, zpotrs);
+        return cholesky_solve(n, nrhs, a, b, zpotrf, zpotrs);
     }
 
     DLLEXPORT MKL_INT s_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, float a[], float b[])
     {
-        return cholesky_solve_factored<float>(n, nrhs, a, b, spotrs);
+        return cholesky_solve_factored(n, nrhs, a, b, spotrs);
     }
 
     DLLEXPORT MKL_INT d_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, double a[], double b[])
     {
-        return cholesky_solve_factored<double>(n, nrhs, a, b, dpotrs);
+        return cholesky_solve_factored(n, nrhs, a, b, dpotrs);
     }
 
     DLLEXPORT MKL_INT c_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex8 a[], MKL_Complex8 b[])
     {
-        return cholesky_solve_factored<MKL_Complex8>(n, nrhs, a, b, cpotrs);
+        return cholesky_solve_factored(n, nrhs, a, b, cpotrs);
     }
 
     DLLEXPORT MKL_INT z_cholesky_solve_factored(MKL_INT n, MKL_INT nrhs, MKL_Complex16 a[], MKL_Complex16 b[])
     {
-        return cholesky_solve_factored<MKL_Complex16>(n, nrhs, a, b, zpotrs);
+        return cholesky_solve_factored(n, nrhs, a, b, zpotrs);
     }
 
     DLLEXPORT MKL_INT s_qr_factor(MKL_INT m, MKL_INT n, float r[], float tau[], float q[], float work[], MKL_INT len)
     {
-        return qr_factor<float>(m, n, r, tau, q, work, len, sgeqrf, sorgqr);
+        return qr_factor(m, n, r, tau, q, work, len, sgeqrf, sorgqr);
     }
 
     DLLEXPORT MKL_INT s_qr_thin_factor(MKL_INT m, MKL_INT n, float q[], float tau[], float r[], float work[], MKL_INT len)
     {
-        return qr_thin_factor<float>(m, n, q, tau, r, work, len, sgeqrf, sorgqr);
+        return qr_thin_factor(m, n, q, tau, r, work, len, sgeqrf, sorgqr);
     }
 
     DLLEXPORT MKL_INT d_qr_factor(MKL_INT m, MKL_INT n, double r[], double tau[], double q[], double work[], MKL_INT len)
     {
-        return qr_factor<double>(m, n, r, tau, q, work, len, dgeqrf, dorgqr);
+        return qr_factor(m, n, r, tau, q, work, len, dgeqrf, dorgqr);
     }
 
     DLLEXPORT MKL_INT d_qr_thin_factor(MKL_INT m, MKL_INT n, double q[], double tau[], double r[], double work[], MKL_INT len)
     {
-        return qr_thin_factor<double>(m, n, q, tau, r, work, len, dgeqrf, dorgqr);
+        return qr_thin_factor(m, n, q, tau, r, work, len, dgeqrf, dorgqr);
     }
 
     DLLEXPORT MKL_INT c_qr_factor(MKL_INT m, MKL_INT n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[], MKL_Complex8 work[], MKL_INT len)
     {
-        return qr_factor<MKL_Complex8>(m, n, r, tau, q, work, len, cgeqrf, cungqr);
+        return qr_factor(m, n, r, tau, q, work, len, cgeqrf, cungqr);
     }
 
     DLLEXPORT MKL_INT c_qr_thin_factor(MKL_INT m, MKL_INT n, MKL_Complex8 q[], MKL_Complex8 tau[], MKL_Complex8 r[], MKL_Complex8 work[], MKL_INT len)
     {
-        return qr_thin_factor<MKL_Complex8>(m, n, q, tau, r, work, len, cgeqrf, cungqr);
+        return qr_thin_factor(m, n, q, tau, r, work, len, cgeqrf, cungqr);
     }
 
     DLLEXPORT MKL_INT z_qr_factor(MKL_INT m, MKL_INT n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[], MKL_Complex16 work[], MKL_INT len)
     {
-        return qr_factor<MKL_Complex16>(m, n, r, tau, q, work, len, zgeqrf, zungqr);
+        return qr_factor(m, n, r, tau, q, work, len, zgeqrf, zungqr);
     }
 
     DLLEXPORT MKL_INT z_qr_thin_factor(MKL_INT m, MKL_INT n, MKL_Complex16 q[], MKL_Complex16 tau[], MKL_Complex16 r[], MKL_Complex16 work[], MKL_INT len)
     {
-        return qr_thin_factor<MKL_Complex16>(m, n, q, tau, r, work, len, zgeqrf, zungqr);
+        return qr_thin_factor(m, n, q, tau, r, work, len, zgeqrf, zungqr);
     }
 
     DLLEXPORT MKL_INT s_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, float a[], float b[], float x[], float work[], MKL_INT len)
     {
-        return qr_solve<float>(m, n, bn, a, b, x, work, len, sgels);
+        return qr_solve(m, n, bn, a, b, x, work, len, sgels);
     }
 
     DLLEXPORT MKL_INT d_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, double a[], double b[], double x[], double work[], MKL_INT len)
     {
-        return qr_solve<double>(m, n, bn, a, b, x, work, len, dgels);
+        return qr_solve(m, n, bn, a, b, x, work, len, dgels);
     }
 
     DLLEXPORT MKL_INT c_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex8 a[], MKL_Complex8 b[], MKL_Complex8 x[], MKL_Complex8 work[], MKL_INT len)
     {
-        return qr_solve<MKL_Complex8>(m, n, bn, a, b, x, work, len, cgels);
+        return qr_solve(m, n, bn, a, b, x, work, len, cgels);
     }
 
     DLLEXPORT MKL_INT z_qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex16 a[], MKL_Complex16 b[], MKL_Complex16 x[], MKL_Complex16 work[], MKL_INT len)
     {
-        return qr_solve<MKL_Complex16>(m, n, bn, a, b, x, work, len, zgels);
+        return qr_solve(m, n, bn, a, b, x, work, len, zgels);
     }
 
     DLLEXPORT MKL_INT s_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, float r[], float b[], float tau[], float x[], float work[], MKL_INT len)
     {
-        return qr_solve_factored<float>(m, n, bn, r, b, tau, x, work, len, sormqr, cblas_strsm);
+        return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, sormqr, cblas_strsm);
     }
 
     DLLEXPORT MKL_INT d_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, double r[], double b[], double tau[], double x[], double work[], MKL_INT len)
     {
-        return qr_solve_factored<double>(m, n, bn, r, b, tau, x, work, len, dormqr, cblas_dtrsm);
+        return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, dormqr, cblas_dtrsm);
     }
 
     DLLEXPORT MKL_INT c_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[], MKL_Complex8 work[], MKL_INT len)
     {
-        return complex_qr_solve_factored<MKL_Complex8>(m, n, bn, r, b, tau, x, work, len, cunmqr, cblas_ctrsm);
+        return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, cunmqr, cblas_ctrsm);
     }
 
     DLLEXPORT MKL_INT z_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[], MKL_Complex16 work[], MKL_INT len)
     {
-        return complex_qr_solve_factored<MKL_Complex16>(m, n, bn, r, b, tau, x, work, len, zunmqr, cblas_ztrsm);
+        return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, zunmqr, cblas_ztrsm);
     }
 
     DLLEXPORT MKL_INT s_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, float a[], float s[], float u[], float v[], float work[], MKL_INT len)
     {
-        return svd_factor<float>(compute_vectors, m, n, a, s, u, v, work, len, sgesvd);
+        return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, sgesvd);
     }
 
     DLLEXPORT MKL_INT d_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, double a[], double s[], double u[], double v[], double work[], MKL_INT len)
     {
-        return svd_factor<double>(compute_vectors, m, n, a, s, u, v, work, len, dgesvd);
+        return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, dgesvd);
     }
 
     DLLEXPORT MKL_INT c_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, MKL_Complex8 a[], MKL_Complex8 s[], MKL_Complex8 u[], MKL_Complex8 v[], MKL_Complex8 work[], MKL_INT len)
@@ -706,23 +668,23 @@ extern "C" {
     {
 		if (isSymmetric)
 		{
-			return sym_eigen_factor<float, float>(n, a, vectors, values, d, LAPACKE_ssyev);
+			return sym_eigen_factor<float>(n, a, vectors, values, d, LAPACKE_ssyev);
 		}
 		else
 		{
-			return eigen_factor<float>(n, a, vectors, values, d, LAPACKE_sgees, LAPACKE_strevc);
+			return eigen_factor(n, a, vectors, values, d, LAPACKE_sgees, LAPACKE_strevc);
 		}
     }
 
-    DLLEXPORT MKL_INT d_eigen(bool isSymmetric, MKL_INT n, double a[], double vectors[], MKL_Complex16 values[], double d[])
+   DLLEXPORT MKL_INT d_eigen(bool isSymmetric, MKL_INT n, double a[], double vectors[], MKL_Complex16 values[], double d[])
     {
         if (isSymmetric)
         {
-            return sym_eigen_factor<double, double>(n, a, vectors, values, d, LAPACKE_dsyev);
+			return sym_eigen_factor<double>(n, a, vectors, values, d, LAPACKE_dsyev);
         }
         else
         {
-            return eigen_factor<double>(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc);
+			return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc);
         }
     }
 
@@ -730,24 +692,23 @@ extern "C" {
     {
 		if (isSymmetric)
 		{
-			return sym_eigen_factor<MKL_Complex8, float>(n, a, vectors, values, d, LAPACKE_cheev);
+			return sym_eigen_factor<float>(n, a, vectors, values, d, LAPACKE_cheev);
 		}
 		else
 		{
-			return -1;
-			//return eigen_factor<MKL_Complex16, LAPACK_Z_SELECT1>(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc);
+			return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc);
 		}
     }
-
+	
     DLLEXPORT MKL_INT z_eigen(bool isSymmetric, MKL_INT n, MKL_Complex16 a[], MKL_Complex16 vectors[], MKL_Complex16 values[], MKL_Complex16 d[])
     {
 		if (isSymmetric)
 		{
-			return sym_eigen_factor<MKL_Complex16, double>(n, a, vectors, values, d, LAPACKE_zheev);
+			return sym_eigen_factor<double>(n, a, vectors, values, d, LAPACKE_zheev);
 		}
 		else
 		{
-			return eigen_complex_factor<MKL_Complex16>(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc);
+			return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc);
 		}
     }
 }
diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs
index 7a4b8e71..442e84c6 100644
--- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs
+++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs
@@ -1172,7 +1172,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
                 throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
             }
 
-            SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length);
+            if (SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0)
+            {
+                throw new NonConvergenceException();
+            }
         }
 
         /// <summary>
@@ -1314,6 +1317,63 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
 
             SafeNativeMethods.z_vector_divide(x.Length, x, y, result);
         }
+
+        /// <summary>
+        /// Computes the eigenvalues and eigenvectors of a matrix.
+        /// </summary>
+        /// <param name="isSymmetric">Wether the matrix is symmetric or not.</param>
+        /// <param name="order">The order of the matrix.</param>
+        /// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
+        /// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
+        /// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
+        /// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
+        public override void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD)
+        {
+            if (matrix == null)
+            {
+                throw new ArgumentNullException("matrix");
+            }
+
+            if (matrix.Length != order * order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix");
+            }
+
+            if (matrixEv == null)
+            {
+                throw new ArgumentNullException("matrixEv");
+            }
+
+            if (matrixEv.Length != order * order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv");
+            }
+
+            if (vectorEv == null)
+            {
+                throw new ArgumentNullException("vectorEv");
+            }
+
+            if (vectorEv.Length != order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv");
+            }
+
+            if (matrixD == null)
+            {
+                throw new ArgumentNullException("matrixD");
+            }
+
+            if (matrixD.Length != order * order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD");
+            }
+
+            if (SafeNativeMethods.z_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0)
+            {
+                throw new NonConvergenceException();
+            }
+        }
     }
 }
 
diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs
index 01322ae4..c903521c 100644
--- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs
+++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs
@@ -31,6 +31,7 @@
 #if NATIVEMKL
 
 using System;
+using System.Numerics;
 using System.Security;
 using MathNet.Numerics.LinearAlgebra.Factorization;
 using MathNet.Numerics.Properties;
@@ -1171,7 +1172,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
                 throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
             }
 
-            SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length);
+            if (SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0)
+            {
+                throw new NonConvergenceException();
+            }
         }
 
         /// <summary>
@@ -1313,6 +1317,63 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
 
             SafeNativeMethods.c_vector_divide(x.Length, x, y, result);
         }
+
+        /// <summary>
+        /// Computes the eigenvalues and eigenvectors of a matrix.
+        /// </summary>
+        /// <param name="isSymmetric">Wether the matrix is symmetric or not.</param>
+        /// <param name="order">The order of the matrix.</param>
+        /// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
+        /// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
+        /// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
+        /// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
+        public override void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD)
+        {
+            if (matrix == null)
+            {
+                throw new ArgumentNullException("matrix");
+            }
+
+            if (matrix.Length != order * order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix");
+            }
+
+            if (matrixEv == null)
+            {
+                throw new ArgumentNullException("matrixEv");
+            }
+
+            if (matrixEv.Length != order * order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv");
+            }
+
+            if (vectorEv == null)
+            {
+                throw new ArgumentNullException("vectorEv");
+            }
+
+            if (vectorEv.Length != order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv");
+            }
+
+            if (matrixD == null)
+            {
+                throw new ArgumentNullException("matrixD");
+            }
+
+            if (matrixD.Length != order * order)
+            {
+                throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD");
+            }
+
+            if (SafeNativeMethods.c_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0)
+            {
+                throw new NonConvergenceException();
+            }
+        }
     }
 }
 
diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs
index 7ba45140..b78b9995 100644
--- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs
+++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs
@@ -1172,7 +1172,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
                 throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
             }
 
-            SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length);
+            if (SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0)
+            {
+                throw new NonConvergenceException();
+            }
         }
 
         /// <summary>
@@ -1366,7 +1369,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
                 throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrixD");
             }
 
-            SafeNativeMethods.d_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD);
+            if (SafeNativeMethods.d_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0)
+            {
+                throw new NonConvergenceException();
+            }
         }
     }
 }
diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs
index 0297d638..3cd86986 100644
--- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs
+++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs
@@ -1172,7 +1172,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
                 throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
             }
 
-            SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length);
+            if (SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0)
+            {
+                throw new NonConvergenceException();
+            }
         }
 
         /// <summary>
@@ -1366,7 +1369,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
                 throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order*order), "matrixD");
             }
 
-            SafeNativeMethods.s_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD);
+            if (SafeNativeMethods.s_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0)
+            {
+                throw new NonConvergenceException();
+            }
         }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs
index 079e1fbd..e120c5f1 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs
@@ -28,6 +28,7 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System;
 using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.LinearAlgebra.Complex;
 using NUnit.Framework;
@@ -91,7 +92,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
         {
             var A = Matrix<Complex>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
-            var factorEvd = A.Evd();
+            var factorEvd = A.Evd(Symmetricity.Hermitian);
             var V = factorEvd.EigenVectors;
             var λ = factorEvd.D;
 
@@ -151,7 +152,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = Matrix<Complex>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = Vector<Complex>.Build.Random(order, 2);
             var bCopy = b.Clone();
@@ -179,7 +180,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = Matrix<Complex>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = Matrix<Complex>.Build.Random(order, order, 2);
             var BCopy = B.Clone();
@@ -212,7 +213,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = Matrix<Complex>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = Vector<Complex>.Build.Random(order, 2);
             var bCopy = b.Clone();
@@ -240,7 +241,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = Matrix<Complex>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = Matrix<Complex>.Build.Random(order, order, 2);
             var BCopy = B.Clone();
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs
index f2f01f8a..f995c472 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs
@@ -117,7 +117,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
         public void CanFactorizeRandomSymmetricMatrix(int order)
         {
             var matrixA = new UserDefinedMatrix(Matrix<Complex>.Build.RandomPositiveDefinite(order, 1).ToArray());
-            var factorEvd = matrixA.Evd();
+            var factorEvd = matrixA.Evd(Symmetricity.Hermitian);
             var eigenVectors = factorEvd.EigenVectors;
             var d = factorEvd.D;
 
@@ -204,7 +204,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = new UserDefinedVector(Vector<Complex>.Build.Random(order, 1).ToArray());
             var bCopy = b.Clone();
@@ -231,7 +231,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = new UserDefinedMatrix(Matrix<Complex>.Build.Random(order, order, 1).ToArray());
             var BCopy = B.Clone();
@@ -264,7 +264,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = new UserDefinedVector(Vector<Complex>.Build.Random(order, 1).ToArray());
             var bCopy = b.Clone();
@@ -292,7 +292,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = new UserDefinedMatrix(Matrix<Complex>.Build.Random(order, order, 1).ToArray());
             var BCopy = B.Clone();
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs
index cf7e08b6..37be44e4 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs
@@ -92,7 +92,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
         {
             var A = Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
-            var factorEvd = A.Evd();
+            var factorEvd = A.Evd(Symmetricity.Hermitian);
             var V = factorEvd.EigenVectors;
             var λ = factorEvd.D;
 
@@ -152,7 +152,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = Vector<Complex32>.Build.Random(order, 2);
             var bCopy = b.Clone();
@@ -179,7 +179,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = Matrix<Complex32>.Build.Random(order, order, 2);
             var BCopy = B.Clone();
@@ -212,7 +212,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = Vector<Complex32>.Build.Random(order, 2);
             var bCopy = b.Clone();
@@ -240,7 +240,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1);
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = Matrix<Complex32>.Build.Random(order, order, 2);
             var BCopy = B.Clone();
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs
index 19f17c91..1c32b514 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs
@@ -115,7 +115,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
         public void CanFactorizeRandomSymmetricMatrix([Values(1, 2, 5, 10, 50, 100)] int order)
         {
             var matrixA = new UserDefinedMatrix(Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1).ToArray());
-            var factorEvd = matrixA.Evd();
+            var factorEvd = matrixA.Evd(Symmetricity.Hermitian);
             var eigenVectors = factorEvd.EigenVectors;
             var d = factorEvd.D;
 
@@ -203,7 +203,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = new UserDefinedVector(Vector<Complex32>.Build.Random(order, 1).ToArray());
             var bCopy = b.Clone();
@@ -230,7 +230,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = new UserDefinedMatrix(Matrix<Complex32>.Build.Random(order, order, 1).ToArray());
             var BCopy = B.Clone();
@@ -263,7 +263,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var b = new UserDefinedVector(Vector<Complex32>.Build.Random(order, 1).ToArray());
             var bCopy = b.Clone();
@@ -291,7 +291,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
             var A = new UserDefinedMatrix(Matrix<Complex32>.Build.RandomPositiveDefinite(order, 1).ToArray());
             MatrixHelpers.ForceHermitian(A);
             var ACopy = A.Clone();
-            var evd = A.Evd();
+            var evd = A.Evd(Symmetricity.Hermitian);
 
             var B = new UserDefinedMatrix(Matrix<Complex32>.Build.Random(order, order, 1).ToArray());
             var BCopy = B.Clone();