Native pull: native Evd, updated provider

Squashed commit of the following:

commit 03291ba72c61e66e02c7a315e401210f96f9f20c
Author: Marcus Cuda <marcus@cuda.org>
Date:   Thu Apr 18 11:18:03 2013 +0300

    Tweaked managed EVD to be native compatible (use 1D arrays) and added naive EvD

commit d7bd0e7c93ba0284a7e8bee22f7d0c56f451febc
Author: Marcus Cuda <marcus@cuda.org>
Date:   Wed Feb 6 12:46:53 2013 +0200

    beginnings go an ATLAS provider

commit c4470fb99d5217a75526fc99e12c9119c12c3ec7
Author: Marcus Cuda <marcus@cuda.org>
Date:   Tue Feb 5 08:13:25 2013 -0800

    tested 32bit version

commit ea58e3bd936a6886130954e7ade63f066d6197cb
Author: Marcus Cuda <marcus@cuda.org>
Date:   Tue Feb 5 05:26:20 2013 -0800

    tweaked code and file name to compile with GCC and to run with mono on linux

src/NativeWrappers/ATLAS/blas.c

@@ -0,0 +1,101 @@
+#include "wrapper_common.h"
+#include "blas.h"
+
+#if GCC 
+extern "C" { 
+#endif
+DLLEXPORT void s_axpy(const int n, const float alpha, const float x[], float y[]){
+	cblas_saxpy(n, alpha, x, 1, y, 1);
+}
+
+DLLEXPORT void d_axpy(const int n, const double alpha, const double x[], double y[]){
+	cblas_daxpy(n, alpha, x, 1, y, 1);
+}
+
+DLLEXPORT void c_axpy(const int n, const Complex8 alpha, const Complex8 x[], Complex8 y[]){
+	cblas_caxpy(n, &alpha, x, 1, y, 1);
+}
+
+DLLEXPORT void z_axpy(const int n, const Complex16 alpha, const Complex16 x[], Complex16 y[]){
+	cblas_zaxpy(n, &alpha, x, 1, y, 1);
+}
+
+DLLEXPORT void s_scale(const int n, const float alpha, float x[]){
+	cblas_sscal(n, alpha, x, 1);
+}
+
+DLLEXPORT void d_scale(const int n, const double alpha, double x[]){
+	cblas_dscal(n, alpha, x, 1);
+}
+
+DLLEXPORT void c_scale(const int n, const Complex8 alpha, Complex8 x[]){
+	cblas_cscal(n, &alpha, x, 1);
+}
+
+DLLEXPORT void z_scale(const int n, const Complex16 alpha, Complex16 x[]){
+	cblas_zscal(n, &alpha, x, 1);
+}
+
+DLLEXPORT float s_dot_product(const int n, const float x[], const float y[]){
+	return cblas_sdot(n, x, 1, y, 1);
+}
+
+DLLEXPORT double d_dot_product(const int n, const double x[], const double y[]){
+	return cblas_ddot(n, x, 1, y, 1);
+}
+
+DLLEXPORT Complex8 c_dot_product(const int n, const Complex8 x[], const Complex8 y[]){
+	Complex8 ret;
+	cblas_cdotu_sub(n, x, 1, y, 1, &ret);
+	return ret;
+}
+
+DLLEXPORT Complex16 z_dot_product(const int n, const Complex16 x[], const Complex16 y[]){
+	Complex16 ret;
+	cblas_zdotu_sub(n, x, 1, y, 1, &ret);
+	return ret;
+}
+
+DLLEXPORT void s_matrix_multiply(const enum CBLAS_TRANSPOSE transA, const enum CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const float alpha, const float x[], const float y[], const float beta, float c[]){
+	int lda = transA == CblasNoTrans ? m : k;
+	int ldb = transB == CblasNoTrans ? k : n;
+
+	cblas_sgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m);
+}
+
+DLLEXPORT void d_matrix_multiply(const enum CBLAS_TRANSPOSE transA, const enum CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const double alpha, const double x[], const double y[], const double beta, double c[]){
+	int lda = transA == CblasNoTrans ? m : k;
+	int ldb = transB == CblasNoTrans ? k : n;
+
+	cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m);
+}
+
+DLLEXPORT void c_matrix_multiply(const enum CBLAS_TRANSPOSE transA, const enum CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const Complex8 alpha, const Complex8 x[], const Complex8 y[], const Complex8 beta, Complex8 c[]){
+	int lda = transA == CblasNoTrans ? m : k;
+	int ldb = transB == CblasNoTrans ? k : n;
+
+	cblas_cgemm(CblasColMajor, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m);
+}
+
+DLLEXPORT void z_matrix_multiply(const enum CBLAS_TRANSPOSE transA, const enum CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const Complex16 alpha, const Complex16 x[], const Complex16 y[], const Complex16 beta, Complex16 c[]){
+	int lda = transA == CblasNoTrans ? m : k;
+	int ldb = transB == CblasNoTrans ? k : n;
+
+	cblas_zgemm(CblasColMajor, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m);
+}
+
+/*char getTransChar(enum TRANSPOSE trans){
+	char cTrans;
+	switch( trans ){
+		case  CblasNoTrans : cTrans = 'N';
+			break;
+		case  CblasTrans : cTrans = 'T';
+			break;
+		case  CblasConjTrans : cTrans = 'C';
+			break;
+	}
+	return cTrans;
+}*/
+#if GCC 
+}
+#endif

src/NativeWrappers/ATLAS/lapack.cpp

@@ -1,64 +1,525 @@
-#include "common.h"
+#include "lapack_common.h"
+#include "wrapper_common.h"
 #include "blas.h"
+#include <algorithm>
+extern "C" {
 #include "clapack.h"
+	// to get atlas to link
+	float _sqrtf(float x) {return sqrt(x);}
+}
 
-extern "C" {
+template<typename T, typename K> 
+inline int lu_factor(int m, T a[], int ipiv[], 
+						 int (*getrf) (CBLAS_ORDER, const int, const int, K*, const int, int*))
+{
+	int info = getrf(CblasColMajor, m, m, a, m, ipiv);
+	shift_ipiv_down(m, ipiv);
+	return info;
+};
 
-	DLLEXPORT int s_cholesky_factor(int n, float a[]){
+template<typename T, typename K> 
-		int info = clapack_spotrf(CblasColMajor, CblasLower, n, a, n);
+inline int lu_inverse(int n, T a[],  
-		for (int i = 0; i < n; ++i)
+						  int (*getrf) (CBLAS_ORDER, const int, const int, K*, const int, int*),
+						  int (*getri) (CBLAS_ORDER, const int, K*, const int, const int*))
 {
-			int index = i * n;
+	int* ipiv = new int[n];
-			for (int j = 0; j < n && i > j; ++j)
+	int info = getrf(CblasColMajor, n, n, a, n, ipiv);
+
+	if (info != 0){
+		delete[] ipiv;
+		return info;
+	}
+
+	info = getri(CblasColMajor, n, a, n, ipiv);
+	delete[] ipiv;
+	return info;
+};
+
+template<typename T, typename K> 
+inline int lu_inverse_factored(int n, T a[], int ipiv[],
+								   int (*getri) (CBLAS_ORDER, const int, K*, const int, const int*))
 {
-				a[index + j] = 0;
+	shift_ipiv_up(n, ipiv);
+	int info = getri(CblasColMajor,n, a, n, ipiv);
+	shift_ipiv_down(n, ipiv);
+	return info;
 }
+
+template<typename T, typename K> 
+inline int lu_solve_factored(int n, int nrhs, T a[], int ipiv[], T b[],
+									int (*getrs) (CBLAS_ORDER, CBLAS_TRANSPOSE, const int, const int, const K*, const int, const int*, K*, const int))
+{
+	shift_ipiv_up(n, ipiv);
+	int info = getrs(CblasColMajor, CblasNoTrans, n, nrhs, a, n, ipiv, b, n);
+	shift_ipiv_down(n, ipiv);
+	return info;
 }
+
+template<typename T, typename K> 
+inline int lu_solve(int n, int nrhs, T a[], T b[], 
+						  int (*getrf) (CBLAS_ORDER, const int, const int, K*, const int, int*),
+						  int (*getrs) (CBLAS_ORDER, CBLAS_TRANSPOSE, const int, const int, const K*, const int, const int*, K*, const int))
+{
+	T* clone = Clone(n, n, a);
+	int* ipiv = new int[n];
+	int info = getrf(CblasColMajor, n, n, clone, n, ipiv);
+
+	if (info != 0){
+		delete[] ipiv;
+		delete[] clone;
 		return info;
 	}
 
-	DLLEXPORT int d_cholesky_factor(int n, double* a){
+	info = getrs(CblasColMajor, CblasNoTrans, n, nrhs, clone, n, ipiv, b, n);
-		int info = clapack_dpotrf(CblasColMajor, CblasLower, n, a, n);
+	delete[] ipiv;
+	delete[] clone;
+	return info;
+}
+
+template<typename T, typename K> 
+inline int cholesky_factor(int n, T* a, int (*potrf) (CBLAS_ORDER, CBLAS_UPLO, const int, K*, const int))
+{
+	int info = potrf(CblasColMajor, CblasLower, n, a, n);
+	T zero = T();
 	for (int i = 0; i < n; ++i)
 	{
 		int index = i * n;
 		for (int j = 0; j < n && i > j; ++j)
 		{
-				a[index + j] = 0;
+			a[index + j] = zero;
 		}
 	}
 	return info;
 }
 
-	DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){
+template<typename T, typename K> 
-		int info = clapack_cpotrf(CblasColMajor, CblasLower, n, a, n);
+inline int cholesky_solve(int n, int nrhs, T a[], T b[],
-		Complex8 zero;
+							  int (*potrf) (CBLAS_ORDER, CBLAS_UPLO, const int, K*, const int),
-		zero.real = 0.0;
+							  int (*potrs) (CBLAS_ORDER, CBLAS_UPLO, const int, const int, const K*, const int, K*, const int))
-		zero.real = 0.0;
-		for (int i = 0; i < n; ++i)
 {
-			int index = i * n;
+	T* clone = Clone(n, n, a);
-			for (int j = 0; j < n && i > j; ++j)
+	int info = potrf(CblasColMajor, CblasLower, n, clone, n);
+
+	if (info != 0){
+		delete[] clone;
+		return info;
+	}
+
+	info = potrs(CblasColMajor, CblasLower, n, nrhs, clone, n, b, n);
+	delete[] clone;
+	return info;
+}
+
+template<typename T, typename K> 
+inline int cholesky_solve_factored(int n, int nrhs, T a[], T b[],
+									   int (*potrs) (CBLAS_ORDER, CBLAS_UPLO, const int, const int, const K*, const int, K*, const int))
 {
-				a[index + j] = zero;
+	return potrs(CblasColMajor, CblasLower, n, nrhs, a, n, b, n);
+}
+
+template<typename T, typename K> 
+inline int qr_factor(int m, int n, T r[], T tau[], T q[], T work[], int len,
+						 int (*geqrf) (const int, const int, K*, const int, T*),
+						 int (*orgqr) (const int, const int, const int, K*, const int, const K*))
+{
+	int info = geqrf(m, n, r, m, tau);
+
+	for (int i = 0; i < m; ++i)
+	{
+		for (int j = 0; j < m && j < n; ++j)
+		{
+			if (i > j)
+			{
+				q[j * m + i] = r[j * m + i];
 			}
 		}
+	}
+
+	//compute the q elements explicitly
+	if (m <= n)
+	{
+		info = orgqr(m, m, m, q, m, tau);
+	}
+	else
+	{
+		info = orgqr(m, m, n, q, m, tau);
+	}
+
 	return info;
 }
 
-	DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){
+template<typename T> 
-		int info = clapack_zpotrf(CblasColMajor, CblasLower, n, a, n);
+inline int qr_thin_factor(int m, int n, T q[], T tau[], T r[], T work[], int len,
-		Complex16 zero;
+						      void (*geqrf) (const int*, const int*, T*, const int*, T*, T*, const int*, int*),
-		zero.real = 0.0;
+						      void (*orgqr) (const int*, const int*, const int*, T*, const int*, const T*, T*, const int*, int*))
-		zero.real = 0.0;
+{
+	int info = 0;
+	geqrf(&m, &n, q, &m, tau, work, &len, &info);
+
 	for (int i = 0; i < n; ++i)
 	{
-			int index = i * n;
+		for (int j = 0; j < n; ++j)
-			for (int j = 0; j < n && i > j; ++j)
 		{
-				a[index + j] = zero;
+			if( i <= j) {
+				r[j * n + i] = q[j * m + i];
+			}
+		}
 	}
+
+	orgqr(&m, &n, &n, q, &m, tau, work, &len, &info);
+
+	return info;
+}
+
+template<typename T> 
+inline int qr_solve(int m, int n, int bn, T a[], T b[], T x[], T work[], int len,
+						void (*gels) (const char*, const int*, const int*, const int*, T*, 
+						const int*, T* b, const int*, T*, const int*, int*))
+{
+	T* clone_a = new T[m*n];
+	std::memcpy(clone_a, a, m*n*sizeof(T));
+
+	T* clone_b = new T[m*bn];
+	std::memcpy(clone_b, b, m*bn*sizeof(T));
+
+	char N = 'N';
+	int info = 0;
+	gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info);
+	copyBtoX(n, n, bn, clone_b, x);
+
+	delete[] clone_a;
+	delete[] clone_b;
+	return info;
+}
+
+template<typename T> 
+inline int qr_solve_factored(int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len,
+								   void (*ormqr) (const char*, const char*, const int*, const int*, const int*, 
+								                  const T*, const int*, const T*, T*, const int*, T*, const int*, int* info),
+								   void (*trsm) (const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG, 
+								                 const int, const int, const T, const T*, const int, T*, const int)) 
+{
+	T* clone_b = new T[m*bn];
+	std::memcpy(clone_b, b, m*bn*sizeof(T));
+
+	char side ='L';
+	char tran = 'T';
+	int info = 0;
+	ormqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info);
+	trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m);
+	copyBtoX(n, n, bn, clone_b, x);
+
+	delete[] clone_b;
+	return info;
+}
+
+template<typename T> 
+inline int complex_qr_solve_factored(int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len,
+								   void (*unmqr) (const char*, const char*, const int*, const int*, const int*, 
+								                  const T*, const int*, const T*, T*, const int*, T*, const int*, int* info),
+								   void (*trsm) (const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG, 
+								                 const int, const int, const void*, const void*, const int, void*, const int ldb)) 
+{
+	T* clone_b = new T[m*bn];
+	std::memcpy(clone_b, b, m*bn*sizeof(T));
+
+	char side ='L';
+	char tran = 'C';
+	int info = 0;
+	unmqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info);
+
+	T one = {1.0f, 0.0f};
+	trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m);
+	copyBtoX(n, n, bn, clone_b, x);
+
+	delete[] clone_b;
+	return info;
+}
+
+template<typename T> 
+inline int svd_factor(bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], T work[], int len,
+								  void (*gesvd) (const char*, const char*, const int*, const int*, T*, const int*, 
+								                 T*, T*, const int*, T*, const int*, T*, const int*, int*))
+{
+	int info = 0;
+	char job = compute_vectors ? 'A' : 'N';
+	gesvd(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info);
+	return info;
+}
+
+
+template<typename T, typename R> 
+inline int complex_svd_factor(bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], T work[], int len,
+								  void (*gesvd) (const char*, const char*, const int*, const int*, T*, const int*, 
+								                 R*, T*, const int*, T*, const int*, T*, const int*, R*, int*))
+{
+	int info = 0;
+	int dim_s = std::min(m,n);
+	R* rwork = new R[5 * dim_s];
+	R* s_local = new R[dim_s];
+	char job = compute_vectors ? 'A' : 'N';
+	gesvd(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info);
+
+	for(int index = 0; index < dim_s; ++index){
+		T value = {s_local[index], 0.0f};
+		s[index] = value;
 	}
+
+	delete[] rwork;
+	delete[] s_local;
 	return info;
 }
+
+extern "C" {
+	DLLEXPORT int s_lu_factor(int m, float a[], int ipiv[]) {
+		return lu_factor<float, float>(m, a, ipiv, clapack_sgetrf);
+	}
+	
+	DLLEXPORT int d_lu_factor(int m, double a[], int ipiv[]) {
+		return lu_factor<double, double>(m, a, ipiv, clapack_dgetrf);
+	}
+	
+	DLLEXPORT int c_lu_factor(int m, Complex8 a[], int ipiv[]) {
+		return lu_factor<Complex8, void>(m, a, ipiv, clapack_cgetrf);
+	}
+	
+	DLLEXPORT int z_lu_factor(int m, Complex16 a[], int ipiv[]) {
+		return lu_factor(m, a, ipiv, clapack_zgetrf);
+	}
+
+	DLLEXPORT int s_lu_inverse(int n, float a[])
+	{
+		return lu_inverse<float, float>(n, a, clapack_sgetrf, clapack_sgetri);
+	}
+
+	DLLEXPORT int d_lu_inverse(int n, double a[])
+	{
+		return lu_inverse<double, double>(n, a, clapack_dgetrf, clapack_dgetri);
+	}
+
+	DLLEXPORT int c_lu_inverse(int n, Complex8 a[])
+	{
+		return lu_inverse<Complex8, void>(n, a, clapack_cgetrf, clapack_cgetri);
+	}
+
+	DLLEXPORT int z_lu_inverse(int n, Complex16 a[])
+	{
+		return lu_inverse<Complex16, void>(n, a, clapack_zgetrf, clapack_zgetri);
+	}
+
+	DLLEXPORT int s_lu_inverse_factored(int n, float a[], int ipiv[], float work[], int lwork)
+	{
+		return lu_inverse_factored<float, float>(n, a, ipiv, clapack_sgetri);
+	}
+
+	DLLEXPORT int d_lu_inverse_factored(int n, double a[], int ipiv[], double work[], int lwork)
+	{
+		return lu_inverse_factored<double, double>(n, a, ipiv, clapack_dgetri);
+	}
+
+	DLLEXPORT int c_lu_inverse_factored(int n, Complex8 a[], int ipiv[], Complex8 work[], int lwork)
+	{
+		return lu_inverse_factored<Complex8, void>(n, a, ipiv, clapack_cgetri);
+	}
+
+	DLLEXPORT int z_lu_inverse_factored(int n, Complex16 a[], int ipiv[], Complex16 work[], int lwork)
+	{
+		return lu_inverse_factored<Complex16, void>(n, a, ipiv, clapack_zgetri);
+	}
+
+	DLLEXPORT int s_lu_solve_factored(int n, int nrhs, float a[], int ipiv[], float b[])
+	{
+		return lu_solve_factored<float, float>(n, nrhs, a, ipiv, b, clapack_sgetrs);
+	}
+
+	DLLEXPORT int  d_lu_solve_factored(int n, int nrhs, double a[], int ipiv[], double b[])
+	{
+		return lu_solve_factored<double, double>(n, nrhs, a, ipiv, b, clapack_dgetrs);
+	}
+
+	DLLEXPORT int c_lu_solve_factored(int n, int nrhs, Complex8 a[], int ipiv[], Complex8 b[])
+	{
+		return lu_solve_factored<Complex8, void>(n, nrhs, a, ipiv, b, clapack_cgetrs);
+	}
+
+	DLLEXPORT int z_lu_solve_factored(int n, int nrhs, Complex16 a[], int ipiv[], Complex16 b[])
+	{
+		return lu_solve_factored<Complex16, void>(n, nrhs, a, ipiv, b, clapack_zgetrs);
+	}
+
+	DLLEXPORT int s_lu_solve(int n, int nrhs, float a[], float b[])
+	{
+		return lu_solve<float, float>(n, nrhs, a, b, clapack_sgetrf, clapack_sgetrs);
+	}
+
+	DLLEXPORT int d_lu_solve(int n, int nrhs, double a[], double b[])
+	{
+		return lu_solve<double, double>(n, nrhs, a, b, clapack_dgetrf, clapack_dgetrs);
+	}
+
+	DLLEXPORT int c_lu_solve(int n, int nrhs, Complex8 a[], Complex8 b[])
+	{
+		return lu_solve<Complex8, void>(n, nrhs, a, b, clapack_cgetrf, clapack_cgetrs);
+	}
+	
+	DLLEXPORT int z_lu_solve(int n, int nrhs, Complex16 a[],  Complex16 b[])
+	{
+		return lu_solve<Complex16, void>(n, nrhs, a, b, clapack_zgetrf, clapack_zgetrs);
+	}
+
+	DLLEXPORT int s_cholesky_factor(int n, float a[]){
+		return cholesky_factor<float, float>(n, a, clapack_spotrf);
+	}
+
+	DLLEXPORT int d_cholesky_factor(int n, double* a){
+		return cholesky_factor<double, double>(n, a, clapack_dpotrf);
+	}
+
+	DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){
+		return cholesky_factor<Complex8, void>(n, a, clapack_cpotrf);
+	}
+
+	DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){
+		return cholesky_factor<Complex16, void>(n, a, clapack_zpotrf);
+	}
+
+	DLLEXPORT int s_cholesky_solve(int n, int nrhs, float a[], float b[])
+	{
+		return cholesky_solve<float, float>(n, nrhs, a, b, clapack_spotrf, clapack_spotrs);
+	}
+
+	DLLEXPORT int d_cholesky_solve(int n, int nrhs, double a[], double b[])
+	{
+		return cholesky_solve<double, double>(n, nrhs, a, b, clapack_dpotrf, clapack_dpotrs);
+	}
+
+	DLLEXPORT int c_cholesky_solve(int n, int nrhs, Complex8 a[], Complex8 b[])
+	{
+		return cholesky_solve<Complex8, void>(n, nrhs, a, b, clapack_cpotrf, clapack_cpotrs);
+	}
+
+	DLLEXPORT int z_cholesky_solve(int n, int nrhs, Complex16 a[], Complex16 b[])
+	{
+		return cholesky_solve<Complex16, void>(n, nrhs, a, b, clapack_zpotrf, clapack_zpotrs);
+	}
+
+	DLLEXPORT int s_cholesky_solve_factored(int n, int nrhs, float a[], float b[])
+	{
+		return cholesky_solve_factored<float, float>(n, nrhs, a, b, clapack_spotrs);
+	}
+
+	DLLEXPORT int d_cholesky_solve_factored(int n, int nrhs, double a[], double b[])
+	{
+		return cholesky_solve_factored<double, double>(n, nrhs, a, b, clapack_dpotrs);
+	}
+
+	DLLEXPORT int c_cholesky_solve_factored(int n, int nrhs, Complex8 a[], Complex8 b[])
+	{
+		return cholesky_solve_factored<Complex8, void>(n, nrhs, a, b, clapack_cpotrs);
+	}
+
+	DLLEXPORT int z_cholesky_solve_factored(int n, int nrhs, Complex16 a[], Complex16 b[])
+	{
+		return cholesky_solve_factored<Complex16, void>(n, nrhs, a, b, clapack_zpotrs);
+	}
+
+	/*DLLEXPORT int s_qr_factor(int m, int n, float r[], float tau[], float q[], float work[], int len)
+	{
+		return qr_factor<float, float>(m, n, r, tau, q, work, len, clapack_sgeqrf, clapack_sorgqr);
+	}
+
+	DLLEXPORT int s_qr_thin_factor(int m, int n, float q[], float tau[], float r[], float work[], int len)
+	{
+		return qr_thin_factor<float>(m, n, q, tau, r, work, len, clapack_sgeqrf, clapack_sorgqr);
+	}
+
+	DLLEXPORT int d_qr_factor(int m, int n, double r[], double tau[], double q[], double work[], int len)
+	{
+		return qr_factor<double>(m, n, r, tau, q, work, len, clapack_dgeqrf, clapack_dorgqr);
+	}
+
+	DLLEXPORT int d_qr_thin_factor(int m, int n, double q[], double tau[], double r[], double work[], int len)
+	{
+		return qr_thin_factor<double>(m, n, q, tau, r, work, len, clapack_dgeqrf, clapack_dorgqr);
+	}
+
+	DLLEXPORT int c_qr_factor(int m, int n, Complex8 r[], Complex8 tau[], Complex8 q[], Complex8 work[], int len)
+	{
+		return qr_factor<Complex8>(m, n, r, tau, q, work, len, clapack_cgeqrf, clapack_cungqr);
+	}
+
+	DLLEXPORT int c_qr_thin_factor(int m, int n, Complex8 q[], Complex8 tau[], Complex8 r[], Complex8 work[], int len)
+	{
+		return qr_thin_factor<Complex8>(m, n, q, tau, r, work, len, clapack_cgeqrf, clapack_cungqr);
+	}
+
+	DLLEXPORT int z_qr_factor(int m, int n, Complex16 r[], Complex16 tau[], Complex16 q[])
+	{
+		return qr_factor<Complex16>(m, n, r, tau, q, work, len, clapack_zgeqrf, clapack_zungqr);
+	}
+
+	DLLEXPORT int z_qr_thin_factor(int m, int n, Complex16 q[], Complex16 tau[], Complex16 r[])
+	{
+		return qr_thin_factor<Complex16>(m, n, q, tau, r, work, len, clapack_zgeqrf, clapack_zungqr);
+	}
+
+	DLLEXPORT int s_qr_solve(int m, int n, int bn, float a[], float b[], float x[], float work[], int len)
+	{
+		return qr_solve<float>(m, n, bn, a, b, x, work, len, sgels);
+	}
+
+	DLLEXPORT int d_qr_solve(int m, int n, int bn, double a[], double b[], double x[], double work[], int len)
+	{
+		return qr_solve<double>(m, n, bn, a, b, x, work, len, dgels);
+	}
+
+	DLLEXPORT int c_qr_solve(int m, int n, int bn, Complex8 a[], Complex8 b[], Complex8 x[], Complex8 work[], int len)
+	{
+		return qr_solve<Complex8>(m, n, bn, a, b, x, work, len, cgels);
+	}
+
+	DLLEXPORT int z_qr_solve(int m, int n, int bn, Complex16 a[], Complex16 b[], Complex16 x[], Complex16 work[], int len)
+	{
+		return qr_solve<Complex16>(m, n, bn, a, b, x, work, len, zgels);
+	}
+
+	DLLEXPORT int s_qr_solve_factored(int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len)
+	{
+		return qr_solve_factored<float>(m, n, bn, r, b, tau, x, work, len, sormqr, cblas_strsm);	
+	}
+
+	DLLEXPORT int d_qr_solve_factored(int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len)
+	{
+		return qr_solve_factored<double>(m, n, bn, r, b, tau, x, work, len, dormqr, cblas_dtrsm);	
+	}
+
+	DLLEXPORT int c_qr_solve_factored(int m, int n, int bn, Complex8 r[], Complex8 b[], Complex8 tau[], Complex8 x[], Complex8 work[], int len)
+	{
+		return complex_qr_solve_factored<Complex8>(m, n, bn, r, b, tau, x, work, len, cunmqr, cblas_ctrsm);	
+	}
+
+	DLLEXPORT int z_qr_solve_factored(int m, int n, int bn, Complex16 r[], Complex16 b[], Complex16 tau[], Complex16 x[], Complex16 work[], int len)
+	{
+		return complex_qr_solve_factored<Complex16>(m, n, bn, r, b, tau, x, work, len, zunmqr, cblas_ztrsm);	
+	}
+
+	DLLEXPORT int s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[], float work[], int len)
+	{
+		return svd_factor<float>(compute_vectors, m, n, a, s, u, v, work, len, sgesvd);
+	}
+
+	DLLEXPORT int d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[], double work[], int len)
+	{
+		return svd_factor<double>(compute_vectors, m, n, a, s, u, v, work, len, dgesvd);
+	}
+
+	DLLEXPORT int c_svd_factor(bool compute_vectors, int m, int n, Complex8 a[], Complex8 s[], Complex8 u[], Complex8 v[], Complex8 work[], int len)
+	{
+		return complex_svd_factor<Complex8, float>(compute_vectors, m, n, a, s, u, v, work, len, cgesvd);
+	}
+	
+	DLLEXPORT int z_svd_factor(bool compute_vectors, int m, int n, Complex16 a[], Complex16 s[], Complex16 u[], Complex16 v[], Complex16 work[], int len)
+	{
+		return complex_svd_factor<Complex16, double>(compute_vectors, m, n, a, s, u, v, work, len, zgesvd);
+	}*/
 }

src/NativeWrappers/Common/lapack_common.h

@@ -18,8 +18,7 @@ inline void shift_ipiv_up(int m, int ipiv[]){
 }
 
 template<typename T> 
-inline T* Clone(const int m, const int n, const T* a)
+inline T* Clone(const int m, const int n, const T* a){
-{
 	T* clone = new T[m*n];
 	memcpy(clone, a, m*n*sizeof(T));
 	return clone;

src/NativeWrappers/Common/resource.rc

@@ -51,8 +51,8 @@ END
 //
 
 VS_VERSION_INFO VERSIONINFO
- FILEVERSION 1,2,1,0
+ FILEVERSION 1,3,0,0
- PRODUCTVERSION 1,2,1,0
+ PRODUCTVERSION 1,3,0,0
  FILEFLAGSMASK 0x17L
 #ifdef _DEBUG
  FILEFLAGS 0x1L
@@ -70,12 +70,12 @@ BEGIN
             VALUE "Comments", "http://numerics.mathdotnet.com/"
             VALUE "CompanyName", "Math.NET"
             VALUE "FileDescription", "MathNET Numerics Native Wrapper"
-            VALUE "FileVersion", "1.2.1.0"
+            VALUE "FileVersion", "1.3.0.0"
             VALUE "InternalName", "Math.NET"
             VALUE "LegalCopyright", "Copyright (C) Math.NET 2009-2013"
             VALUE "OriginalFilename", "MathNet.Numerics"
             VALUE "ProductName", "Math.NET Numerics"
-            VALUE "ProductVersion", "1.2.1.0"
+            VALUE "ProductVersion", "1.3.0.0"
         END
     END
     BLOCK "VarFileInfo"

src/NativeWrappers/Linux/build.sh

@@ -1,6 +1,11 @@
 export INTEL=/opt/intel
 export MKL=$INTEL/mkl
+export OPENMP=$INTEL/composerxe/lib
 
-g++ --shared -fPIC -o ./x64/MathNet.Numerics.MKL.so -I$MKL/include -I../Common  ../MKL/vector_functions.c ../MKL/blas.c ../MKL/lapack.cpp $INTEL/lib/intel64/libiomp5.a $MKL/lib/intel64/libmkl_intel_lp64.a $MKL/lib/intel64/libmkl_intel_thread.a $MKL/lib/intel64/libmkl_core.a 
+g++ -DGCC -m64 --shared -fPIC -o ../../../../MKL/Linux/x64/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common ../MKL/vector_functions.c ../MKL/blas.c ../MKL/lapack.cpp -Wl,--start-group  $MKL/lib/intel64/libmkl_intel_lp64.a $MKL/lib/intel64/libmkl_intel_thread.a $MKL/lib/intel64/libmkl_core.a -Wl,--end-group -L$OPENMP/intel64 -liomp5 -lpthread -lm  
 
-g++ -m32 --shared -fPIC -o ./x86/MathNet.Numerics.MKL.so -I$MKL/include -I../Common  ../MKL/vector_functions.c ../MKL/blas.c ../MKL/lapack.cpp $INTEL/lib/ia32/libiomp5.a $MKL/lib/ia32/libmkl_intel.a $MKL/lib/ia32/libmkl_intel_thread.a $MKL/lib/ia32/libmkl_core.a 
+cp $OPENMP/intel64/libiomp5.so  ../../../../MKL/Linux/x64/
+
+g++ -DGCC -m32 --shared -fPIC -o ../../../../MKL/Linux/x86/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common ../MKL/vector_functions.c ../MKL/blas.c ../MKL/lapack.cpp -Wl,--start-group $MKL/lib/ia32/libmkl_intel.a $MKL/lib/ia32/libmkl_intel_thread.a $MKL/lib/ia32/libmkl_core.a -Wl,--end-group -L$OPENMP/ia32 -liomp5 -lpthread -lm  
+
+cp $OPENMP/ia32/libiomp5.so  ../../../../MKL/Linux/x86/

src/NativeWrappers/MKL/blas.c

@@ -1,6 +1,9 @@
 #include "mkl_cblas.h"
 #include "wrapper_common.h"
 
+#if GCC 
+extern "C" { 
+#endif
 DLLEXPORT void s_axpy(const MKL_INT n, const float alpha, const float x[], float y[]){
 	cblas_saxpy(n, alpha, x, 1, y, 1);
 }
@@ -80,3 +83,7 @@ DLLEXPORT void z_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB,
 
 	cblas_zgemm(CblasColMajor, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m);
 }
+
+#if GCC 
+} 
+#endif

src/NativeWrappers/MKL/lapack.cpp

@@ -1,13 +1,19 @@
-#include "mkl_lapack.h"
+#include <algorithm>
+#include <complex>
+#define MKL_Complex8 std::complex<float>
+#define MKL_Complex16 std::complex<double>
+
+#include "mkl_lapack.h"
 #include "mkl_cblas.h"
 #include "lapack_common.h"
 #include "wrapper_common.h"
-#include <algorithm>
+#include "mkl_lapacke.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*))
 {
+    std::complex<double> x = 5;
 	MKL_INT info = 0;
     getrf(&m, &m, a, &m, ipiv, &info);
     shift_ipiv_down(m, ipiv);
@@ -23,7 +29,8 @@ inline MKL_INT lu_inverse(MKL_INT n, T a[], T work[], MKL_INT lwork,
     MKL_INT info = 0;
     getrf(&n, &n, a, &n, ipiv, &info);
 
-	if (info != 0){
+    if (info != 0)
+    {
         delete[] ipiv;
         return info;
     }
@@ -66,7 +73,8 @@ inline MKL_INT lu_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
     MKL_INT info = 0;
     getrf(&n, &n, clone, &n, ipiv, &info);
 
-	if (info != 0){
+    if (info != 0)
+    {
         delete[] ipiv;
         delete[] clone;
         return info;
@@ -88,14 +96,17 @@ inline MKL_INT cholesky_factor(MKL_INT n, T* a,
     MKL_INT info = 0;
     potrf(&uplo, &n, a, &n, &info);
     T zero = T();
+
     for (MKL_INT i = 0; i < n; ++i)
     {
         MKL_INT index = i * n;
+
         for (MKL_INT j = 0; j < n && i > j; ++j)
         {
             a[index + j] = zero;
         }
     }
+
     return info;
 }
 
@@ -109,7 +120,8 @@ inline MKL_INT cholesky_solve(MKL_INT n, MKL_INT nrhs, T a[], T b[],
     MKL_INT info = 0;
     potrf(&uplo, &n, clone, &n, &info);
 
-	if (info != 0){
+    if (info != 0)
+    {
         delete[] clone;
         return info;
     }
@@ -173,14 +185,14 @@ inline MKL_INT qr_thin_factor(MKL_INT m, MKL_INT n, T q[], T tau[], T r[], T wor
     {
         for (MKL_INT j = 0; j < n; ++j)
         {
-			if( i <= j) {
+            if (i <= j)
+            {
                 r[j * n + i] = q[j * m + i];
             }
         }
     }
 
     orgqr(&m, &n, &n, q, &m, tau, work, &len, &info);
-
     return info;
 }
 
@@ -191,19 +203,10 @@ inline MKL_INT qr_solve(MKL_INT m, MKL_INT n, MKL_INT bn, T a[], T b[], T x[], T
 {
     T* clone_a = Clone(m, n, a);
     T* clone_b = Clone(m, bn, b);
-
     char N = 'N';
     MKL_INT info = 0;
     gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info);
-
+    copyBtoX(m, n, bn, clone_b, x);
-	for (MKL_INT i = 0; i < n; ++i)
-	{
-		for (MKL_INT j = 0; j < bn; ++j)
-		{
-			x[j * n + i] = clone_b[j * m + i];
-		}
-	}
-
     delete[] clone_a;
     delete[] clone_b;
     return info;
@@ -222,14 +225,7 @@ inline MKL_INT qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[], T b[],
     MKL_INT info = 0;
     ormqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info);
     trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m);
-	for (MKL_INT i = 0; i < n; ++i)
+    copyBtoX(m, n, bn, clone_b, x);
-	{
-		for (MKL_INT j = 0; j < bn; ++j)
-		{
-			x[j * n + i] = clone_b[j * m + i];
-		}
-	}
-
     delete[] clone_b;
     return info;
 }
@@ -246,17 +242,9 @@ inline MKL_INT complex_qr_solve_factored(MKL_INT m, MKL_INT n, MKL_INT bn, T r[]
     char tran = 'C';
     MKL_INT info = 0;
     unmqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info);
-
+	T one = 1.0f;
-	T one = {1.0f, 0.0f};
     trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m);
-	for (MKL_INT i = 0; i < n; ++i)
+    copyBtoX(m, n, bn, clone_b, x);
-	{
-		for (MKL_INT j = 0; j < bn; ++j)
-		{
-			x[j * n + i] = clone_b[j * m + i];
-		}
-	}
-
     delete[] clone_b;
     return info;
 }
@@ -284,9 +272,9 @@ inline MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T
     char job = compute_vectors ? 'A' : 'N';
     gesvd(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info);
 
-	for(MKL_INT index = 0; index < dim_s; ++index){
+    for (MKL_INT index = 0; index < dim_s; ++index)
-		T value = {s_local[index], 0.0f};
+    {
-		s[index] = value;
+        s[index] = s_local[index];
     }
 
     delete[] rwork;
@@ -294,6 +282,142 @@ 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*))
+{
+    T* clone_a = Clone(n, n, a);
+    T* wr = new T[n];
+    T* wi = new T[n];
+
+	MKL_INT sdim;
+    MKL_INT info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, wr, wi, vectors, n);
+    if (info != 0)
+    {
+        delete[] clone_a;
+        delete[] wr;
+        delete[] wi;
+        return info;
+    }
+
+    MKL_INT m;
+    info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m);
+    if (info != 0)
+    {
+        delete[] clone_a;
+        delete[] wr;
+        delete[] wi;
+        return info;
+    }
+
+    for (MKL_INT index = 0; index < n; ++index)
+    {
+        values[index] = MKL_Complex16(wr[index], wi[index]);
+    }
+
+    for (MKL_INT  i = 0; i < n; ++i)
+    {
+        MKL_INT in = i * n;
+        d[in + i] = wr[i];
+
+        if (wi[i] > 0)
+        {
+            d[in + n + i] = wi[i];
+        }
+        else if (wi[i] < 0)
+        {
+            d[in - n + i] = wi[i];
+        }
+    }
+
+    delete[] clone_a;
+    delete[] wr;
+    delete[] wi;
+    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*))
+{
+	T* clone_a = Clone(n, n, a);
+	T* w = new T[n];
+
+	MKL_INT sdim;
+	MKL_INT info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, w, vectors, n);
+	if (info != 0)
+	{
+		delete[] clone_a;
+		delete[] w;
+		return info;
+	}
+
+	MKL_INT m;
+	info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m);
+	if (info != 0)
+	{
+		delete[] clone_a;
+		delete[] w;
+		return info;
+	}
+
+	for (MKL_INT i = 0; i < n; ++i)
+	{
+		values[i] = w[i];
+		d[i * n + i] = w[i];
+	}
+
+	delete[] clone_a;
+	delete[] w;
+	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*))
+{
+    T* clone_a = Clone(n, n, a);
+    R* w = new R[n];
+
+	MKL_INT info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w);
+	if (info != 0)
+	{
+		delete[] clone_a;
+		delete[] w;
+		return info;
+	}
+
+	memcpy(vectors, clone_a, n*n*sizeof(T));
+
+	for (MKL_INT index = 0; index < n; ++index)
+	{
+		values[index] = MKL_Complex16(w[index]);
+	}
+    
+	for (MKL_INT j = 0; j < n; j++)
+    {
+        MKL_INT jn = j*n;
+
+        for (MKL_INT i = 0; i < n; ++i)
+        {
+            if (i == j)
+            {
+                d[jn + i] = w[i];
+            }
+        }
+    }
+
+    delete[] clone_a;
+    delete[] w;
+    return info;
+}
+
 extern "C" {
 
     DLLEXPORT float s_matrix_norm(char norm, MKL_INT m, MKL_INT n, float a[], float work[])
@@ -316,19 +440,23 @@ extern "C" {
         return zlange(&norm, &m, &n, a, &m, work);
     }
 
-	DLLEXPORT MKL_INT s_lu_factor(MKL_INT m, float a[], MKL_INT ipiv[]) {
+    DLLEXPORT MKL_INT s_lu_factor(MKL_INT m, float a[], MKL_INT ipiv[])
+    {
         return lu_factor<float>(m, a, ipiv, sgetrf);
     }
 
-	DLLEXPORT MKL_INT d_lu_factor(MKL_INT m, double a[], MKL_INT ipiv[]) {
+    DLLEXPORT MKL_INT d_lu_factor(MKL_INT m, double a[], MKL_INT ipiv[])
+    {
         return lu_factor<double>(m, a, ipiv, dgetrf);
     }
 
-	DLLEXPORT MKL_INT c_lu_factor(MKL_INT m, MKL_Complex8 a[], MKL_INT ipiv[]) {
+    DLLEXPORT MKL_INT c_lu_factor(MKL_INT m, MKL_Complex8 a[], MKL_INT ipiv[])
+    {
         return lu_factor<MKL_Complex8>(m, a, ipiv, cgetrf);
     }
 
-	DLLEXPORT MKL_INT z_lu_factor(MKL_INT m, MKL_Complex16 a[], MKL_INT ipiv[]) {
+    DLLEXPORT MKL_INT z_lu_factor(MKL_INT m, MKL_Complex16 a[], MKL_INT ipiv[])
+    {
         return lu_factor<MKL_Complex16>(m, a, ipiv, zgetrf);
     }
 
@@ -412,19 +540,23 @@ extern "C" {
         return lu_solve<MKL_Complex16>(n, nrhs, a, b, zgetrf, zgetrs);
     }
 
-	DLLEXPORT MKL_INT s_cholesky_factor(MKL_INT n, float a[]){
+    DLLEXPORT MKL_INT s_cholesky_factor(MKL_INT n, float a[])
+    {
         return cholesky_factor<float>(n, a, spotrf);
     }
 
-	DLLEXPORT MKL_INT d_cholesky_factor(MKL_INT n, double* a){
+    DLLEXPORT MKL_INT d_cholesky_factor(MKL_INT n, double* a)
+    {
         return cholesky_factor<double>(n, a, dpotrf);
     }
 
-	DLLEXPORT MKL_INT c_cholesky_factor(MKL_INT n, MKL_Complex8 a[]){
+    DLLEXPORT MKL_INT c_cholesky_factor(MKL_INT n, MKL_Complex8 a[])
+    {
         return cholesky_factor<MKL_Complex8>(n, a, cpotrf);
     }
 
-	DLLEXPORT MKL_INT z_cholesky_factor(MKL_INT n, MKL_Complex16 a[]){
+    DLLEXPORT MKL_INT z_cholesky_factor(MKL_INT n, MKL_Complex16 a[])
+    {
         return cholesky_factor<MKL_Complex16>(n, a, zpotrf);
     }
 
@@ -567,4 +699,53 @@ extern "C" {
     {
         return complex_svd_factor<MKL_Complex16, double>(compute_vectors, m, n, a, s, u, v, work, len, zgesvd);
     }
+
+    DLLEXPORT MKL_INT s_eigen(bool isSymmetric, MKL_INT n, float a[], float vectors[], MKL_Complex16 values[], float d[])
+    {
+		if (isSymmetric)
+		{
+			return sym_eigen_factor<float, float>(n, a, vectors, values, d, LAPACKE_ssyev);
+		}
+		else
+		{
+			return eigen_factor<float>(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[])
+    {
+        if (isSymmetric)
+        {
+            return sym_eigen_factor<double, double>(n, a, vectors, values, d, LAPACKE_dsyev);
+        }
+        else
+        {
+            return eigen_factor<double>(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc);
+        }
+    }
+
+    DLLEXPORT MKL_INT c_eigen(bool isSymmetric, MKL_INT n, MKL_Complex8 a[], MKL_Complex8 vectors[], MKL_Complex16 values[], MKL_Complex8 d[])
+    {
+		if (isSymmetric)
+		{
+			return sym_eigen_factor<MKL_Complex8, 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);
+		}
+    }
+
+    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);
+		}
+		else
+		{
+			return eigen_complex_factor<MKL_Complex16>(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc);
+		}
+    }
 }

src/NativeWrappers/MKL/vector_functions.c

@@ -1,7 +1,9 @@
 #include "mkl_vml.h"
 #include "wrapper_common.h"
 
-
+#if GCC
+extern "C" {
+#endif
 DLLEXPORT void s_vector_add( const int n, const float x[], const float y[], float result[] ){
 	vsAdd( n, x, y, result );
 }
@@ -65,3 +67,6 @@ DLLEXPORT void z_vector_multiply( const int n, const MKL_Complex16 x[], const MK
 DLLEXPORT void z_vector_divide( const int n, const MKL_Complex16 x[], const MKL_Complex16 y[], MKL_Complex16 result[] ){
 	vzDiv( n, x, y, result );
 }
+#if GCC
+}
+#endif

src/NativeWrappers/Windows/ATLASWrapper/ATLASWrapper.vcxproj

@@ -0,0 +1,160 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project DefaultTargets="Build" ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup Label="ProjectConfigurations">
+    <ProjectConfiguration Include="Debug|Win32">
+      <Configuration>Debug</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Debug|x64">
+      <Configuration>Debug</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|Win32">
+      <Configuration>Release</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|x64">
+      <Configuration>Release</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+  </ItemGroup>
+  <PropertyGroup Label="Globals">
+    <ProjectGuid>{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}</ProjectGuid>
+    <Keyword>Win32Proj</Keyword>
+    <RootNamespace>ATLASWrapper</RootNamespace>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'" Label="Configuration">
+    <ConfigurationType>DynamicLibrary</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v110</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
+    <ConfigurationType>DynamicLibrary</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v110</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
+    <ConfigurationType>DynamicLibrary</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v110</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
+    <ConfigurationType>DynamicLibrary</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v110</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
+  <ImportGroup Label="ExtensionSettings">
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="PropertySheets">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="PropertySheets">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <PropertyGroup Label="UserMacros" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <LinkIncremental>true</LinkIncremental>
+    <IncludePath>D:\Source\ATLAS\include;..\..\Common;$(IncludePath)</IncludePath>
+    <OutDir>$(SolutionDir)..\..\..\..\ATLAS\Windows\x86\</OutDir>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <LinkIncremental>true</LinkIncremental>
+    <IncludePath>D:\Source\ATLAS\include;..\..\Common;$(IncludePath)</IncludePath>
+    <OutDir>$(SolutionDir)..\..\..\..\ATLAS\Windows\x64\</OutDir>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <LinkIncremental>false</LinkIncremental>
+    <IncludePath>D:\Source\ATLAS\include;..\..\Common;$(IncludePath)</IncludePath>
+    <OutDir>$(SolutionDir)..\..\..\..\ATLAS\Windows\x86\</OutDir>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <LinkIncremental>false</LinkIncremental>
+    <IncludePath>D:\Source\ATLAS\include;..\..\Common;$(IncludePath)</IncludePath>
+    <OutDir>$(SolutionDir)..\..\..\..\ATLAS\Windows\x64\</OutDir>
+  </PropertyGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <ClCompile>
+      <PrecompiledHeader>NotUsing</PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>Disabled</Optimization>
+      <PreprocessorDefinitions>WIN32;_DEBUG;_WINDOWS;_USRDLL;ATLASWRAPPER_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+    </ClCompile>
+    <Link>
+      <SubSystem>Windows</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalDependencies>libptcblas.a;libatlas.a;libptlapack.a;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <ClCompile>
+      <PrecompiledHeader>NotUsing</PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>Disabled</Optimization>
+      <PreprocessorDefinitions>WIN32;_DEBUG;_WINDOWS;_USRDLL;ATLASWRAPPER_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+    </ClCompile>
+    <Link>
+      <SubSystem>Windows</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalDependencies>libptcblas.a;libatlas.a;libptlapack.a;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <PrecompiledHeader>NotUsing</PrecompiledHeader>
+      <Optimization>MaxSpeed</Optimization>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <PreprocessorDefinitions>WIN32;NDEBUG;_WINDOWS;_USRDLL;ATLASWRAPPER_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+    </ClCompile>
+    <Link>
+      <SubSystem>Windows</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <AdditionalDependencies>libptcblas.a;libatlas.a;libptlapack.a;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <PrecompiledHeader>NotUsing</PrecompiledHeader>
+      <Optimization>MaxSpeed</Optimization>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <PreprocessorDefinitions>WIN32;NDEBUG;_WINDOWS;_USRDLL;ATLASWRAPPER_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+    </ClCompile>
+    <Link>
+      <SubSystem>Windows</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <AdditionalDependencies>libptcblas.a;libatlas.a;libptlapack.a;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemGroup>
+    <ResourceCompile Include="..\..\Common\resource.rc" />
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="..\..\ATLAS\blas.c" />
+    <ClCompile Include="..\..\ATLAS\lapack.cpp" />
+    <ClCompile Include="..\..\Common\WindowsDLL.cpp" />
+  </ItemGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
+  <ImportGroup Label="ExtensionTargets">
+  </ImportGroup>
+</Project>

src/NativeWrappers/Windows/ATLASWrapper/ATLASWrapper.vcxproj.filters

@@ -0,0 +1,33 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup>
+    <Filter Include="Source Files">
+      <UniqueIdentifier>{4FC737F1-C7A5-4376-A066-2A32D752A2FF}</UniqueIdentifier>
+      <Extensions>cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx</Extensions>
+    </Filter>
+    <Filter Include="Header Files">
+      <UniqueIdentifier>{93995380-89BD-4b04-88EB-625FBE52EBFB}</UniqueIdentifier>
+      <Extensions>h;hpp;hxx;hm;inl;inc;xsd</Extensions>
+    </Filter>
+    <Filter Include="Resource Files">
+      <UniqueIdentifier>{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}</UniqueIdentifier>
+      <Extensions>rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms</Extensions>
+    </Filter>
+  </ItemGroup>
+  <ItemGroup>
+    <ResourceCompile Include="..\..\Common\resource.rc">
+      <Filter>Resource Files</Filter>
+    </ResourceCompile>
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="..\..\Common\WindowsDLL.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="..\..\ATLAS\blas.c">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="..\..\ATLAS\lapack.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+  </ItemGroup>
+</Project>

src/NativeWrappers/Windows/MKL/MKLWrapper.vcxproj

@@ -241,6 +241,7 @@
       <WarningLevel>Level3</WarningLevel>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
       <CompileAs>Default</CompileAs>
+      <AdditionalOptions>/Qvec-report:1 %(AdditionalOptions)</AdditionalOptions>
     </ClCompile>
     <Link>
       <AdditionalDependencies>libiomp5md.lib;mkl_intel_c.lib;mkl_intel_thread.lib;mkl_core.lib;%(AdditionalDependencies)</AdditionalDependencies>
@@ -271,6 +272,7 @@
       <WarningLevel>Level3</WarningLevel>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
       <CompileAs>Default</CompileAs>
+      <AdditionalOptions>/Qvec-report:1 %(AdditionalOptions)</AdditionalOptions>
     </ClCompile>
     <Link>
       <AdditionalDependencies>libiomp5md.lib;mkl_intel_lp64.lib;mkl_intel_thread.lib;mkl_core.lib;%(AdditionalDependencies)</AdditionalDependencies>
@@ -292,9 +294,6 @@
     <ClCompile Include="..\..\MKL\lapack.cpp" />
     <ClCompile Include="..\..\MKL\vector_functions.c" />
   </ItemGroup>
-  <ItemGroup>
-    <ClInclude Include="..\..\Common\wrapper_common.h" />
-  </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="..\..\Common\resource.rc" />
   </ItemGroup>

src/NativeWrappers/Windows/MKL/MKLWrapper.vcxproj.filters

@@ -28,11 +28,6 @@
       <Filter>Source Files</Filter>
     </ClCompile>
   </ItemGroup>
-  <ItemGroup>
-    <ClInclude Include="..\..\Common\wrapper_common.h">
-      <Filter>Header Files</Filter>
-    </ClInclude>
-  </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="..\..\Common\resource.rc">
       <Filter>Resource Files</Filter>

src/NativeWrappers/Windows/NativeWrappers.sln

@@ -12,22 +12,40 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Common", "Common", "{5A0892
 EndProject
 Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "MKLWrapper", "MKL\MKLWrapper.vcxproj", "{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}"
 EndProject
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "ATLASWrapper", "ATLASWrapper\ATLASWrapper.vcxproj", "{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}"
+EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|Mixed Platforms = Debug|Mixed Platforms
 		Debug|Win32 = Debug|Win32
 		Debug|x64 = Debug|x64
+		Release|Mixed Platforms = Release|Mixed Platforms
 		Release|Win32 = Release|Win32
 		Release|x64 = Release|x64
 	EndGlobalSection
 	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Debug|Mixed Platforms.ActiveCfg = Debug|Win32
+		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Debug|Mixed Platforms.Build.0 = Debug|Win32
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Debug|Win32.ActiveCfg = Debug|Win32
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Debug|Win32.Build.0 = Debug|Win32
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Debug|x64.ActiveCfg = Debug|x64
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Debug|x64.Build.0 = Debug|x64
+		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Mixed Platforms.ActiveCfg = Release|Win32
+		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Mixed Platforms.Build.0 = Release|Win32
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Win32.ActiveCfg = Release|Win32
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Win32.Build.0 = Release|Win32
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.ActiveCfg = Release|x64
 		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.Build.0 = Release|x64
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Debug|Mixed Platforms.ActiveCfg = Debug|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Debug|Mixed Platforms.Build.0 = Debug|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Debug|Win32.ActiveCfg = Debug|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Debug|Win32.Build.0 = Debug|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Debug|x64.ActiveCfg = Debug|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|Mixed Platforms.ActiveCfg = Release|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|Mixed Platforms.Build.0 = Release|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|Win32.ActiveCfg = Release|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|Win32.Build.0 = Release|Win32
+		{2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|x64.ActiveCfg = Release|x64
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE

src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProvider.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -93,5 +93,49 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// The requested <see cref="Norm"/> of the matrix.
         /// </returns>
         Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix, double[] work);
+
+        /// <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>
+        void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD);
+
+        /// <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>
+        void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD);
+
+        /// <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>
+        void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD);
+
+        /// <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>
+        void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD);
     }
 }

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs

@@ -3,7 +3,7 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -24,14 +24,14 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
-
 namespace MathNet.Numerics.Algorithms.LinearAlgebra
 {
     using System;
     using System.Numerics;
     using Properties;
     using Threading;
+    using Numerics.LinearAlgebra.Complex.Factorization;
+    using Numerics.LinearAlgebra.Generic.Factorization;
 
     /// <summary>
     /// The managed linear algebra provider.
@@ -720,7 +720,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="constN">The constant number of columns of matrix op(B) and of the matrix C.</param>
         /// <param name="constK">The constant number of columns of matrix op(A) and the rows of the matrix op(B).</param>
         /// <param name="first">Indicates if this is the first recursion.</param>
-        private static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] matrixA, int shiftArow, int shiftAcol, Complex[] matrixB, int shiftBrow, int shiftBcol, Complex[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first)
+        static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] matrixA, int shiftArow, int shiftAcol, Complex[] matrixB, int shiftBrow, int shiftBcol, Complex[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first)
         {
             if (m + n <= Control.ParallelizeOrder)
             {
@@ -1345,7 +1345,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="colLimit">Total columns</param>
         /// <param name="multipliers">Multipliers calculated previously</param>
         /// <param name="availableCores">Number of available processors</param>
-        private static void DoCholeskyStep(Complex[] data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores)
+        static void DoCholeskyStep(Complex[] data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores)
         {
             var tmpColCount = colLimit - firstCol;
 
@@ -1453,7 +1453,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="orderA">The number of rows and columns in A.</param>
         /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
         /// <param name="index">The column to solve for.</param>
-        private static void DoCholeskySolve(Complex[] a, int orderA, Complex[] b, int index)
+        static void DoCholeskySolve(Complex[] a, int orderA, Complex[] b, int index)
         {
             var cindex = index*orderA;
 
@@ -1757,7 +1757,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnStart">The first column</param>
         /// <param name="columnCount">The last column</param>
         /// <param name="availableCores">Number of available CPUs</param>
-        private static void ComputeQR(Complex[] work, int workIndex, Complex[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
+        static void ComputeQR(Complex[] work, int workIndex, Complex[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
         {
             if (rowStart > rowCount || columnStart > columnCount)
             {
@@ -1801,7 +1801,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="rowCount">The number of rows in matrix</param>
         /// <param name="row">The first row</param>
         /// <param name="column">Column index</param>
-        private static void GenerateColumn(Complex[] work, Complex[] a, int rowCount, int row, int column)
+        static void GenerateColumn(Complex[] work, Complex[] a, int rowCount, int row, int column)
         {
             var tmp = column*rowCount;
             var index = tmp + row;
@@ -2987,5 +2987,85 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 }
             }
         }
+
+        /// <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 virtual 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");
+            }
+            var matrixCopy = new Complex[matrix.Length];
+            Array.Copy(matrix, matrixCopy, matrix.Length);
+
+            if (isSymmetric)
+            {
+                var tau = new Complex[order];
+                var d = new double[order];
+                var e = new double[order];
+
+                DenseEvd.SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
+                DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order);
+                DenseEvd.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order);
+
+                for (var i = 0; i < order; i++)
+                {
+                    vectorEv[i] = new Complex(d[i], e[i]);
+                }
+            }
+            else
+            {
+                DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order);
+                DenseEvd.NonsymmetricReduceHessenberToRealSchur(vectorEv, matrixEv, matrixCopy, order);
+            }
+
+            for (var i = 0; i < order; i ++)
+            {
+                matrixD[i * order + i] = vectorEv[i];
+            }
+        }
     }
 }

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs

@@ -3,7 +3,7 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -24,14 +24,15 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using System.Numerics;
-using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
-
 namespace MathNet.Numerics.Algorithms.LinearAlgebra
 {
     using System;
     using Properties;
     using Threading;
+    using System.Numerics;
+    using Numerics.LinearAlgebra.Complex32.Factorization;
+    using Numerics.LinearAlgebra.Generic.Factorization;
+    using Numerics.LinearAlgebra.Complex32;
 
     /// <summary>
     /// The managed linear algebra provider.
@@ -2984,5 +2985,87 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 }
             }
         }
+
+        /// <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 virtual 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");
+            }
+
+            var matrixCopy = new Complex32[matrix.Length];
+            Array.Copy(matrix, matrixCopy, matrix.Length);
+            var v = new DenseVector(order);
+            if (isSymmetric)
+            {
+                var tau = new Complex32[order];
+                var d = new float[order];
+                var e = new float[order];
+
+                DenseEvd.SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
+                DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order);
+                DenseEvd.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order);
+
+                for (var i = 0; i < order; i++)
+                {
+                    vectorEv[i] = new Complex(d[i], e[i]);
+                    matrixD[i * order + i] = new Complex32(d[i], e[i]);
+                }
+            }
+            else
+            {
+                DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order);
+                DenseEvd.NonsymmetricReduceHessenberToRealSchur(v.Values, matrixEv, matrixCopy, order);
+                for (var i = 0; i < order; i++)
+                {
+                    vectorEv[i] = new Complex(v[i].Real, v[i].Imaginary);
+                    matrixD[i * order + i] = v[i];
+                }
+            }
+        }
     }
 }

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs

@@ -3,7 +3,7 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -24,11 +24,11 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
-
 namespace MathNet.Numerics.Algorithms.LinearAlgebra
 {
     using System;
+    using System.Numerics;
+    using Numerics.LinearAlgebra.Generic.Factorization;
     using Properties;
     using Threading;
 
@@ -2933,5 +2933,98 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 }
             }
         }
+
+        /// <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 virtual void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] 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");
+            }
+
+            var d = new double[order];
+            var e = new double[order];
+
+            if (isSymmetric)
+            {
+                Buffer.BlockCopy(matrix, 0, matrixEv, 0, matrix.Length * Constants.SizeOfDouble);
+                var om1 = order - 1;
+                for (var i = 0; i < order; i++)
+                {
+                    d[i] = matrixEv[i*order + om1];
+                }
+
+                Numerics.LinearAlgebra.Double.Factorization.DenseEvd.SymmetricTridiagonalize(matrixEv, d, e, order);
+                Numerics.LinearAlgebra.Double.Factorization.DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order);
+            }
+            else
+            {
+                var matrixH = new double[matrix.Length];
+                Buffer.BlockCopy(matrix, 0, matrixH, 0, matrix.Length * Constants.SizeOfDouble);
+                Numerics.LinearAlgebra.Double.Factorization.DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order);
+                Numerics.LinearAlgebra.Double.Factorization.DenseEvd.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order);
+            }
+
+            for (var i = 0; i < order; i++)
+            {
+                vectorEv[i] = new Complex(d[i], e[i]);
+
+                var io = i * order;
+                matrixD[io + i] = d[i];
+
+                if (e[i] > 0)
+                {
+                    matrixD[io + order + i] = e[i];
+                    matrixD[(i+1) * order + i] = e[i];
+                }
+                else if (e[i] < 0)
+                {
+                    matrixD[io - order + i] = e[i];
+                }
+            }
+        }
     }
 }

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs

@@ -3,7 +3,7 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -24,11 +24,11 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
-
 namespace MathNet.Numerics.Algorithms.LinearAlgebra
 {
     using System;
+    using System.Numerics;
+    using Numerics.LinearAlgebra.Generic.Factorization;
     using Properties;
     using Threading;
 
@@ -2936,5 +2936,97 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 }
             }
         }
+
+        /// <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 virtual void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] 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");
+            }
+
+            var d = new float[order];
+            var e = new float[order];
+
+            if (isSymmetric)
+            {
+                Buffer.BlockCopy(matrix, 0, matrixEv, 0, matrix.Length * Constants.SizeOfFloat);
+                var om1 = order - 1;
+                for (var i = 0; i < order; i++)
+                {
+                    d[i] = matrixEv[i * order + om1];
+                }
+
+                Numerics.LinearAlgebra.Single.Factorization.DenseEvd.SymmetricTridiagonalize(matrixEv, d, e, order);
+                Numerics.LinearAlgebra.Single.Factorization.DenseEvd.SymmetricDiagonalize(matrixEv, d, e, order);
+            }
+            else
+            {
+                var matrixH = new float[matrix.Length];
+                Buffer.BlockCopy(matrix, 0, matrixH, 0, matrix.Length * Constants.SizeOfFloat);
+                Numerics.LinearAlgebra.Single.Factorization.DenseEvd.NonsymmetricReduceToHessenberg(matrixEv, matrixH, order);
+                Numerics.LinearAlgebra.Single.Factorization.DenseEvd.NonsymmetricReduceHessenberToRealSchur(matrixEv, matrixH, d, e, order);
+            }
+
+            for (var i = 0; i < order; i++)
+            {
+                vectorEv[i] = new Complex(d[i], e[i]);
+
+                var io = i * order;
+                matrixD[io + i] = d[i];
+
+                if (e[i] > 0)
+                {
+                    matrixD[io + order + i] = e[i];
+                }
+                else if (e[i] < 0)
+                {
+                    matrixD[io - order + i] = e[i];
+                }
+            }
+        }
     }
 }

src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2011 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation

src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2011 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,12 +28,11 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
-
 namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
 {
     using System;
     using System.Security;
+    using Numerics.LinearAlgebra.Generic.Factorization;
     using Properties;
 
     /// <summary>

src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.double.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2011 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,12 +28,12 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
-
 namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
 {
     using System;
+    using System.Numerics;
     using System.Security;
+    using Numerics.LinearAlgebra.Generic.Factorization;
     using Properties;
 
     /// <summary>
@@ -700,7 +700,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
 
             var work = new double[columnsA*Control.BlockSize];
             SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length);
-
         }
 
 
@@ -1277,5 +1276,59 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
 
             SafeNativeMethods.d_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, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] 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");
+            }
+
+            SafeNativeMethods.d_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD);
+        }
     }
 }

src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.float.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2011 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,16 +28,12 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-/* This file is automatically generated - do not modify it. 
-   Last generated on UTC 2011-04-17 06:45:23Z
-*/
-
-using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
-
 namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
 {
     using System;
+    using System.Numerics;
     using System.Security;
+    using Numerics.LinearAlgebra.Generic.Factorization;
     using Properties;
 
     /// <summary>
@@ -1177,5 +1173,59 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
 
             SafeNativeMethods.s_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, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] 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");
+            }
+
+            SafeNativeMethods.s_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD);
+        }
     }
 }

src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs

@@ -2,7 +2,7 @@
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -266,6 +266,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
         internal static extern int z_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v, [In, Out] Complex[] work, int len);
 
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int s_eigen(bool isSymmetric, int n, [In] float[] a, [In, Out] float[] vectors, [In, Out] Complex[] values, [In, Out] float[] d);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int d_eigen(bool isSymmetric, int n, [In] double[] a, [In, Out] double[] vectors, [In, Out]  Complex[] values, [In, Out] double[] d);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int c_eigen(bool isSymmetric, int n, [In] Complex32[] a, [In, Out] Complex32[] vectors, [In, Out]  Complex[] values, [In, Out] Complex32[] d);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int z_eigen(bool isSymmetric, int n, [In] Complex[] a, [In, Out] Complex[] vectors, [In, Out]  Complex[] values, [In, Out] Complex[] d);
+        
         #endregion LAPACK
 
         #region Vector Functions

src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -27,10 +27,10 @@
 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
+
 namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 {
     using System;
-    using System.Numerics;
     using Generic;
     using Properties;
     
@@ -72,45 +72,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 
             var order = matrix.RowCount;
 
-            // Initialize matricies for eigenvalues and eigenvectors
+            // Initialize matrices for eigenvalues and eigenvectors
             MatrixEv = DenseMatrix.Identity(order);
             MatrixD = matrix.CreateMatrix(order, order);
             VectorEv = new DenseVector(order);
 
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                 }
             }
 
-            if (IsSymmetric)
+            Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) MatrixEv).Values,
-            {
+                ((DenseVector) VectorEv).Values, ((DenseMatrix) MatrixD).Values);
-                var matrixCopy = matrix.ToArray();
-                var tau = new Complex[order];
-                var d = new double[order];
-                var e = new double[order];
-
-                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
-                SymmetricDiagonalize(((DenseMatrix)MatrixEv).Values, d, e, order);
-                SymmetricUntridiagonalize(((DenseMatrix)MatrixEv).Values, matrixCopy, tau, order);
-
-                for (var i = 0; i < order; i++)
-                {
-                    VectorEv[i] = new Complex(d[i], e[i]);
-                }
-            }
-            else
-            {
-                var matrixH = matrix.ToArray();
-                NonsymmetricReduceToHessenberg(((DenseMatrix)MatrixEv).Values, matrixH, order);
-                NonsymmetricReduceHessenberToRealSchur(((DenseVector)VectorEv).Values, ((DenseMatrix)MatrixEv).Values, matrixH, order);
-            }
-
-            MatrixD.SetDiagonal(VectorEv);
         }
 
         /// <summary>
@@ -125,14 +103,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for 
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricTridiagonalize(Complex[,] matrixA, double[] d, double[] e, Complex[] tau, int order)
+        internal static void SymmetricTridiagonalize(System.Numerics.Complex[] matrixA, double[] d, double[] e, System.Numerics.Complex[] tau, int order)
         {
             double hh;
-            tau[order - 1] = Complex.One;
+            tau[order - 1] = System.Numerics.Complex.One;
 
             for (var i = 0; i < order; i++)
             {
-                d[i] = matrixA[i, i].Real;
+                d[i] = matrixA[i*order + i].Real;
             }
 
             // Householder reduction to tridiagonal form.
@@ -144,61 +122,62 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 
                 for (var k = 0; k < i; k++)
                 {
-                    scale = scale + Math.Abs(matrixA[i, k].Real) + Math.Abs(matrixA[i, k].Imaginary);
+                    scale = scale + Math.Abs(matrixA[k*order + i].Real) + Math.Abs(matrixA[k*order + i].Imaginary);
                 }
 
                 if (scale == 0.0)
                 {
-                    tau[i - 1] = Complex.One;
+                    tau[i - 1] = System.Numerics.Complex.One;
                     e[i] = 0.0;
                 }
                 else
                 {
                     for (var k = 0; k < i; k++)
                     {
-                        matrixA[i, k] /= scale;
+                        matrixA[k*order + i] /= scale;
-                        h += matrixA[i, k].MagnitudeSquared();
+                        h += matrixA[k*order + i].MagnitudeSquared();
                     }
 
-                    Complex g = Math.Sqrt(h);
+                    System.Numerics.Complex g = Math.Sqrt(h);
                     e[i] = scale*g.Real;
 
-                    Complex temp;
+                    System.Numerics.Complex temp;
-                    var f = matrixA[i, i - 1];
+                    var im1Oi = (i - 1)*order + i;
+                    var f = matrixA[im1Oi];
                     if (f.Magnitude != 0)
                     {
-                        temp = -(matrixA[i, i - 1].Conjugate() * tau[i].Conjugate()) / f.Magnitude;
+                        temp = -(matrixA[im1Oi].Conjugate()*tau[i].Conjugate())/f.Magnitude;
                         h += f.Magnitude*g.Real;
                         g = 1.0 + (g/f.Magnitude);
-                        matrixA[i, i - 1] *= g;
+                        matrixA[im1Oi] *= g;
                     }
                     else
                     {
                         temp = -tau[i].Conjugate();
-                        matrixA[i, i - 1] = g;
+                        matrixA[im1Oi] = g;
                     }
 
                     if ((f.Magnitude == 0) || (i != 1))
                     {
-                        f = Complex.Zero;
+                        f = System.Numerics.Complex.Zero;
                         for (var j = 0; j < i; j++)
                         {
-                            var tmp = Complex.Zero;
+                            var tmp = System.Numerics.Complex.Zero;
-
+                            var jO = j*order;
                             // Form element of A*U.
                             for (var k = 0; k <= j; k++)
                             {
-                                tmp += matrixA[j, k] * matrixA[i, k].Conjugate();
+                                tmp += matrixA[k*order + j]*matrixA[k*order + i].Conjugate();
                             }
 
                             for (var k = j + 1; k <= i - 1; k++)
                             {
-                                tmp += matrixA[k, j].Conjugate() * matrixA[i, k].Conjugate();
+                                tmp += matrixA[jO + k].Conjugate()*matrixA[k*order + i].Conjugate();
                             }
 
                             // Form element of P
                             tau[j] = tmp/h;
-                            f += (tmp / h) * matrixA[i, j];
+                            f += (tmp/h)*matrixA[jO + i];
                         }
 
                         hh = f.Real/(h + h);
@@ -206,33 +185,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                         // Form the reduced A.
                         for (var j = 0; j < i; j++)
                         {
-                            f = matrixA[i, j].Conjugate();
+                            f = matrixA[j*order + i].Conjugate();
                             g = tau[j] - (hh*f);
                             tau[j] = g.Conjugate();
 
                             for (var k = 0; k <= j; k++)
                             {
-                                matrixA[j, k] -= (f * tau[k]) + (g * matrixA[i, k]);
+                                matrixA[k*order + j] -= (f*tau[k]) + (g*matrixA[k*order + i]);
                             }
                         }
                     }
 
                     for (var k = 0; k < i; k++)
                     {
-                        matrixA[i, k] *= scale;
+                        matrixA[k*order + i] *= scale;
                     }
 
                     tau[i - 1] = temp.Conjugate();
                 }
 
                 hh = d[i];
-                d[i] = matrixA[i, i].Real;
+                d[i] = matrixA[i*order + i].Real;
-                matrixA[i, i] = new Complex(hh, scale * Math.Sqrt(h));
+                matrixA[i*order + i] = new System.Numerics.Complex(hh, scale*Math.Sqrt(h));
             }
 
             hh = d[0];
-            d[0] = matrixA[0, 0].Real;
+            d[0] = matrixA[0].Real;
-            matrixA[0, 0] = hh;
+            matrixA[0] = hh;
             e[0] = 0.0;
         }
 
@@ -247,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricDiagonalize(Complex[] dataEv, double[] d, double[] e, int order)
+        internal static void SymmetricDiagonalize(System.Numerics.Complex[] dataEv, double[] d, double[] e, int order)
         {
             const int Maxiter = 1000;
 
@@ -347,8 +326,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                         {
                             throw new ArgumentException(Resources.ConvergenceFailed);
                         }
-                    }
+                    } while (Math.Abs(e[l]) > eps*tst1);
-                    while (Math.Abs(e[l]) > eps * tst1);
                 }
 
                 d[l] = d[l] + f;
@@ -394,7 +372,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// by Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricUntridiagonalize(Complex[] dataEv, Complex[,] matrixA, Complex[] tau, int order)
+        internal static void SymmetricUntridiagonalize(System.Numerics.Complex[] dataEv, System.Numerics.Complex[] matrixA, System.Numerics.Complex[] tau, int order)
         {
             for (var i = 0; i < order; i++)
             {
@@ -407,22 +385,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             // Recover and apply the Householder matrices.
             for (var i = 1; i < order; i++)
             {
-                var h = matrixA[i, i].Imaginary;
+                var h = matrixA[i*order + i].Imaginary;
                 if (h != 0)
                 {
                     for (var j = 0; j < order; j++)
                     {
-                        var s = Complex.Zero;
+                        var s = System.Numerics.Complex.Zero;
                         for (var k = 0; k < i; k++)
                         {
-                            s += dataEv[(j * order) + k] * matrixA[i, k];
+                            s += dataEv[(j*order) + k]*matrixA[k*order + i];
                         }
 
                         s = (s/h)/h;
 
                         for (var k = 0; k < i; k++)
                         {
-                            dataEv[(j * order) + k] -= s * matrixA[i, k].Conjugate();
+                            dataEv[(j*order) + k] -= s*matrixA[k*order + i].Conjugate();
                         }
                     }
                 }
@@ -439,17 +417,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutines in EISPACK.</remarks>
-        private static void NonsymmetricReduceToHessenberg(Complex[] dataEv, Complex[,] matrixH, int order)
+        internal static void NonsymmetricReduceToHessenberg(System.Numerics.Complex[] dataEv, System.Numerics.Complex[] matrixH, int order)
         {
-            var ort = new Complex[order];
+            var ort = new System.Numerics.Complex[order];
 
             for (var m = 1; m < order - 1; m++)
             {
                 // Scale column.
                 var scale = 0.0;
+                var mm1O = (m - 1)*order;
                 for (var i = m; i < order; i++)
                 {
-                    scale += Math.Abs(matrixH[i, m - 1].Real) + Math.Abs(matrixH[i, m - 1].Imaginary);
+                    scale += Math.Abs(matrixH[mm1O + i].Real) + Math.Abs(matrixH[mm1O + i].Imaginary);
                 }
 
                 if (scale != 0.0)
@@ -458,7 +437,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     var h = 0.0;
                     for (var i = order - 1; i >= m; i--)
                     {
-                        ort[i] = matrixH[i, m - 1] / scale;
+                        ort[i] = matrixH[mm1O + i]/scale;
                         h += ort[i].MagnitudeSquared();
                     }
 
@@ -472,43 +451,44 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     else
                     {
                         ort[m] = g;
-                        matrixH[m, m - 1] = scale;
+                        matrixH[mm1O + m] = scale;
                     }
 
                     // Apply Householder similarity transformation
                     // H = (I-u*u'/h)*H*(I-u*u')/h)
                     for (var j = m; j < order; j++)
                     {
-                        var f = Complex.Zero;
+                        var f = System.Numerics.Complex.Zero;
+                        var jO = j*order;
                         for (var i = order - 1; i >= m; i--)
                         {
-                            f += ort[i].Conjugate() * matrixH[i, j];
+                            f += ort[i].Conjugate()*matrixH[jO + i];
                         }
 
                         f = f/h;
                         for (var i = m; i < order; i++)
                         {
-                            matrixH[i, j] -= f * ort[i];
+                            matrixH[jO + i] -= f*ort[i];
                         }
                     }
 
                     for (var i = 0; i < order; i++)
                     {
-                        var f = Complex.Zero;
+                        var f = System.Numerics.Complex.Zero;
                         for (var j = order - 1; j >= m; j--)
                         {
-                            f += ort[j] * matrixH[i, j];
+                            f += ort[j]*matrixH[j*order + i];
                         }
 
                         f = f/h;
                         for (var j = m; j < order; j++)
                         {
-                            matrixH[i, j] -= f * ort[j].Conjugate();
+                            matrixH[j*order + i] -= f*ort[j].Conjugate();
                         }
                     }
 
                     ort[m] = scale*ort[m];
-                    matrixH[m, m - 1] *= -g;
+                    matrixH[mm1O + m] *= -g;
                 }
             }
 
@@ -517,24 +497,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             {
                 for (var j = 0; j < order; j++)
                 {
-                    dataEv[(j * order) + i] = i == j ? Complex.One : Complex.Zero;
+                    dataEv[(j*order) + i] = i == j ? System.Numerics.Complex.One : System.Numerics.Complex.Zero;
                 }
             }
 
             for (var m = order - 2; m >= 1; m--)
             {
-                if (matrixH[m, m - 1] != Complex.Zero && ort[m] != Complex.Zero)
+                var mm1O = (m - 1)*order;
+                var mm1Om = mm1O + m;
+                if (matrixH[mm1Om] != System.Numerics.Complex.Zero && ort[m] != System.Numerics.Complex.Zero)
                 {
-                    var norm = (matrixH[m, m - 1].Real * ort[m].Real) + (matrixH[m, m - 1].Imaginary * ort[m].Imaginary);
+                    var norm = (matrixH[mm1Om].Real*ort[m].Real) + (matrixH[mm1Om].Imaginary*ort[m].Imaginary);
 
                     for (var i = m + 1; i < order; i++)
                     {
-                        ort[i] = matrixH[i, m - 1];
+                        ort[i] = matrixH[mm1O + i];
                     }
 
                     for (var j = m; j < order; j++)
                     {
-                        var g = Complex.Zero;
+                        var g = System.Numerics.Complex.Zero;
                         for (var i = m; i < order; i++)
                         {
                             g += ort[i].Conjugate()*dataEv[(j*order) + i];
@@ -553,18 +535,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             // Create real subdiagonal elements.
             for (var i = 1; i < order; i++)
             {
-                if (matrixH[i, i - 1].Imaginary != 0.0)
+                var im1 = i - 1;
+                var im1O = im1*order;
+                var im1Oi = im1O + i;
+                var iO = i*order;
+                if (matrixH[im1Oi].Imaginary != 0.0)
                 {
-                    var y = matrixH[i, i - 1] / matrixH[i, i - 1].Magnitude;
+                    var y = matrixH[im1Oi]/matrixH[im1Oi].Magnitude;
-                    matrixH[i, i - 1] = matrixH[i, i - 1].Magnitude;
+                    matrixH[im1Oi] = matrixH[im1Oi].Magnitude;
                     for (var j = i; j < order; j++)
                     {
-                        matrixH[i, j] *= y.Conjugate();
+                        matrixH[j*order + i] *= y.Conjugate();
                     }
 
                     for (var j = 0; j <= Math.Min(i + 1, order - 1); j++)
                     {
-                        matrixH[j, i] *= y;
+                        matrixH[iO + j] *= y;
                     }
 
                     for (var j = 0; j < order; j++)
@@ -586,14 +572,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, Complex[] dataEv, Complex[,] matrixH, int order)
+        internal static void NonsymmetricReduceHessenberToRealSchur(System.Numerics.Complex[] vectorV, System.Numerics.Complex[] dataEv, System.Numerics.Complex[] matrixH, int order)
         {
             // Initialize
             var n = order - 1;
             var eps = Precision.DoubleMachinePrecision;
 
             double norm;
-            Complex x, y, z, exshift = Complex.Zero;
+            System.Numerics.Complex x, y, z, exshift = System.Numerics.Complex.Zero;
 
             // Outer loop over eigenvalue index
             var iter = 0;
@@ -603,8 +589,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 var l = n;
                 while (l > 0)
                 {
-                    var tst1 = Math.Abs(matrixH[l - 1, l - 1].Real) + Math.Abs(matrixH[l - 1, l - 1].Imaginary) + Math.Abs(matrixH[l, l].Real) + Math.Abs(matrixH[l, l].Imaginary);
+                    var lm1 = l - 1;
-                    if (Math.Abs(matrixH[l, l - 1].Real) < eps * tst1)
+                    var lm1O = lm1*order;
+                    var lO = l*order;
+                    var tst1 = Math.Abs(matrixH[lm1O + lm1].Real) + Math.Abs(matrixH[lm1O + lm1].Imaginary) + Math.Abs(matrixH[lO + l].Real) + Math.Abs(matrixH[lO + l].Imaginary);
+                    if (Math.Abs(matrixH[lm1O + l].Real) < eps*tst1)
                     {
                         break;
                     }
@@ -612,27 +601,31 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     l--;
                 }
 
+                var nm1 = n - 1;
+                var nm1O = nm1*order;
+                var nO = n*order;
+                var nOn = nO + n;
                 // Check for convergence
                 // One root found
                 if (l == n)
                 {
-                    matrixH[n, n] += exshift;
+                    matrixH[nOn] += exshift;
-                    vectorV[n] = matrixH[n, n];
+                    vectorV[n] = matrixH[nOn];
                     n--;
                     iter = 0;
                 }
                 else
                 {
                     // Form shift
-                    Complex s;
+                    System.Numerics.Complex s;
                     if (iter != 10 && iter != 20)
                     {
-                        s = matrixH[n, n];
+                        s = matrixH[nOn];
-                        x = matrixH[n - 1, n] * matrixH[n, n - 1].Real;
+                        x = matrixH[nO + nm1]*matrixH[nm1O + n].Real;
 
                         if (x.Real != 0.0 || x.Imaginary != 0.0)
                         {
-                            y = (matrixH[n - 1, n - 1] - s) / 2.0;
+                            y = (matrixH[nm1O + nm1] - s)/2.0;
                             z = ((y*y) + x).SquareRoot();
                             if ((y.Real*z.Real) + (y.Imaginary*z.Imaginary) < 0.0)
                             {
@@ -646,12 +639,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     else
                     {
                         // Form exceptional shift
-                        s = Math.Abs(matrixH[n, n - 1].Real) + Math.Abs(matrixH[n - 1, n - 2].Real);
+                        s = Math.Abs(matrixH[nm1O + n].Real) + Math.Abs(matrixH[(n - 2)*order + nm1].Real);
                     }
 
                     for (var i = 0; i <= n; i++)
                     {
-                        matrixH[i, i] -= s;
+                        matrixH[i*order + i] -= s;
                     }
 
                     exshift += s;
@@ -660,31 +653,35 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     // Reduce to triangle (rows)
                     for (var i = l + 1; i <= n; i++)
                     {
-                        s = matrixH[i, i - 1].Real;
+                        var im1 = i - 1;
-                        norm = SpecialFunctions.Hypotenuse(matrixH[i - 1, i - 1].Magnitude, s.Real);
+                        var im1O = im1*order;
-                        x = matrixH[i - 1, i - 1] / norm;
+                        var im1Oim1 = im1O + im1;
+                        s = matrixH[im1O + i].Real;
+                        norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real);
+                        x = matrixH[im1Oim1]/norm;
                         vectorV[i - 1] = x;
-                        matrixH[i - 1, i - 1] = norm;
+                        matrixH[im1Oim1] = norm;
-                        matrixH[i, i - 1] = new Complex(0.0, s.Real / norm);
+                        matrixH[im1O + i] = new System.Numerics.Complex(0.0, s.Real/norm);
 
                         for (var j = i; j < order; j++)
                         {
-                            y = matrixH[i - 1, j];
+                            var jO = j*order;
-                            z = matrixH[i, j];
+                            y = matrixH[jO + im1];
-                            matrixH[i - 1, j] = (x.Conjugate() * y) + (matrixH[i, i - 1].Imaginary * z);
+                            z = matrixH[jO + i];
-                            matrixH[i, j] = (x * z) - (matrixH[i, i - 1].Imaginary * y);
+                            matrixH[jO + im1] = (x.Conjugate()*y) + (matrixH[im1O + i].Imaginary*z);
+                            matrixH[jO + i] = (x*z) - (matrixH[im1O + i].Imaginary*y);
                         }
                     }
 
-                    s = matrixH[n, n];
+                    s = matrixH[nOn];
                     if (s.Imaginary != 0.0)
                     {
-                        s /= matrixH[n, n].Magnitude;
+                        s /= matrixH[nOn].Magnitude;
-                        matrixH[n, n] = matrixH[n, n].Magnitude;
+                        matrixH[nOn] = matrixH[nOn].Magnitude;
 
                         for (var j = n + 1; j < order; j++)
                         {
-                            matrixH[n, j] *= s.Conjugate();
+                            matrixH[j*order + n] *= s.Conjugate();
                         }
                     }
 
@@ -692,29 +689,34 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     for (var j = l + 1; j <= n; j++)
                     {
                         x = vectorV[j - 1];
+                        var jO = j*order;
+                        var jm1 = j - 1;
+                        var jm1O = jm1*order;
+                        var jm1Oj = jm1O + j;
                         for (var i = 0; i <= j; i++)
                         {
-                            z = matrixH[i, j];
+                            var jm1Oi = jm1O + i;
+                            z = matrixH[jO + i];
                             if (i != j)
                             {
-                                y = matrixH[i, j - 1];
+                                y = matrixH[jm1Oi];
-                                matrixH[i, j - 1] = (x * y) + (matrixH[j, j - 1].Imaginary * z);
+                                matrixH[jm1Oi] = (x*y) + (matrixH[jm1O + j].Imaginary*z);
                             }
                             else
                             {
-                                y = matrixH[i, j - 1].Real;
+                                y = matrixH[jm1Oi].Real;
-                                matrixH[i, j - 1] = new Complex((x.Real * y.Real) - (x.Imaginary * y.Imaginary) + (matrixH[j, j - 1].Imaginary * z.Real), matrixH[i, j - 1].Imaginary);
+                                matrixH[jm1Oi] = new System.Numerics.Complex((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[jm1O + j].Imaginary*z.Real), matrixH[jm1Oi].Imaginary);
                             }
 
-                            matrixH[i, j] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y);
+                            matrixH[jO + i] = (x.Conjugate()*z) - (matrixH[jm1O + j].Imaginary*y);
                         }
 
                         for (var i = 0; i < order; i++)
                         {
                             y = dataEv[((j - 1)*order) + i];
                             z = dataEv[(j*order) + i];
-                            dataEv[((j - 1) * order) + i] = (x * y) + (matrixH[j, j - 1].Imaginary * z);
+                            dataEv[jm1O + i] = (x*y) + (matrixH[jm1Oj].Imaginary*z);
-                            dataEv[(j * order) + i] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y);
+                            dataEv[jO + i] = (x.Conjugate()*z) - (matrixH[jm1Oj].Imaginary*y);
                         }
                     }
 
@@ -722,12 +724,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     {
                         for (var i = 0; i <= n; i++)
                         {
-                            matrixH[i, n] *= s;
+                            matrixH[nO + i] *= s;
                         }
 
                         for (var i = 0; i < order; i++)
                         {
-                            dataEv[(n * order) + i] *= s;
+                            dataEv[nO + i] *= s;
                         }
                     }
                 }
@@ -740,7 +742,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             {
                 for (var j = i; j < order; j++)
                 {
-                    norm = Math.Max(norm, Math.Abs(matrixH[i, j].Real) + Math.Abs(matrixH[i, j].Imaginary));
+                    norm = Math.Max(norm, Math.Abs(matrixH[j*order + i].Real) + Math.Abs(matrixH[j*order + i].Imaginary));
                 }
             }
 
@@ -756,15 +758,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 
             for (n = order - 1; n > 0; n--)
             {
+                var nO = n*order;
+                var nOn = nO + n;
                 x = vectorV[n];
-                matrixH[n, n] = 1.0;
+                matrixH[nOn] = 1.0;
 
                 for (var i = n - 1; i >= 0; i--)
                 {
                     z = 0.0;
                     for (var j = i + 1; j <= n; j++)
                     {
-                        z += matrixH[i, j] * matrixH[j, n];
+                        z += matrixH[j*order + i]*matrixH[nO + j];
                     }
 
                     y = x - vectorV[i];
@@ -773,15 +777,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                         y = eps*norm;
                     }
 
-                    matrixH[i, n] = z / y;
+                    matrixH[nO + i] = z/y;
 
                     // Overflow control
-                    var tr = Math.Abs(matrixH[i, n].Real) + Math.Abs(matrixH[i, n].Imaginary);
+                    var tr = Math.Abs(matrixH[nO + i].Real) + Math.Abs(matrixH[nO + i].Imaginary);
                     if ((eps*tr)*tr > 1)
                     {
                         for (var j = i; j <= n; j++)
                         {
-                            matrixH[j, n] = matrixH[j, n] / tr;
+                            matrixH[nO + j] = matrixH[nO + j]/tr;
                         }
                     }
                 }
@@ -790,15 +794,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             // Back transformation to get eigenvectors of original matrix
             for (var j = order - 1; j > 0; j--)
             {
+                var jO = j*order;
                 for (var i = 0; i < order; i++)
                 {
-                    z = Complex.Zero;
+                    z = System.Numerics.Complex.Zero;
                     for (var k = 0; k <= j; k++)
                     {
-                        z += dataEv[(k * order) + i] * matrixH[k, j];
+                        z += dataEv[(k*order) + i]*matrixH[jO + k];
                     }
 
-                    dataEv[(j * order) + i] = z;
+                    dataEv[jO + i] = z;
                 }
             }
         }
@@ -808,7 +813,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// </summary>
         /// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
         /// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
-        public override void Solve(Matrix<Complex> input, Matrix<Complex> result)
+        public override void Solve(Matrix<System.Numerics.Complex> input, Matrix<System.Numerics.Complex> result)
         {
             // Check for proper arguments.
             if (input == null)
@@ -842,13 +847,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             if (IsSymmetric)
             {
                 var order = VectorEv.Count;
-                var tmp = new Complex[order];
+                var tmp = new System.Numerics.Complex[order];
 
                 for (var k = 0; k < order; k++)
                 {
                     for (var j = 0; j < order; j++)
                     {
-                        Complex value = 0.0;
+                        System.Numerics.Complex value = 0.0;
                         if (j < order)
                         {
                             for (var i = 0; i < order; i++)
@@ -864,7 +869,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 
                     for (var j = 0; j < order; j++)
                     {
-                        Complex value = 0.0;
+                        System.Numerics.Complex value = 0.0;
                         for (var i = 0; i < order; i++)
                         {
                             value += ((DenseMatrix) MatrixEv).Values[(i*order) + j]*tmp[i];
@@ -885,7 +890,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// </summary>
         /// <param name="input">The right hand side vector, <b>b</b>.</param>
         /// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
-        public override void Solve(Vector<Complex> input, Vector<Complex> result)
+        public override void Solve(Vector<System.Numerics.Complex> input, Vector<System.Numerics.Complex> result)
         {
             if (input == null)
             {
@@ -914,8 +919,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             {
                 // Symmetric case -> x = V * inv(λ) * VH * b;
                 var order = VectorEv.Count;
-                var tmp = new Complex[order];
+                var tmp = new System.Numerics.Complex[order];
-                Complex value;
+                System.Numerics.Complex value;
 
                 for (var j = 0; j < order; j++)
                 {
@@ -936,7 +941,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 for (var j = 0; j < order; j++)
                 {
                     value = 0;
-                    for (int i = 0; i < order; i++)
+                    for (var i = 0; i < order; i++)
                     {
                         value += ((DenseMatrix) MatrixEv).Values[(i*order) + j]*tmp[i];
                     }

src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs

@@ -79,9 +79,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
            
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                 }

src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -27,12 +27,12 @@
 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
+
+
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 {
     using System;
-    using System.Numerics;
     using Generic;
-    using Numerics;
     using Properties;
 
     /// <summary>
@@ -73,48 +73,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 
             var order = matrix.RowCount;
 
-            // Initialize matricies for eigenvalues and eigenvectors
+            // Initialize matrices for eigenvalues and eigenvectors
             MatrixEv = DenseMatrix.Identity(order);
             MatrixD = matrix.CreateMatrix(order, order);
-            VectorEv = new LinearAlgebra.Complex.DenseVector(order);
+            VectorEv = new Complex.DenseVector(order);
 
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                 }
             }
 
-            if (IsSymmetric)
+            Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) MatrixEv).Values,
-            {
+                ((Complex.DenseVector) VectorEv).Values, ((DenseMatrix) MatrixD).Values);
-                var matrixCopy = matrix.ToArray();
-                var tau = new Complex32[order];
-                var d = new float[order];
-                var e = new float[order];
-
-                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
-                SymmetricDiagonalize(((DenseMatrix)MatrixEv).Values, d, e, order);
-                SymmetricUntridiagonalize(((DenseMatrix)MatrixEv).Values, matrixCopy, tau, order);
-
-                for (var i = 0; i < order; i++)
-                {
-                    VectorEv[i] = new Complex(d[i], e[i]);
-                }
-            }
-            else
-            {
-                var matrixH = matrix.ToArray();
-                NonsymmetricReduceToHessenberg(((DenseMatrix)MatrixEv).Values, matrixH, order);
-                NonsymmetricReduceHessenberToRealSchur(((LinearAlgebra.Complex.DenseVector)VectorEv).Values, ((DenseMatrix)MatrixEv).Values, matrixH, order);
-            }
-
-            for (var i = 0; i < VectorEv.Count; i++)
-            {
-                MatrixD.At(i, i, (Complex32)VectorEv[i]);
-            }
         }
 
         /// <summary>
@@ -129,14 +104,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for 
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricTridiagonalize(Complex32[,] matrixA, float[] d, float[] e, Complex32[] tau, int order)
+        internal static void SymmetricTridiagonalize(Numerics.Complex32[] matrixA, float[] d, float[] e, Numerics.Complex32[] tau, int order)
         {
             float hh;
-            tau[order - 1] = Complex32.One;
+            tau[order - 1] = Numerics.Complex32.One;
 
             for (var i = 0; i < order; i++)
             {
-                d[i] = matrixA[i, i].Real;
+                d[i] = matrixA[i*order + i].Real;
             }
 
             // Householder reduction to tridiagonal form.
@@ -148,61 +123,62 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 
                 for (var k = 0; k < i; k++)
                 {
-                    scale = scale + Math.Abs(matrixA[i, k].Real) + Math.Abs(matrixA[i, k].Imaginary);
+                    scale = scale + Math.Abs(matrixA[k*order + i].Real) + Math.Abs(matrixA[k*order + i].Imaginary);
                 }
 
                 if (scale == 0.0f)
                 {
-                    tau[i - 1] = Complex32.One;
+                    tau[i - 1] = Numerics.Complex32.One;
                     e[i] = 0.0f;
                 }
                 else
                 {
                     for (var k = 0; k < i; k++)
                     {
-                        matrixA[i, k] /= scale;
+                        matrixA[k*order + i] /= scale;
-                        h += matrixA[i, k].MagnitudeSquared;
+                        h += matrixA[k*order + i].MagnitudeSquared;
                     }
 
-                    Complex32 g = (float)Math.Sqrt(h);
+                    Numerics.Complex32 g = (float) Math.Sqrt(h);
                     e[i] = scale*g.Real;
 
-                    Complex32 temp;
+                    Numerics.Complex32 temp;
-                    var f = matrixA[i, i - 1];
+                    var im1Oi = (i - 1)*order + i;
-                    if (f.Magnitude != 0)
+                    var f = matrixA[im1Oi];
+                    if (f.Magnitude != 0.0f)
                     {
-                        temp = -(matrixA[i, i - 1].Conjugate() * tau[i].Conjugate()) / f.Magnitude;
+                        temp = -(matrixA[im1Oi].Conjugate()*tau[i].Conjugate())/f.Magnitude;
                         h += f.Magnitude*g.Real;
                         g = 1.0f + (g/f.Magnitude);
-                        matrixA[i, i - 1] *= g;
+                        matrixA[im1Oi] *= g;
                     }
                     else
                     {
                         temp = -tau[i].Conjugate();
-                        matrixA[i, i - 1] = g;
+                        matrixA[im1Oi] = g;
                     }
 
-                    if ((f.Magnitude == 0) || (i != 1))
+                    if ((f.Magnitude == 0.0f) || (i != 1))
                     {
-                        f = Complex32.Zero;
+                        f = Numerics.Complex32.Zero;
                         for (var j = 0; j < i; j++)
                         {
-                            var tmp = Complex32.Zero;
+                            var tmp = Numerics.Complex32.Zero;
-
+                            var jO = j*order;
                             // Form element of A*U.
                             for (var k = 0; k <= j; k++)
                             {
-                                tmp += matrixA[j, k] * matrixA[i, k].Conjugate();
+                                tmp += matrixA[k*order + j]*matrixA[k*order + i].Conjugate();
                             }
 
                             for (var k = j + 1; k <= i - 1; k++)
                             {
-                                tmp += matrixA[k, j].Conjugate() * matrixA[i, k].Conjugate();
+                                tmp += matrixA[jO + k].Conjugate()*matrixA[k*order + i].Conjugate();
                             }
 
                             // Form element of P
                             tau[j] = tmp/h;
-                            f += (tmp / h) * matrixA[i, j];
+                            f += (tmp/h)*matrixA[jO + i];
                         }
 
                         hh = f.Real/(h + h);
@@ -210,33 +186,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                         // Form the reduced A.
                         for (var j = 0; j < i; j++)
                         {
-                            f = matrixA[i, j].Conjugate();
+                            f = matrixA[j*order + i].Conjugate();
                             g = tau[j] - (hh*f);
                             tau[j] = g.Conjugate();
 
                             for (var k = 0; k <= j; k++)
                             {
-                                matrixA[j, k] -= (f * tau[k]) + (g * matrixA[i, k]);
+                                matrixA[k*order + j] -= (f*tau[k]) + (g*matrixA[k*order + i]);
                             }
                         }
                     }
 
                     for (var k = 0; k < i; k++)
                     {
-                        matrixA[i, k] *= scale;
+                        matrixA[k*order + i] *= scale;
                     }
 
                     tau[i - 1] = temp.Conjugate();
                 }
 
                 hh = d[i];
-                d[i] = matrixA[i, i].Real;
+                d[i] = matrixA[i*order + i].Real;
-                matrixA[i, i] = new Complex32(hh, scale * (float)Math.Sqrt(h));
+                matrixA[i*order + i] = new Numerics.Complex32(hh, scale*(float) Math.Sqrt(h));
             }
 
             hh = d[0];
-            d[0] = matrixA[0, 0].Real;
+            d[0] = matrixA[0].Real;
-            matrixA[0, 0] = hh;
+            matrixA[0] = hh;
             e[0] = 0.0f;
         }
 
@@ -251,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricDiagonalize(Complex32[] dataEv, float[] d, float[] e, int order)
+        internal static void SymmetricDiagonalize(Numerics.Complex32[] dataEv, float[] d, float[] e, int order)
         {
             const int Maxiter = 1000;
 
@@ -351,8 +327,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                         {
                             throw new ArgumentException(Resources.ConvergenceFailed);
                         }
-                    }
+                    } while (Math.Abs(e[l]) > eps*tst1);
-                    while (Math.Abs(e[l]) > eps * tst1);
                 }
 
                 d[l] = d[l] + f;
@@ -398,7 +373,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// by Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricUntridiagonalize(Complex32[] dataEv, Complex32[,] matrixA, Complex32[] tau, int order)
+        internal static void SymmetricUntridiagonalize(Numerics.Complex32[] dataEv, Numerics.Complex32[] matrixA, Numerics.Complex32[] tau, int order)
         {
             for (var i = 0; i < order; i++)
             {
@@ -411,22 +386,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             // Recover and apply the Householder matrices.
             for (var i = 1; i < order; i++)
             {
-                var h = matrixA[i, i].Imaginary;
+                var h = matrixA[i*order + i].Imaginary;
                 if (h != 0)
                 {
                     for (var j = 0; j < order; j++)
                     {
-                        var s = Complex32.Zero;
+                        var s = Numerics.Complex32.Zero;
                         for (var k = 0; k < i; k++)
                         {
-                            s += dataEv[(j * order) + k] * matrixA[i, k];
+                            s += dataEv[(j*order) + k]*matrixA[k*order + i];
                         }
 
                         s = (s/h)/h;
 
                         for (var k = 0; k < i; k++)
                         {
-                            dataEv[(j * order) + k] -= s * matrixA[i, k].Conjugate();
+                            dataEv[(j*order) + k] -= s*matrixA[k*order + i].Conjugate();
                         }
                     }
                 }
@@ -443,17 +418,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutines in EISPACK.</remarks>
-        private static void NonsymmetricReduceToHessenberg(Complex32[] dataEv, Complex32[,] matrixH, int order)
+        internal static void NonsymmetricReduceToHessenberg(Numerics.Complex32[] dataEv, Numerics.Complex32[] matrixH, int order)
         {
-            var ort = new Complex32[order];
+            var ort = new Numerics.Complex32[order];
 
             for (var m = 1; m < order - 1; m++)
             {
                 // Scale column.
                 var scale = 0.0f;
+                var mm1O = (m - 1)*order;
                 for (var i = m; i < order; i++)
                 {
-                    scale += Math.Abs(matrixH[i, m - 1].Real) + Math.Abs(matrixH[i, m - 1].Imaginary);
+                    scale += Math.Abs(matrixH[mm1O + i].Real) + Math.Abs(matrixH[mm1O + i].Imaginary);
                 }
 
                 if (scale != 0.0f)
@@ -462,7 +438,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                     var h = 0.0f;
                     for (var i = order - 1; i >= m; i--)
                     {
-                        ort[i] = matrixH[i, m - 1] / scale;
+                        ort[i] = matrixH[mm1O + i]/scale;
                         h += ort[i].MagnitudeSquared;
                     }
 
@@ -476,43 +452,44 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                     else
                     {
                         ort[m] = g;
-                        matrixH[m, m - 1] = scale;
+                        matrixH[mm1O + m] = scale;
                     }
 
                     // Apply Householder similarity transformation
                     // H = (I-u*u'/h)*H*(I-u*u')/h)
                     for (var j = m; j < order; j++)
                     {
-                        var f = Complex32.Zero;
+                        var f = Numerics.Complex32.Zero;
+                        var jO = j*order;
                         for (var i = order - 1; i >= m; i--)
                         {
-                            f += ort[i].Conjugate() * matrixH[i, j];
+                            f += ort[i].Conjugate()*matrixH[jO + i];
                         }
 
                         f = f/h;
                         for (var i = m; i < order; i++)
                         {
-                            matrixH[i, j] -= f * ort[i];
+                            matrixH[jO + i] -= f*ort[i];
                         }
                     }
 
                     for (var i = 0; i < order; i++)
                     {
-                        var f = Complex32.Zero;
+                        var f = Numerics.Complex32.Zero;
                         for (var j = order - 1; j >= m; j--)
                         {
-                            f += ort[j] * matrixH[i, j];
+                            f += ort[j]*matrixH[j*order + i];
                         }
 
                         f = f/h;
                         for (var j = m; j < order; j++)
                         {
-                            matrixH[i, j] -= f * ort[j].Conjugate();
+                            matrixH[j*order + i] -= f*ort[j].Conjugate();
                         }
                     }
 
                     ort[m] = scale*ort[m];
-                    matrixH[m, m - 1] *= -g;
+                    matrixH[mm1O + m] *= -g;
                 }
             }
 
@@ -521,24 +498,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             {
                 for (var j = 0; j < order; j++)
                 {
-                    dataEv[(j * order) + i] = i == j ? Complex32.One : Complex32.Zero;
+                    dataEv[(j*order) + i] = i == j ? Numerics.Complex32.One : Numerics.Complex32.Zero;
                 }
             }
 
             for (var m = order - 2; m >= 1; m--)
             {
-                if (matrixH[m, m - 1] != Complex32.Zero && ort[m] != Complex32.Zero)
+                var mm1O = (m - 1)*order;
+                var mm1Om = mm1O + m;
+                if (matrixH[mm1Om] != Numerics.Complex32.Zero && ort[m] != Numerics.Complex32.Zero)
                 {
-                    var norm = (matrixH[m, m - 1].Real * ort[m].Real) + (matrixH[m, m - 1].Imaginary * ort[m].Imaginary);
+                    var norm = (matrixH[mm1Om].Real*ort[m].Real) + (matrixH[mm1Om].Imaginary*ort[m].Imaginary);
 
                     for (var i = m + 1; i < order; i++)
                     {
-                        ort[i] = matrixH[i, m - 1];
+                        ort[i] = matrixH[mm1O + i];
                     }
 
                     for (var j = m; j < order; j++)
                     {
-                        var g = Complex32.Zero;
+                        var g = Numerics.Complex32.Zero;
                         for (var i = m; i < order; i++)
                         {
                             g += ort[i].Conjugate()*dataEv[(j*order) + i];
@@ -557,18 +536,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             // Create real subdiagonal elements.
             for (var i = 1; i < order; i++)
             {
-                if (matrixH[i, i - 1].Imaginary != 0.0f)
+                var im1 = i - 1;
+                var im1O = im1*order;
+                var im1Oi = im1O + i;
+                var iO = i*order;
+                if (matrixH[im1Oi].Imaginary != 0.0f)
                 {
-                    var y = matrixH[i, i - 1] / matrixH[i, i - 1].Magnitude;
+                    var y = matrixH[im1Oi]/matrixH[im1Oi].Magnitude;
-                    matrixH[i, i - 1] = matrixH[i, i - 1].Magnitude;
+                    matrixH[im1Oi] = matrixH[im1Oi].Magnitude;
                     for (var j = i; j < order; j++)
                     {
-                        matrixH[i, j] *= y.Conjugate();
+                        matrixH[j*order + i] *= y.Conjugate();
                     }
 
                     for (var j = 0; j <= Math.Min(i + 1, order - 1); j++)
                     {
-                        matrixH[j, i] *= y;
+                        matrixH[iO + j] *= y;
                     }
 
                     for (var j = 0; j < order; j++)
@@ -590,14 +573,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, Complex32[] dataEv, Complex32[,] matrixH, int order)
+        internal static void NonsymmetricReduceHessenberToRealSchur(Numerics.Complex32[] vectorV, Numerics.Complex32[] dataEv, Numerics.Complex32[] matrixH, int order)
         {
             // Initialize
             var n = order - 1;
             var eps = (float) Precision.SingleMachinePrecision;
 
             float norm;
-            Complex32 x, y, z, exshift = Complex32.Zero;
+            Numerics.Complex32 x, y, z, exshift = Numerics.Complex32.Zero;
 
             // Outer loop over eigenvalue index
             var iter = 0;
@@ -607,8 +590,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 var l = n;
                 while (l > 0)
                 {
-                    var tst1 = Math.Abs(matrixH[l - 1, l - 1].Real) + Math.Abs(matrixH[l - 1, l - 1].Imaginary) + Math.Abs(matrixH[l, l].Real) + Math.Abs(matrixH[l, l].Imaginary);
+                    var lm1 = l - 1;
-                    if (Math.Abs(matrixH[l, l - 1].Real) < eps * tst1)
+                    var lm1O = lm1*order;
+                    var lO = l*order;
+                    var tst1 = Math.Abs(matrixH[lm1O + lm1].Real) + Math.Abs(matrixH[lm1O + lm1].Imaginary) + Math.Abs(matrixH[lO + l].Real) + Math.Abs(matrixH[lO + l].Imaginary);
+                    if (Math.Abs(matrixH[lm1O + l].Real) < eps*tst1)
                     {
                         break;
                     }
@@ -616,29 +602,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                     l--;
                 }
 
+                var nm1 = n - 1;
+                var nm1O = nm1*order;
+                var nO = n*order;
+                var nOn = nO + n;
                 // Check for convergence
                 // One root found
                 if (l == n)
                 {
-                    matrixH[n, n] += exshift;
+                    matrixH[nOn] += exshift;
-                    vectorV[n] = matrixH[n, n].ToComplex();
+                    vectorV[n] = matrixH[nOn];
                     n--;
                     iter = 0;
                 }
                 else
                 {
                     // Form shift
-                    Complex32 s;
+                    Numerics.Complex32 s;
                     if (iter != 10 && iter != 20)
                     {
-                        s = matrixH[n, n];
+                        s = matrixH[nOn];
-                        x = matrixH[n - 1, n] * matrixH[n, n - 1].Real;
+                        x = matrixH[nO + nm1]*matrixH[nm1O + n].Real;
 
                         if (x.Real != 0.0f || x.Imaginary != 0.0f)
                         {
-                            y = (matrixH[n - 1, n - 1] - s) / 2.0f;
+                            y = (matrixH[nm1O + nm1] - s)/2.0f;
                             z = ((y*y) + x).SquareRoot();
-                            if ((y.Real * z.Real) + (y.Imaginary * z.Imaginary) < 0.0f)
+                            if ((y.Real*z.Real) + (y.Imaginary*z.Imaginary) < 0.0)
                             {
                                 z *= -1.0f;
                             }
@@ -650,12 +640,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                     else
                     {
                         // Form exceptional shift
-                        s = Math.Abs(matrixH[n, n - 1].Real) + Math.Abs(matrixH[n - 1, n - 2].Real);
+                        s = Math.Abs(matrixH[nm1O + n].Real) + Math.Abs(matrixH[(n - 2)*order + nm1].Real);
                     }
 
                     for (var i = 0; i <= n; i++)
                     {
-                        matrixH[i, i] -= s;
+                        matrixH[i*order + i] -= s;
                     }
 
                     exshift += s;
@@ -664,61 +654,70 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                     // Reduce to triangle (rows)
                     for (var i = l + 1; i <= n; i++)
                     {
-                        s = matrixH[i, i - 1].Real;
+                        var im1 = i - 1;
-                        norm = SpecialFunctions.Hypotenuse(matrixH[i - 1, i - 1].Magnitude, s.Real);
+                        var im1O = im1*order;
-                        x = matrixH[i - 1, i - 1] / norm;
+                        var im1Oim1 = im1O + im1;
-                        vectorV[i - 1] = x.ToComplex();
+                        s = matrixH[im1O + i].Real;
-                        matrixH[i - 1, i - 1] = norm;
+                        norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real);
-                        matrixH[i, i - 1] = new Complex32(0.0f, s.Real / norm);
+                        x = matrixH[im1Oim1]/norm;
+                        vectorV[i - 1] = x;
+                        matrixH[im1Oim1] = norm;
+                        matrixH[im1O + i] = new Numerics.Complex32(0.0f, s.Real/norm);
 
                         for (var j = i; j < order; j++)
                         {
-                            y = matrixH[i - 1, j];
+                            var jO = j*order;
-                            z = matrixH[i, j];
+                            y = matrixH[jO + im1];
-                            matrixH[i - 1, j] = (x.Conjugate() * y) + (matrixH[i, i - 1].Imaginary * z);
+                            z = matrixH[jO + i];
-                            matrixH[i, j] = (x * z) - (matrixH[i, i - 1].Imaginary * y);
+                            matrixH[jO + im1] = (x.Conjugate()*y) + (matrixH[im1O + i].Imaginary*z);
+                            matrixH[jO + i] = (x*z) - (matrixH[im1O + i].Imaginary*y);
                         }
                     }
 
-                    s = matrixH[n, n];
+                    s = matrixH[nOn];
                     if (s.Imaginary != 0.0f)
                     {
-                        s /= matrixH[n, n].Magnitude;
+                        s /= matrixH[nOn].Magnitude;
-                        matrixH[n, n] = matrixH[n, n].Magnitude;
+                        matrixH[nOn] = matrixH[nOn].Magnitude;
 
                         for (var j = n + 1; j < order; j++)
                         {
-                            matrixH[n, j] *= s.Conjugate();
+                            matrixH[j*order + n] *= s.Conjugate();
                         }
                     }
 
                     // Inverse operation (columns).
                     for (var j = l + 1; j <= n; j++)
                     {
-                        x = (Complex32)vectorV[j - 1];
+                        x = vectorV[j - 1];
+                        var jO = j*order;
+                        var jm1 = j - 1;
+                        var jm1O = jm1*order;
+                        var jm1Oj = jm1O + j;
                         for (var i = 0; i <= j; i++)
                         {
-                            z = matrixH[i, j];
+                            var jm1Oi = jm1O + i;
+                            z = matrixH[jO + i];
                             if (i != j)
                             {
-                                y = matrixH[i, j - 1];
+                                y = matrixH[jm1Oi];
-                                matrixH[i, j - 1] = (x * y) + (matrixH[j, j - 1].Imaginary * z);
+                                matrixH[jm1Oi] = (x*y) + (matrixH[jm1O + j].Imaginary*z);
                             }
                             else
                             {
-                                y = matrixH[i, j - 1].Real;
+                                y = matrixH[jm1Oi].Real;
-                                matrixH[i, j - 1] = new Complex32((x.Real * y.Real) - (x.Imaginary * y.Imaginary) + (matrixH[j, j - 1].Imaginary * z.Real), matrixH[i, j - 1].Imaginary);
+                                matrixH[jm1Oi] = new Numerics.Complex32((x.Real*y.Real) - (x.Imaginary*y.Imaginary) + (matrixH[jm1O + j].Imaginary*z.Real), matrixH[jm1Oi].Imaginary);
                             }
 
-                            matrixH[i, j] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y);
+                            matrixH[jO + i] = (x.Conjugate()*z) - (matrixH[jm1O + j].Imaginary*y);
                         }
 
                         for (var i = 0; i < order; i++)
                         {
                             y = dataEv[((j - 1)*order) + i];
                             z = dataEv[(j*order) + i];
-                            dataEv[((j - 1) * order) + i] = (x * y) + (matrixH[j, j - 1].Imaginary * z);
+                            dataEv[jm1O + i] = (x*y) + (matrixH[jm1Oj].Imaginary*z);
-                            dataEv[(j * order) + i] = (x.Conjugate() * z) - (matrixH[j, j - 1].Imaginary * y);
+                            dataEv[jO + i] = (x.Conjugate()*z) - (matrixH[jm1Oj].Imaginary*y);
                         }
                     }
 
@@ -726,12 +725,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                     {
                         for (var i = 0; i <= n; i++)
                         {
-                            matrixH[i, n] *= s;
+                            matrixH[nO + i] *= s;
                         }
 
                         for (var i = 0; i < order; i++)
                         {
-                            dataEv[(n * order) + i] *= s;
+                            dataEv[nO + i] *= s;
                         }
                     }
                 }
@@ -744,7 +743,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             {
                 for (var j = i; j < order; j++)
                 {
-                    norm = Math.Max(norm, Math.Abs(matrixH[i, j].Real) + Math.Abs(matrixH[i, j].Imaginary));
+                    norm = Math.Max(norm, Math.Abs(matrixH[j*order + i].Real) + Math.Abs(matrixH[j*order + i].Imaginary));
                 }
             }
 
@@ -753,39 +752,41 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 return;
             }
 
-            if (norm == 0.0f)
+            if (norm == 0.0)
             {
                 return;
             }
 
             for (n = order - 1; n > 0; n--)
             {
-                x = (Complex32)vectorV[n];
+                var nO = n*order;
-                matrixH[n, n] = 1.0f;
+                var nOn = nO + n;
+                x = vectorV[n];
+                matrixH[nOn] = 1.0f;
 
                 for (var i = n - 1; i >= 0; i--)
                 {
                     z = 0.0f;
                     for (var j = i + 1; j <= n; j++)
                     {
-                        z += matrixH[i, j] * matrixH[j, n];
+                        z += matrixH[j*order + i]*matrixH[nO + j];
                     }
 
-                    y = x - (Complex32)vectorV[i];
+                    y = x - vectorV[i];
                     if (y.Real == 0.0f && y.Imaginary == 0.0f)
                     {
                         y = eps*norm;
                     }
 
-                    matrixH[i, n] = z / y;
+                    matrixH[nO + i] = z/y;
 
                     // Overflow control
-                    var tr = Math.Abs(matrixH[i, n].Real) + Math.Abs(matrixH[i, n].Imaginary);
+                    var tr = Math.Abs(matrixH[nO + i].Real) + Math.Abs(matrixH[nO + i].Imaginary);
                     if ((eps*tr)*tr > 1)
                     {
                         for (var j = i; j <= n; j++)
                         {
-                            matrixH[j, n] = matrixH[j, n] / tr;
+                            matrixH[nO + j] = matrixH[nO + j]/tr;
                         }
                     }
                 }
@@ -794,15 +795,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             // Back transformation to get eigenvectors of original matrix
             for (var j = order - 1; j > 0; j--)
             {
+                var jO = j*order;
                 for (var i = 0; i < order; i++)
                 {
-                    z = Complex32.Zero;
+                    z = Numerics.Complex32.Zero;
                     for (var k = 0; k <= j; k++)
                     {
-                        z += dataEv[(k * order) + i] * matrixH[k, j];
+                        z += dataEv[(k*order) + i]*matrixH[jO + k];
                     }
 
-                    dataEv[(j * order) + i] = z;
+                    dataEv[jO + i] = z;
                 }
             }
         }
@@ -812,7 +814,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// </summary>
         /// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
         /// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
-        public override void Solve(Matrix<Complex32> input, Matrix<Complex32> result)
+        public override void Solve(Matrix<Numerics.Complex32> input, Matrix<Numerics.Complex32> result)
         {
             // Check for proper arguments.
             if (input == null)
@@ -846,13 +848,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             if (IsSymmetric)
             {
                 var order = VectorEv.Count;
-                var tmp = new Complex32[order];
+                var tmp = new Numerics.Complex32[order];
 
                 for (var k = 0; k < order; k++)
                 {
                     for (var j = 0; j < order; j++)
                     {
-                        Complex32 value = 0.0f;
+                        Numerics.Complex32 value = 0.0f;
                         if (j < order)
                         {
                             for (var i = 0; i < order; i++)
@@ -868,7 +870,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 
                     for (var j = 0; j < order; j++)
                     {
-                        Complex32 value = 0.0f;
+                        Numerics.Complex32 value = 0.0f;
                         for (var i = 0; i < order; i++)
                         {
                             value += ((DenseMatrix) MatrixEv).Values[(i*order) + j]*tmp[i];
@@ -889,7 +891,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// </summary>
         /// <param name="input">The right hand side vector, <b>b</b>.</param>
         /// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
-        public override void Solve(Vector<Complex32> input, Vector<Complex32> result)
+        public override void Solve(Vector<Numerics.Complex32> input, Vector<Numerics.Complex32> result)
         {
             if (input == null)
             {
@@ -918,8 +920,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             {
                 // Symmetric case -> x = V * inv(λ) * VH * b;
                 var order = VectorEv.Count;
-                var tmp = new Complex32[order];
+                var tmp = new Numerics.Complex32[order];
-                Complex32 value;
+                Numerics.Complex32 value;
 
                 for (var j = 0; j < order; j++)
                 {
@@ -940,7 +942,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 for (var j = 0; j < order; j++)
                 {
                     value = 0;
-                    for (int i = 0; i < order; i++)
+                    for (var i = 0; i < order; i++)
                     {
                         value += ((DenseMatrix) MatrixEv).Values[(i*order) + j]*tmp[i];
                     }

src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs

@@ -80,9 +80,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
            
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                 }

src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -27,6 +27,7 @@
 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
+
 namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 {
     using System;
@@ -72,58 +73,23 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 
             var order = matrix.RowCount;
 
-            // Initialize matricies for eigenvalues and eigenvectors
+            // Initialize matrices for eigenvalues and eigenvectors
             MatrixEv = matrix.CreateMatrix(order, order);
             MatrixD = matrix.CreateMatrix(order, order);
             VectorEv = new LinearAlgebra.Complex.DenseVector(order);
 
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i);
                 }
             }
 
-            var d = new double[order];
+            Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) MatrixEv).Values,
-            var e = new double[order];
+                ((LinearAlgebra.Complex.DenseVector)VectorEv).Values, ((DenseMatrix)MatrixD).Values);
-
-            if (IsSymmetric)
-            {
-                matrix.CopyTo(MatrixEv);
-                d = MatrixEv.Row(order - 1).ToArray();
-
-                SymmetricTridiagonalize(((DenseMatrix)MatrixEv).Values, d, e, order);
-                SymmetricDiagonalize(((DenseMatrix)MatrixEv).Values, d, e, order);
-            }
-            else
-            {
-                var matrixH = matrix.ToArray();
-
-                NonsymmetricReduceToHessenberg(((DenseMatrix)MatrixEv).Values, matrixH, order);
-                NonsymmetricReduceHessenberToRealSchur(((DenseMatrix)MatrixEv).Values, matrixH, d, e, order);
-            }
-
-            for (var i = 0; i < order; i++)
-            {
-                MatrixD.At(i, i, d[i]);
-
-                if (e[i] > 0)
-                {
-                    MatrixD.At(i, i + 1, e[i]);
-                }
-                else if (e[i] < 0)
-                {
-                    MatrixD.At(i, i - 1, e[i]);
-                }
-            }
-
-            for (var i = 0; i < order; i++)
-            {
-                VectorEv[i] = new Complex(d[i], e[i]);
-            }
         }
 
         /// <summary>
@@ -137,7 +103,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricTridiagonalize(double[] a, double[] d, double[] e, int order)
+        internal static void SymmetricTridiagonalize(double[] a, double[] d, double[] e, int order)
         {
             // Householder reduction to tridiagonal form.
             for (var i = order - 1; i > 0; i--)
@@ -290,9 +256,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricDiagonalize(double[] a, double[] d, double[] e, int order)
+        internal static void SymmetricDiagonalize(double[] a, double[] d, double[] e, int order)
         {
-            const int Maxiter = 1000;
+            const int maxiter = 1000;
 
             for (var i = 1; i < order; i++)
             {
@@ -386,12 +352,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 
                         // Check for convergence. If too many iterations have been performed, 
                         // throw exception that Convergence Failed
-                        if (iter >= Maxiter)
+                        if (iter >= maxiter)
                         {
                             throw new ArgumentException(Resources.ConvergenceFailed);
                         }
-                    }
+                    } while (Math.Abs(e[l]) > eps*tst1);
-                    while (Math.Abs(e[l]) > eps * tst1);
                 }
 
                 d[l] = d[l] + f;
@@ -436,26 +401,28 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutines in EISPACK.</remarks>
-        private static void NonsymmetricReduceToHessenberg(double[] a, double[,] matrixH, int order)
+        internal static void NonsymmetricReduceToHessenberg(double[] a, double[] matrixH, int order)
         {
             var ort = new double[order];
-
+            var high = order - 1;
-            for (var m = 1; m < order - 1; m++)
+            for (var m = 1; m <= high - 1; m++)
             {
+                var mm1 = m - 1;
+                var mm1O = mm1*order;
                 // Scale column.
                 var scale = 0.0;
-                for (var i = m; i < order; i++)
+                for (var i = m; i <= high; i++)
                 {
-                    scale = scale + Math.Abs(matrixH[i, m - 1]);
+                    scale += Math.Abs(matrixH[mm1O + i]);
                 }
 
                 if (scale != 0.0)
                 {
                     // Compute Householder transformation.
                     var h = 0.0;
-                    for (var i = order - 1; i >= m; i--)
+                    for (var i = high; i >= m; i--)
                     {
-                        ort[i] = matrixH[i, m - 1] / scale;
+                        ort[i] = matrixH[mm1O + i]/scale;
                         h += ort[i]*ort[i];
                     }
 
@@ -472,36 +439,38 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                     // H = (I-u*u'/h)*H*(I-u*u')/h)
                     for (var j = m; j < order; j++)
                     {
+                        var jO = j*order;
                         var f = 0.0;
                         for (var i = order - 1; i >= m; i--)
                         {
-                            f += ort[i] * matrixH[i, j];
+                            f += ort[i]*matrixH[jO + i];
                         }
 
                         f = f/h;
-                        for (var i = m; i < order; i++)
+
+                        for (var i = m; i <= high; i++)
                         {
-                            matrixH[i, j] -= f * ort[i];
+                            matrixH[jO + i] -= f*ort[i];
                         }
                     }
 
-                    for (var i = 0; i < order; i++)
+                    for (var i = 0; i <= high; i++)
                     {
                         var f = 0.0;
-                        for (var j = order - 1; j >= m; j--)
+                        for (var j = high; j >= m; j--)
                         {
-                            f += ort[j] * matrixH[i, j];
+                            f += ort[j]*matrixH[j*order + i];
                         }
-
                         f = f/h;
-                        for (var j = m; j < order; j++)
+
+                        for (var j = m; j <= high; j++)
                         {
-                            matrixH[i, j] -= f * ort[j];
+                            matrixH[j*order + i] -= f*ort[j];
                         }
                     }
 
                     ort[m] = scale*ort[m];
-                    matrixH[m, m - 1] = scale * g;
+                    matrixH[mm1O + m] = scale*g;
                 }
             }
 
@@ -514,28 +483,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 }
             }
 
-            for (var m = order - 2; m >= 1; m--)
+            for (var m = high - 1; m >= 1; m--)
             {
-                if (matrixH[m, m - 1] != 0.0)
+                var mm1 = m - 1;
+                var mm1O = mm1*order;
+                var mm1Om = mm1O + m;
+                if (matrixH[mm1Om] != 0.0)
                 {
-                    for (var i = m + 1; i < order; i++)
+                    for (var i = m + 1; i <= high; i++)
                     {
-                        ort[i] = matrixH[i, m - 1];
+                        ort[i] = matrixH[mm1O + i];
                     }
 
-                    for (var j = m; j < order; j++)
+                    for (var j = m; j <= high; j++)
                     {
                         var g = 0.0;
-                        for (var i = m; i < order; i++)
+                        var jO = j*order;
+                        for (var i = m; i <= high; i++)
                         {
-                            g += ort[i] * a[(j * order) + i];
+                            g += ort[i]*a[jO + i];
                         }
 
                         // Double division avoids possible underflow
-                        g = (g / ort[m]) / matrixH[m, m - 1];
+                        g = (g/ort[m])/matrixH[mm1Om];
-                        for (var i = m; i < order; i++)
+
+                        for (var i = m; i <= high; i++)
                         {
-                            a[(j * order) + i] += g * ort[i];
+                            a[jO + i] += g*ort[i];
                         }
                     }
                 }
@@ -554,13 +528,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void NonsymmetricReduceHessenberToRealSchur(double[] a, double[,] matrixH, double[] d, double[] e, int order)
+        internal static void NonsymmetricReduceHessenberToRealSchur(double[] a, double[] matrixH, double[] d, double[] e, int order)
         {
             // Initialize
             var n = order - 1;
-            var eps = Precision.DoubleMachinePrecision;
+            var eps = Math.Pow(2.0, -52.0);
             var exshift = 0.0;
-            double p = 0, q = 0, r = 0, s = 0, z = 0, w, x, y;
+            double p = 0, q = 0, r = 0, s = 0, z = 0;
+            double w, x, y;
 
             // Store roots isolated by balanc and compute matrix norm
             var norm = 0.0;
@@ -568,7 +543,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             {
                 for (var j = Math.Max(i - 1, 0); j < order; j++)
                 {
-                    norm = norm + Math.Abs(matrixH[i, j]);
+                    norm = norm + Math.Abs(matrixH[j*order + i]);
                 }
             }
 
@@ -580,14 +555,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 var l = n;
                 while (l > 0)
                 {
-                    s = Math.Abs(matrixH[l - 1, l - 1]) + Math.Abs(matrixH[l, l]);
+                    var lm1 = l - 1;
+                    var lm1O = lm1*order;
+                    s = Math.Abs(matrixH[lm1O + lm1]) + Math.Abs(matrixH[l*order + l]);
 
                     if (s == 0.0)
                     {
                         s = norm;
                     }
 
-                    if (Math.Abs(matrixH[l, l - 1]) < eps * s)
+                    if (Math.Abs(matrixH[lm1O + l]) < eps*s)
                     {
                         break;
                     }
@@ -599,8 +576,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 // One root found
                 if (l == n)
                 {
-                    matrixH[n, n] = matrixH[n, n] + exshift;
+                    var index = n*order + n;
-                    d[n] = matrixH[n, n];
+                    matrixH[index] += exshift;
+                    d[n] = matrixH[index];
                     e[n] = 0.0;
                     n--;
                     iter = 0;
@@ -609,13 +587,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 }
                 else if (l == n - 1)
                 {
-                    w = matrixH[n, n - 1] * matrixH[n - 1, n];
+                    var nO = n*order;
-                    p = (matrixH[n - 1, n - 1] - matrixH[n, n]) / 2.0;
+                    var nm1 = n - 1;
+                    var nm1O = nm1*order;
+                    var nOn = nO + n;
+
+                    w = matrixH[nm1O + n]*matrixH[nO + nm1];
+                    p = (matrixH[nm1O + nm1] - matrixH[nOn])/2.0;
                     q = (p*p) + w;
                     z = Math.Sqrt(Math.Abs(q));
-                    matrixH[n, n] = matrixH[n, n] + exshift;
+
-                    matrixH[n - 1, n - 1] = matrixH[n - 1, n - 1] + exshift;
+                    matrixH[nOn] += exshift;
-                    x = matrixH[n, n];
+                    matrixH[nm1O + nm1] += exshift;
+                    x = matrixH[nOn];
 
                     // Real pair
                     if (q >= 0)
@@ -629,9 +613,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                             z = p - z;
                         }
 
-                        d[n - 1] = x + z;
+                        d[nm1] = x + z;
 
-                        d[n] = d[n - 1];
+                        d[n] = d[nm1];
                         if (z != 0.0)
                         {
                             d[n] = x - (w/z);
@@ -639,7 +623,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 
                         e[n - 1] = 0.0;
                         e[n] = 0.0;
-                        x = matrixH[n, n - 1];
+                        x = matrixH[nm1O + n];
                         s = Math.Abs(x) + Math.Abs(z);
                         p = x/s;
                         q = z/s;
@@ -650,25 +634,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                         // Row modification
                         for (var j = n - 1; j < order; j++)
                         {
-                            z = matrixH[n - 1, j];
+                            var jO = j*order;
-                            matrixH[n - 1, j] = (q * z) + (p * matrixH[n, j]);
+                            var jOn = jO + n;
-                            matrixH[n, j] = (q * matrixH[n, j]) - (p * z);
+                            z = matrixH[jO + nm1];
+                            matrixH[jO + nm1] = (q*z) + (p*matrixH[jOn]);
+                            matrixH[jOn] = (q*matrixH[jOn]) - (p*z);
                         }
 
                         // Column modification
                         for (var i = 0; i <= n; i++)
                         {
-                            z = matrixH[i, n - 1];
+                            var nOi = nO + i;
-                            matrixH[i, n - 1] = (q * z) + (p * matrixH[i, n]);
+                            z = matrixH[nm1O + i];
-                            matrixH[i, n] = (q * matrixH[i, n]) - (p * z);
+                            matrixH[nm1O + i] = (q*z) + (p*matrixH[nOi]);
+                            matrixH[nOi] = (q*matrixH[nOi]) - (p*z);
                         }
 
                         // Accumulate transformations
                         for (var i = 0; i < order; i++)
                         {
-                            z = a[((n - 1) * order) + i];
+                            var nOi = nO + i;
-                            a[((n - 1) * order) + i] = (q * z) + (p * a[(n * order) + i]);
+                            z = a[nm1O + i];
-                            a[(n * order) + i] = (q * a[(n * order) + i]) - (p * z);
+                            a[nm1O + i] = (q*z) + (p*a[nOi]);
+                            a[nOi] = (q*a[nOi]) - (p*z);
                         }
 
                         // Complex pair
@@ -688,14 +676,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 }
                 else
                 {
+                    var nO = n*order;
+                    var nm1 = n - 1;
+                    var nm1O = nm1*order;
+                    var nOn = nO + n;
+
                     // Form shift
-                    x = matrixH[n, n];
+                    x = matrixH[nOn];
                     y = 0.0;
                     w = 0.0;
                     if (l < n)
                     {
-                        y = matrixH[n - 1, n - 1];
+                        y = matrixH[nm1O + nm1];
-                        w = matrixH[n, n - 1] * matrixH[n - 1, n];
+                        w = matrixH[nm1O + n]*matrixH[nO + nm1];
                     }
 
                     // Wilkinson's original ad hoc shift
@@ -704,10 +697,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                         exshift += x;
                         for (var i = 0; i <= n; i++)
                         {
-                            matrixH[i, i] -= x;
+                            matrixH[i*order + i] -= x;
                         }
 
-                        s = Math.Abs(matrixH[n, n - 1]) + Math.Abs(matrixH[n - 1, n - 2]);
+                        s = Math.Abs(matrixH[nm1O + n]) + Math.Abs(matrixH[(n - 2)*order + nm1]);
                         x = y = 0.75*s;
                         w = (-0.4375)*s*s;
                     }
@@ -728,7 +721,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                             s = x - (w/(((y - x)/2.0) + s));
                             for (var i = 0; i <= n; i++)
                             {
-                                matrixH[i, i] -= s;
+                                matrixH[i*order + i] -= s;
                             }
 
                             exshift += s;
@@ -736,18 +729,28 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                         }
                     }
 
-                    iter = iter + 1; // (Could check iteration count here.)
+                    iter = iter + 1;
+                    if (iter >= 30*order)
+                    {
+                        throw new ArgumentException(Resources.ConvergenceFailed);
+                    }
 
                     // Look for two consecutive small sub-diagonal elements
                     var m = n - 2;
                     while (m >= l)
                     {
-                        z = matrixH[m, m];
+                        var mp1 = m + 1;
+                        var mm1 = m - 1;
+                        var mO = m*order;
+                        var mp1O = mp1*order;
+                        var mm1O = mm1*order;
+
+                        z = matrixH[mO + m];
                         r = x - z;
                         s = y - z;
-                        p = (((r * s) - w) / matrixH[m + 1, m]) + matrixH[m, m + 1];
+                        p = (((r*s) - w)/matrixH[mO + mp1]) + matrixH[mp1O + m];
-                        q = matrixH[m + 1, m + 1] - z - r - s;
+                        q = matrixH[mp1O + mp1] - z - r - s;
-                        r = matrixH[m + 2, m + 1];
+                        r = matrixH[mp1O + (m + 2)];
                         s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
                         p = p/s;
                         q = q/s;
@@ -758,7 +761,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                             break;
                         }
 
-                        if (Math.Abs(matrixH[m, m - 1]) * (Math.Abs(q) + Math.Abs(r)) < eps * (Math.Abs(p) * (Math.Abs(matrixH[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(matrixH[m + 1, m + 1]))))
+                        if (Math.Abs(matrixH[mm1O + m])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[mm1O + mm1]) + Math.Abs(z) + Math.Abs(matrixH[mp1O + mp1]))))
                         {
                             break;
                         }
@@ -766,40 +769,45 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                         m--;
                     }
 
-                    for (var i = m + 2; i <= n; i++)
+                    var mp2 = m + 2;
+                    for (var i = mp2; i <= n; i++)
                     {
-                        matrixH[i, i - 2] = 0.0;
+                        matrixH[(i - 2)*order + i] = 0.0;
-                        if (i > m + 2)
+                        if (i > mp2)
                         {
-                            matrixH[i, i - 3] = 0.0;
+                            matrixH[(i - 3)*order + i] = 0.0;
                         }
                     }
 
                     // Double QR step involving rows l:n and columns m:n
                     for (var k = m; k <= n - 1; k++)
                     {
-                        bool notlast = k != n - 1;
+                        var notlast = k != n - 1;
-
+                        var kO = k*order;
+                        var km1 = k - 1;
+                        var kp1 = k + 1;
+                        var kp2 = k + 2;
+                        var kp1O = kp1*order;
+                        var kp2O = kp2*order;
+                        var km1O = km1*order;
                         if (k != m)
                         {
-                            p = matrixH[k, k - 1];
+                            p = matrixH[km1O + k];
-                            q = matrixH[k + 1, k - 1];
+                            q = matrixH[km1O + kp1];
-                            r = notlast ? matrixH[k + 2, k - 1] : 0.0;
+                            r = notlast ? matrixH[km1O + kp2] : 0.0;
                             x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
-                            if (x != 0.0)
+                            if (x == 0.0)
                             {
+                                continue;
+                            }
+
                             p = p/x;
                             q = q/x;
                             r = r/x;
                         }
-                        }
-
-                        if (x == 0.0)
-                        {
-                            break;
-                        }
 
                         s = Math.Sqrt((p*p) + (q*q) + (r*r));
+
                         if (p < 0)
                         {
                             s = -s;
@@ -809,11 +817,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                         {
                             if (k != m)
                             {
-                                matrixH[k, k - 1] = (-s) * x;
+                                matrixH[km1O + k] = (-s)*x;
                             }
                             else if (l != m)
                             {
-                                matrixH[k, k - 1] = -matrixH[k, k - 1];
+                                matrixH[km1O + k] = -matrixH[km1O + k];
                             }
 
                             p = p + s;
@@ -826,46 +834,49 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                             // Row modification
                             for (var j = k; j < order; j++)
                             {
-                                p = matrixH[k, j] + (q * matrixH[k + 1, j]);
+                                var jO = j*order;
-
+                                var jOk = jO + k;
+                                var jOkp1 = jO + kp1;
+                                var jOkp2 = jO + kp2;
+                                p = matrixH[jOk] + (q*matrixH[jOkp1]);
                                 if (notlast)
                                 {
-                                    p = p + (r * matrixH[k + 2, j]);
+                                    p = p + (r*matrixH[jOkp2]);
-                                    matrixH[k + 2, j] = matrixH[k + 2, j] - (p * z);
+                                    matrixH[jOkp2] -= (p*z);
                                 }
 
-                                matrixH[k, j] = matrixH[k, j] - (p * x);
+                                matrixH[jOk] -= (p*x);
-                                matrixH[k + 1, j] = matrixH[k + 1, j] - (p * y);
+                                matrixH[jOkp1] -= (p*y);
                             }
 
                             // Column modification
                             for (var i = 0; i <= Math.Min(n, k + 3); i++)
                             {
-                                p = (x * matrixH[i, k]) + (y * matrixH[i, k + 1]);
+                                p = (x*matrixH[kO + i]) + (y*matrixH[kp1O + i]);
 
                                 if (notlast)
                                 {
-                                    p = p + (z * matrixH[i, k + 2]);
+                                    p = p + (z*matrixH[kp2O + i]);
-                                    matrixH[i, k + 2] = matrixH[i, k + 2] - (p * r);
+                                    matrixH[kp2O + i] -= (p*r);
                                 }
 
-                                matrixH[i, k] = matrixH[i, k] - p;
+                                matrixH[kO + i] -= p;
-                                matrixH[i, k + 1] = matrixH[i, k + 1] - (p * q);
+                                matrixH[kp1O + i] -= (p*q);
                             }
 
                             // Accumulate transformations
                             for (var i = 0; i < order; i++)
                             {
-                                p = (x * a[(k * order) + i]) + (y * a[((k + 1) * order) + i]);
+                                p = (x*a[kO + i]) + (y*a[kp1O + i]);
 
                                 if (notlast)
                                 {
-                                    p = p + (z * a[((k + 2) * order) + i]);
+                                    p = p + (z*a[kp2O + i]);
-                                    a[((k + 2) * order) + i] -= p * r;
+                                    a[kp2O + i] -= p*r;
                                 }
 
-                                a[(k * order) + i] -= p;
+                                a[kO + i] -= p;
-                                a[((k + 1) * order) + i] -= p * q;
+                                a[kp1O + i] -= p*q;
                             }
                         } // (s != 0)
                     } // k loop
@@ -880,23 +891,31 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 
             for (n = order - 1; n >= 0; n--)
             {
-                double t;
+                var nO = n*order;
+                var nm1 = n - 1;
+                var nm1O = nm1*order;
 
                 p = d[n];
                 q = e[n];
 
+
                 // Real vector
+                double t;
                 if (q == 0.0)
                 {
                     var l = n;
-                    matrixH[n, n] = 1.0;
+                    matrixH[nO + n] = 1.0;
                     for (var i = n - 1; i >= 0; i--)
                     {
-                        w = matrixH[i, i] - p;
+                        var ip1 = i + 1;
+                        var iO = i*order;
+                        var ip1O = ip1*order;
+
+                        w = matrixH[iO + i] - p;
                         r = 0.0;
                         for (var j = l; j <= n; j++)
                         {
-                            r = r + (matrixH[i, j] * matrixH[j, n]);
+                            r = r + (matrixH[j*order + i]*matrixH[nO + j]);
                         }
 
                         if (e[i] < 0.0)
@@ -911,39 +930,39 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                             {
                                 if (w != 0.0)
                                 {
-                                    matrixH[i, n] = (-r) / w;
+                                    matrixH[nO + i] = (-r)/w;
                                 }
                                 else
                                 {
-                                    matrixH[i, n] = (-r) / (eps * norm);
+                                    matrixH[nO + i] = (-r)/(eps*norm);
                                 }
 
                                 // Solve real equations
                             }
                             else
                             {
-                                x = matrixH[i, i + 1];
+                                x = matrixH[ip1O + i];
-                                y = matrixH[i + 1, i];
+                                y = matrixH[iO + ip1];
                                 q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]);
                                 t = ((x*s) - (z*r))/q;
-                                matrixH[i, n] = t;
+                                matrixH[nO + i] = t;
                                 if (Math.Abs(x) > Math.Abs(z))
                                 {
-                                    matrixH[i + 1, n] = (-r - (w * t)) / x;
+                                    matrixH[nO + ip1] = (-r - (w*t))/x;
                                 }
                                 else
                                 {
-                                    matrixH[i + 1, n] = (-s - (y * t)) / z;
+                                    matrixH[nO + ip1] = (-s - (y*t))/z;
                                 }
                             }
 
                             // Overflow control
-                            t = Math.Abs(matrixH[i, n]);
+                            t = Math.Abs(matrixH[nO + i]);
                             if ((eps*t)*t > 1)
                             {
                                 for (var j = i; j <= n; j++)
                                 {
-                                    matrixH[j, n] = matrixH[j, n] / t;
+                                    matrixH[nO + j] /= t;
                                 }
                             }
                         }
@@ -956,31 +975,36 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                     var l = n - 1;
 
                     // Last vector component imaginary so matrix is triangular
-                    if (Math.Abs(matrixH[n, n - 1]) > Math.Abs(matrixH[n - 1, n]))
+                    if (Math.Abs(matrixH[nm1O + n]) > Math.Abs(matrixH[nO + nm1]))
                     {
-                        matrixH[n - 1, n - 1] = q / matrixH[n, n - 1];
+                        matrixH[nm1O + nm1] = q/matrixH[nm1O + n];
-                        matrixH[n - 1, n] = (-(matrixH[n, n] - p)) / matrixH[n, n - 1];
+                        matrixH[nO + nm1] = (-(matrixH[nO + n] - p))/matrixH[nm1O + n];
                     }
                     else
                     {
-                        var res = Cdiv(0.0, -matrixH[n - 1, n], matrixH[n - 1, n - 1] - p, q);
+                        var res = Cdiv(0.0, -matrixH[nO + nm1], matrixH[nm1O + nm1] - p, q);
-                        matrixH[n - 1, n - 1] = res.Real;
+                        matrixH[nm1O + nm1] = res.Real;
-                        matrixH[n - 1, n] = res.Imaginary;
+                        matrixH[nO + nm1] = res.Imaginary;
                     }
 
-                    matrixH[n, n - 1] = 0.0;
+                    matrixH[nm1O + n] = 0.0;
-                    matrixH[n, n] = 1.0;
+                    matrixH[nO + n] = 1.0;
                     for (var i = n - 2; i >= 0; i--)
                     {
-                        double ra = 0.0;
+                        var ip1 = i + 1;
-                        double sa = 0.0;
+                        var iO = i*order;
+                        var ip1O = ip1*order;
+                        var ra = 0.0;
+                        var sa = 0.0;
                         for (var j = l; j <= n; j++)
                         {
-                            ra = ra + (matrixH[i, j] * matrixH[j, n - 1]);
+                            var jO = j*order;
-                            sa = sa + (matrixH[i, j] * matrixH[j, n]);
+                            var jOi = jO + i;
+                            ra = ra + (matrixH[jOi]*matrixH[nm1O + j]);
+                            sa = sa + (matrixH[jOi]*matrixH[nO + j]);
                         }
 
-                        w = matrixH[i, i] - p;
+                        w = matrixH[iO + i] - p;
 
                         if (e[i] < 0.0)
                         {
@@ -994,46 +1018,46 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                             if (e[i] == 0.0)
                             {
                                 var res = Cdiv(-ra, -sa, w, q);
-                                matrixH[i, n - 1] = res.Real;
+                                matrixH[nm1O + i] = res.Real;
-                                matrixH[i, n] = res.Imaginary;
+                                matrixH[nO + i] = res.Imaginary;
                             }
                             else
                             {
                                 // Solve complex equations
-                                x = matrixH[i, i + 1];
+                                x = matrixH[ip1O + i];
-                                y = matrixH[i + 1, i];
+                                y = matrixH[iO + ip1];
 
-                                double vr = ((d[i] - p) * (d[i] - p)) + (e[i] * e[i]) - (q * q);
+                                var vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q);
-                                double vi = (d[i] - p) * 2.0 * q;
+                                var vi = (d[i] - p)*2.0*q;
                                 if ((vr == 0.0) && (vi == 0.0))
                                 {
                                     vr = eps*norm*(Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z));
                                 }
 
                                 var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi);
-                                matrixH[i, n - 1] = res.Real;
+                                matrixH[nm1O + i] = res.Real;
-                                matrixH[i, n] = res.Imaginary;
+                                matrixH[nO + i] = res.Imaginary;
                                 if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q)))
                                 {
-                                    matrixH[i + 1, n - 1] = (-ra - (w * matrixH[i, n - 1]) + (q * matrixH[i, n])) / x;
+                                    matrixH[nm1O + ip1] = (-ra - (w*matrixH[nm1O + i]) + (q*matrixH[nO + i]))/x;
-                                    matrixH[i + 1, n] = (-sa - (w * matrixH[i, n]) - (q * matrixH[i, n - 1])) / x;
+                                    matrixH[nO + ip1] = (-sa - (w*matrixH[nO + i]) - (q*matrixH[nm1O + i]))/x;
                                 }
                                 else
                                 {
-                                    res = Cdiv(-r - (y * matrixH[i, n - 1]), -s - (y * matrixH[i, n]), z, q);
+                                    res = Cdiv(-r - (y*matrixH[nm1O + i]), -s - (y*matrixH[nO + i]), z, q);
-                                    matrixH[i + 1, n - 1] = res.Real;
+                                    matrixH[nm1O + ip1] = res.Real;
-                                    matrixH[i + 1, n] = res.Imaginary;
+                                    matrixH[nO + ip1] = res.Imaginary;
                                 }
                             }
 
                             // Overflow control
-                            t = Math.Max(Math.Abs(matrixH[i, n - 1]), Math.Abs(matrixH[i, n]));
+                            t = Math.Max(Math.Abs(matrixH[nm1O + i]), Math.Abs(matrixH[nO + i]));
                             if ((eps*t)*t > 1)
                             {
                                 for (var j = i; j <= n; j++)
                                 {
-                                    matrixH[j, n - 1] = matrixH[j, n - 1] / t;
+                                    matrixH[nm1O + j] /= t;
-                                    matrixH[j, n] = matrixH[j, n] / t;
+                                    matrixH[nO + j] /= t;
                                 }
                             }
                         }
@@ -1044,15 +1068,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             // Back transformation to get eigenvectors of original matrix
             for (var j = order - 1; j >= 0; j--)
             {
+                var jO = j*order;
                 for (var i = 0; i < order; i++)
                 {
                     z = 0.0;
                     for (var k = 0; k <= j; k++)
                     {
-                        z = z + (a[(k * order) + i] * matrixH[k, j]);
+                        z = z + (a[k*order + i]*matrixH[jO + k]);
                     }
 
-                    a[(j * order) + i] = z;
+                    a[jO + i] = z;
                 }
             }
         }
@@ -1065,14 +1090,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// <param name="yreal">Real part of Y</param>
         /// <param name="yimag">Imaginary part of Y</param>
         /// <returns>Division result as a <see cref="Complex"/> number.</returns>
-        private static Complex Cdiv(double xreal, double ximag, double yreal, double yimag)
+        static System.Numerics.Complex Cdiv(double xreal, double ximag, double yreal, double yimag)
         {
             if (Math.Abs(yimag) < Math.Abs(yreal))
             {
-                return new Complex((xreal + (ximag * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal))), (ximag - (xreal * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal))));
+                return new System.Numerics.Complex((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))));
             }
 
-            return new Complex((ximag + (xreal * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag))), (-xreal + (ximag * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag))));
+            return new System.Numerics.Complex((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))));
         }
 
         /// <summary>
@@ -1208,7 +1233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 for (var j = 0; j < order; j++)
                 {
                     value = 0;
-                    for (int i = 0; i < order; i++)
+                    for (var i = 0; i < order; i++)
                     {
                         value += ((DenseMatrix) MatrixEv).Values[(i*order) + j]*tmp[i];
                     }

src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs

@@ -79,9 +79,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
            
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i);
                 }

src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -27,12 +27,12 @@
 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
+
 namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 {
     using System;
     using System.Numerics;
     using Generic;
-    using Numerics;
     using Properties;
 
     /// <summary>
@@ -73,58 +73,23 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 
             var order = matrix.RowCount;
 
-            // Initialize matricies for eigenvalues and eigenvectors
+            // Initialize matrices for eigenvalues and eigenvectors
             MatrixEv = matrix.CreateMatrix(order, order);
             MatrixD = matrix.CreateMatrix(order, order);
             VectorEv = new LinearAlgebra.Complex.DenseVector(order);
 
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i);
                 }
             }
 
-            var d = new float[order];
+            Control.LinearAlgebraProvider.EigenDecomp(IsSymmetric, order, matrix.Values, ((DenseMatrix) MatrixEv).Values,
-            var e = new float[order];
+                ((LinearAlgebra.Complex.DenseVector) VectorEv).Values, ((DenseMatrix) MatrixD).Values);
-
-            if (IsSymmetric)
-            {
-                matrix.CopyTo(MatrixEv);
-                d = MatrixEv.Row(order - 1).ToArray();
-
-                SymmetricTridiagonalize(((DenseMatrix)MatrixEv).Values, d, e, order);
-                SymmetricDiagonalize(((DenseMatrix)MatrixEv).Values, d, e, order);
-            }
-            else
-            {
-                var matrixH = matrix.ToArray();
-
-                NonsymmetricReduceToHessenberg(((DenseMatrix)MatrixEv).Values, matrixH, order);
-                NonsymmetricReduceHessenberToRealSchur(((DenseMatrix)MatrixEv).Values, matrixH, d, e, order);
-            }
-
-            for (var i = 0; i < order; i++)
-            {
-                MatrixD.At(i, i, d[i]);
-
-                if (e[i] > 0)
-                {
-                    MatrixD.At(i, i + 1, e[i]);
-                }
-                else if (e[i] < 0)
-                {
-                    MatrixD.At(i, i - 1, e[i]);
-                }
-            }
-
-            for (var i = 0; i < order; i++)
-            {
-                VectorEv[i] = new Complex(d[i], e[i]);
-            }
         }
 
         /// <summary>
@@ -138,7 +103,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricTridiagonalize(float[] a, float[] d, float[] e, int order)
+        internal static void SymmetricTridiagonalize(float[] a, float[] d, float[] e, int order)
         {
             // Householder reduction to tridiagonal form.
             for (var i = order - 1; i > 0; i--)
@@ -291,7 +256,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
         /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void SymmetricDiagonalize(float[] a, float[] d, float[] e, int order)
+        internal static void SymmetricDiagonalize(float[] a, float[] d, float[] e, int order)
         {
             const int Maxiter = 1000;
 
@@ -304,7 +269,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 
             var f = 0.0f;
             var tst1 = 0.0f;
-            var eps = Precision.DoubleMachinePrecision;
+            var eps = Precision.SingleMachinePrecision;
             for (var l = 0; l < order; l++)
             {
                 // Find small subdiagonal element
@@ -391,8 +356,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                         {
                             throw new ArgumentException(Resources.ConvergenceFailed);
                         }
-                    }
+                    } while (Math.Abs(e[l]) > eps*tst1);
-                    while (Math.Abs(e[l]) > eps * tst1);
                 }
 
                 d[l] = d[l] + f;
@@ -437,26 +401,28 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutines in EISPACK.</remarks>
-        private static void NonsymmetricReduceToHessenberg(float[] a, float[,] matrixH, int order)
+        internal static void NonsymmetricReduceToHessenberg(float[] a, float[] matrixH, int order)
         {
             var ort = new float[order];
-
+            var high = order - 1;
-            for (var m = 1; m < order - 1; m++)
+            for (var m = 1; m <= high - 1; m++)
             {
+                var mm1 = m - 1;
+                var mm1O = mm1*order;
                 // Scale column.
                 var scale = 0.0f;
-                for (var i = m; i < order; i++)
+                for (var i = m; i <= high; i++)
                 {
-                    scale = scale + Math.Abs(matrixH[i, m - 1]);
+                    scale += Math.Abs(matrixH[mm1O + i]);
                 }
 
                 if (scale != 0.0f)
                 {
                     // Compute Householder transformation.
                     var h = 0.0f;
-                    for (var i = order - 1; i >= m; i--)
+                    for (var i = high; i >= m; i--)
                     {
-                        ort[i] = matrixH[i, m - 1] / scale;
+                        ort[i] = matrixH[mm1O + i]/scale;
                         h += ort[i]*ort[i];
                     }
 
@@ -473,36 +439,38 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                     // H = (I-u*u'/h)*H*(I-u*u')/h)
                     for (var j = m; j < order; j++)
                     {
+                        var jO = j*order;
                         var f = 0.0f;
                         for (var i = order - 1; i >= m; i--)
                         {
-                            f += ort[i] * matrixH[i, j];
+                            f += ort[i]*matrixH[jO + i];
                         }
 
                         f = f/h;
-                        for (var i = m; i < order; i++)
+
+                        for (var i = m; i <= high; i++)
                         {
-                            matrixH[i, j] -= f * ort[i];
+                            matrixH[jO + i] -= f*ort[i];
                         }
                     }
 
-                    for (var i = 0; i < order; i++)
+                    for (var i = 0; i <= high; i++)
                     {
                         var f = 0.0f;
-                        for (var j = order - 1; j >= m; j--)
+                        for (var j = high; j >= m; j--)
                         {
-                            f += ort[j] * matrixH[i, j];
+                            f += ort[j]*matrixH[j*order + i];
                         }
-
                         f = f/h;
-                        for (var j = m; j < order; j++)
+
+                        for (var j = m; j <= high; j++)
                         {
-                            matrixH[i, j] -= f * ort[j];
+                            matrixH[j*order + i] -= f*ort[j];
                         }
                     }
 
                     ort[m] = scale*ort[m];
-                    matrixH[m, m - 1] = scale * g;
+                    matrixH[mm1O + m] = scale*g;
                 }
             }
 
@@ -515,28 +483,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 }
             }
 
-            for (var m = order - 2; m >= 1; m--)
+            for (var m = high - 1; m >= 1; m--)
             {
-                if (matrixH[m, m - 1] != 0.0f)
+                var mm1 = m - 1;
+                var mm1O = mm1*order;
+                var mm1Om = mm1O + m;
+                if (matrixH[mm1Om] != 0.0)
                 {
-                    for (var i = m + 1; i < order; i++)
+                    for (var i = m + 1; i <= high; i++)
                     {
-                        ort[i] = matrixH[i, m - 1];
+                        ort[i] = matrixH[mm1O + i];
                     }
 
-                    for (var j = m; j < order; j++)
+                    for (var j = m; j <= high; j++)
                     {
                         var g = 0.0f;
-                        for (var i = m; i < order; i++)
+                        var jO = j*order;
+                        for (var i = m; i <= high; i++)
                         {
-                            g += ort[i] * a[(j * order) + i];
+                            g += ort[i]*a[jO + i];
                         }
 
                         // Double division avoids possible underflow
-                        g = (g / ort[m]) / matrixH[m, m - 1];
+                        g = (g/ort[m])/matrixH[mm1Om];
-                        for (var i = m; i < order; i++)
+
+                        for (var i = m; i <= high; i++)
                         {
-                            a[(j * order) + i] += g * ort[i];
+                            a[jO + i] += g*ort[i];
                         }
                     }
                 }
@@ -555,13 +528,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// by Martin and Wilkinson, Handbook for Auto. Comp.,
         /// Vol.ii-Linear Algebra, and the corresponding
         /// Fortran subroutine in EISPACK.</remarks>
-        private static void NonsymmetricReduceHessenberToRealSchur(float[] a, float[,] matrixH, float[] d, float[] e, int order)
+        internal static void NonsymmetricReduceHessenberToRealSchur(float[] a, float[] matrixH, float[] d, float[] e, int order)
         {
             // Initialize
             var n = order - 1;
             var eps = (float) Precision.SingleMachinePrecision;
             var exshift = 0.0f;
-            float p = 0, q = 0, r = 0, s = 0, z = 0, w, x, y;
+            float p = 0, q = 0, r = 0, s = 0, z = 0;
+            float w, x, y;
 
             // Store roots isolated by balanc and compute matrix norm
             var norm = 0.0f;
@@ -569,7 +543,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             {
                 for (var j = Math.Max(i - 1, 0); j < order; j++)
                 {
-                    norm = norm + Math.Abs(matrixH[i, j]);
+                    norm = norm + Math.Abs(matrixH[j*order + i]);
                 }
             }
 
@@ -581,14 +555,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 var l = n;
                 while (l > 0)
                 {
-                    s = Math.Abs(matrixH[l - 1, l - 1]) + Math.Abs(matrixH[l, l]);
+                    var lm1 = l - 1;
+                    var lm1O = lm1*order;
+                    s = Math.Abs(matrixH[lm1O + lm1]) + Math.Abs(matrixH[l*order + l]);
 
-                    if (s == 0.0f)
+                    if (s == 0.0)
                     {
                         s = norm;
                     }
 
-                    if (Math.Abs(matrixH[l, l - 1]) < eps * s)
+                    if (Math.Abs(matrixH[lm1O + l]) < eps*s)
                     {
                         break;
                     }
@@ -600,8 +576,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 // One root found
                 if (l == n)
                 {
-                    matrixH[n, n] = matrixH[n, n] + exshift;
+                    var index = n*order + n;
-                    d[n] = matrixH[n, n];
+                    matrixH[index] += exshift;
+                    d[n] = matrixH[index];
                     e[n] = 0.0f;
                     n--;
                     iter = 0;
@@ -610,13 +587,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 }
                 else if (l == n - 1)
                 {
-                    w = matrixH[n, n - 1] * matrixH[n - 1, n];
+                    var nO = n*order;
-                    p = (matrixH[n - 1, n - 1] - matrixH[n, n]) / 2.0f;
+                    var nm1 = n - 1;
+                    var nm1O = nm1*order;
+                    var nOn = nO + n;
+
+                    w = matrixH[nm1O + n]*matrixH[nO + nm1];
+                    p = (matrixH[nm1O + nm1] - matrixH[nOn])/2.0f;
                     q = (p*p) + w;
                     z = (float) Math.Sqrt(Math.Abs(q));
-                    matrixH[n, n] = matrixH[n, n] + exshift;
+
-                    matrixH[n - 1, n - 1] = matrixH[n - 1, n - 1] + exshift;
+                    matrixH[nOn] += exshift;
-                    x = matrixH[n, n];
+                    matrixH[nm1O + nm1] += exshift;
+                    x = matrixH[nOn];
 
                     // Real pair
                     if (q >= 0)
@@ -630,17 +613,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                             z = p - z;
                         }
 
-                        d[n - 1] = x + z;
+                        d[nm1] = x + z;
 
-                        d[n] = d[n - 1];
+                        d[n] = d[nm1];
-                        if (z != 0.0f)
+                        if (z != 0.0)
                         {
                             d[n] = x - (w/z);
                         }
 
                         e[n - 1] = 0.0f;
                         e[n] = 0.0f;
-                        x = matrixH[n, n - 1];
+                        x = matrixH[nm1O + n];
                         s = Math.Abs(x) + Math.Abs(z);
                         p = x/s;
                         q = z/s;
@@ -651,25 +634,29 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                         // Row modification
                         for (var j = n - 1; j < order; j++)
                         {
-                            z = matrixH[n - 1, j];
+                            var jO = j*order;
-                            matrixH[n - 1, j] = (q * z) + (p * matrixH[n, j]);
+                            var jOn = jO + n;
-                            matrixH[n, j] = (q * matrixH[n, j]) - (p * z);
+                            z = matrixH[jO + nm1];
+                            matrixH[jO + nm1] = (q*z) + (p*matrixH[jOn]);
+                            matrixH[jOn] = (q*matrixH[jOn]) - (p*z);
                         }
 
                         // Column modification
                         for (var i = 0; i <= n; i++)
                         {
-                            z = matrixH[i, n - 1];
+                            var nOi = nO + i;
-                            matrixH[i, n - 1] = (q * z) + (p * matrixH[i, n]);
+                            z = matrixH[nm1O + i];
-                            matrixH[i, n] = (q * matrixH[i, n]) - (p * z);
+                            matrixH[nm1O + i] = (q*z) + (p*matrixH[nOi]);
+                            matrixH[nOi] = (q*matrixH[nOi]) - (p*z);
                         }
 
                         // Accumulate transformations
                         for (var i = 0; i < order; i++)
                         {
-                            z = a[((n - 1) * order) + i];
+                            var nOi = nO + i;
-                            a[((n - 1) * order) + i] = (q * z) + (p * a[(n * order) + i]);
+                            z = a[nm1O + i];
-                            a[(n * order) + i] = (q * a[(n * order) + i]) - (p * z);
+                            a[nm1O + i] = (q*z) + (p*a[nOi]);
+                            a[nOi] = (q*a[nOi]) - (p*z);
                         }
 
                         // Complex pair
@@ -689,14 +676,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 }
                 else
                 {
+                    var nO = n*order;
+                    var nm1 = n - 1;
+                    var nm1O = nm1*order;
+                    var nOn = nO + n;
+
                     // Form shift
-                    x = matrixH[n, n];
+                    x = matrixH[nOn];
                     y = 0.0f;
                     w = 0.0f;
                     if (l < n)
                     {
-                        y = matrixH[n - 1, n - 1];
+                        y = matrixH[nm1O + nm1];
-                        w = matrixH[n, n - 1] * matrixH[n - 1, n];
+                        w = matrixH[nm1O + n]*matrixH[nO + nm1];
                     }
 
                     // Wilkinson's original ad hoc shift
@@ -705,10 +697,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                         exshift += x;
                         for (var i = 0; i <= n; i++)
                         {
-                            matrixH[i, i] -= x;
+                            matrixH[i*order + i] -= x;
                         }
 
-                        s = Math.Abs(matrixH[n, n - 1]) + Math.Abs(matrixH[n - 1, n - 2]);
+                        s = Math.Abs(matrixH[nm1O + n]) + Math.Abs(matrixH[(n - 2)*order + nm1]);
                         x = y = 0.75f*s;
                         w = (-0.4375f)*s*s;
                     }
@@ -729,7 +721,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                             s = x - (w/(((y - x)/2.0f) + s));
                             for (var i = 0; i <= n; i++)
                             {
-                                matrixH[i, i] -= s;
+                                matrixH[i*order + i] -= s;
                             }
 
                             exshift += s;
@@ -737,18 +729,28 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                         }
                     }
 
-                    iter = iter + 1; // (Could check iteration count here.)
+                    iter = iter + 1;
+                    if (iter >= 30*order)
+                    {
+                        throw new ArgumentException(Resources.ConvergenceFailed);
+                    }
 
                     // Look for two consecutive small sub-diagonal elements
                     var m = n - 2;
                     while (m >= l)
                     {
-                        z = matrixH[m, m];
+                        var mp1 = m + 1;
+                        var mm1 = m - 1;
+                        var mO = m*order;
+                        var mp1O = mp1*order;
+                        var mm1O = mm1*order;
+
+                        z = matrixH[mO + m];
                         r = x - z;
                         s = y - z;
-                        p = (((r * s) - w) / matrixH[m + 1, m]) + matrixH[m, m + 1];
+                        p = (((r*s) - w)/matrixH[mO + mp1]) + matrixH[mp1O + m];
-                        q = matrixH[m + 1, m + 1] - z - r - s;
+                        q = matrixH[mp1O + mp1] - z - r - s;
-                        r = matrixH[m + 2, m + 1];
+                        r = matrixH[mp1O + (m + 2)];
                         s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
                         p = p/s;
                         q = q/s;
@@ -759,7 +761,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                             break;
                         }
 
-                        if (Math.Abs(matrixH[m, m - 1]) * (Math.Abs(q) + Math.Abs(r)) < eps * (Math.Abs(p) * (Math.Abs(matrixH[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(matrixH[m + 1, m + 1]))))
+                        if (Math.Abs(matrixH[mm1O + m])*(Math.Abs(q) + Math.Abs(r)) < eps*(Math.Abs(p)*(Math.Abs(matrixH[mm1O + mm1]) + Math.Abs(z) + Math.Abs(matrixH[mp1O + mp1]))))
                         {
                             break;
                         }
@@ -767,40 +769,45 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                         m--;
                     }
 
-                    for (var i = m + 2; i <= n; i++)
+                    var mp2 = m + 2;
+                    for (var i = mp2; i <= n; i++)
                     {
-                        matrixH[i, i - 2] = 0.0f;
+                        matrixH[(i - 2)*order + i] = 0.0f;
-                        if (i > m + 2)
+                        if (i > mp2)
                         {
-                            matrixH[i, i - 3] = 0.0f;
+                            matrixH[(i - 3)*order + i] = 0.0f;
                         }
                     }
 
                     // Double QR step involving rows l:n and columns m:n
                     for (var k = m; k <= n - 1; k++)
                     {
-                        bool notlast = k != n - 1;
+                        var notlast = k != n - 1;
-
+                        var kO = k*order;
+                        var km1 = k - 1;
+                        var kp1 = k + 1;
+                        var kp2 = k + 2;
+                        var kp1O = kp1*order;
+                        var kp2O = kp2*order;
+                        var km1O = km1*order;
                         if (k != m)
                         {
-                            p = matrixH[k, k - 1];
+                            p = matrixH[km1O + k];
-                            q = matrixH[k + 1, k - 1];
+                            q = matrixH[km1O + kp1];
-                            r = notlast ? matrixH[k + 2, k - 1] : 0.0f;
+                            r = notlast ? matrixH[km1O + kp2] : 0.0f;
                             x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
-                            if (x != 0.0f)
+                            if (x == 0.0f)
                             {
+                                continue;
+                            }
+
                             p = p/x;
                             q = q/x;
                             r = r/x;
                         }
-                        }
-
-                        if (x == 0.0f)
-                        {
-                            break;
-                        }
 
                         s = (float) Math.Sqrt((p*p) + (q*q) + (r*r));
+
                         if (p < 0)
                         {
                             s = -s;
@@ -810,11 +817,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                         {
                             if (k != m)
                             {
-                                matrixH[k, k - 1] = (-s) * x;
+                                matrixH[km1O + k] = (-s)*x;
                             }
                             else if (l != m)
                             {
-                                matrixH[k, k - 1] = -matrixH[k, k - 1];
+                                matrixH[km1O + k] = -matrixH[km1O + k];
                             }
 
                             p = p + s;
@@ -827,46 +834,49 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                             // Row modification
                             for (var j = k; j < order; j++)
                             {
-                                p = matrixH[k, j] + (q * matrixH[k + 1, j]);
+                                var jO = j*order;
-
+                                var jOk = jO + k;
+                                var jOkp1 = jO + kp1;
+                                var jOkp2 = jO + kp2;
+                                p = matrixH[jOk] + (q*matrixH[jOkp1]);
                                 if (notlast)
                                 {
-                                    p = p + (r * matrixH[k + 2, j]);
+                                    p = p + (r*matrixH[jOkp2]);
-                                    matrixH[k + 2, j] = matrixH[k + 2, j] - (p * z);
+                                    matrixH[jOkp2] -= (p*z);
                                 }
 
-                                matrixH[k, j] = matrixH[k, j] - (p * x);
+                                matrixH[jOk] -= (p*x);
-                                matrixH[k + 1, j] = matrixH[k + 1, j] - (p * y);
+                                matrixH[jOkp1] -= (p*y);
                             }
 
                             // Column modification
                             for (var i = 0; i <= Math.Min(n, k + 3); i++)
                             {
-                                p = (x * matrixH[i, k]) + (y * matrixH[i, k + 1]);
+                                p = (x*matrixH[kO + i]) + (y*matrixH[kp1O + i]);
 
                                 if (notlast)
                                 {
-                                    p = p + (z * matrixH[i, k + 2]);
+                                    p = p + (z*matrixH[kp2O + i]);
-                                    matrixH[i, k + 2] = matrixH[i, k + 2] - (p * r);
+                                    matrixH[kp2O + i] -= (p*r);
                                 }
 
-                                matrixH[i, k] = matrixH[i, k] - p;
+                                matrixH[kO + i] -= p;
-                                matrixH[i, k + 1] = matrixH[i, k + 1] - (p * q);
+                                matrixH[kp1O + i] -= (p*q);
                             }
 
                             // Accumulate transformations
                             for (var i = 0; i < order; i++)
                             {
-                                p = (x * a[(k * order) + i]) + (y * a[((k + 1) * order) + i]);
+                                p = (x*a[kO + i]) + (y*a[kp1O + i]);
 
                                 if (notlast)
                                 {
-                                    p = p + (z * a[((k + 2) * order) + i]);
+                                    p = p + (z*a[kp2O + i]);
-                                    a[((k + 2) * order) + i] -= p * r;
+                                    a[kp2O + i] -= p*r;
                                 }
 
-                                a[(k * order) + i] -= p;
+                                a[kO + i] -= p;
-                                a[((k + 1) * order) + i] -= p * q;
+                                a[kp1O + i] -= p*q;
                             }
                         } // (s != 0)
                     } // k loop
@@ -881,26 +891,34 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 
             for (n = order - 1; n >= 0; n--)
             {
-                float t;
+                var nO = n*order;
+                var nm1 = n - 1;
+                var nm1O = nm1*order;
 
                 p = d[n];
                 q = e[n];
 
+
                 // Real vector
+                float t;
                 if (q == 0.0f)
                 {
                     var l = n;
-                    matrixH[n, n] = 1.0f;
+                    matrixH[nO + n] = 1.0f;
                     for (var i = n - 1; i >= 0; i--)
                     {
-                        w = matrixH[i, i] - p;
+                        var ip1 = i + 1;
+                        var iO = i*order;
+                        var ip1O = ip1*order;
+
+                        w = matrixH[iO + i] - p;
                         r = 0.0f;
                         for (var j = l; j <= n; j++)
                         {
-                            r = r + (matrixH[i, j] * matrixH[j, n]);
+                            r = r + (matrixH[j*order + i]*matrixH[nO + j]);
                         }
 
-                        if (e[i] < 0.0f)
+                        if (e[i] < 0.0)
                         {
                             z = w;
                             s = r;
@@ -912,39 +930,39 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                             {
                                 if (w != 0.0f)
                                 {
-                                    matrixH[i, n] = (-r) / w;
+                                    matrixH[nO + i] = (-r)/w;
                                 }
                                 else
                                 {
-                                    matrixH[i, n] = (-r) / (eps * norm);
+                                    matrixH[nO + i] = (-r)/(eps*norm);
                                 }
 
                                 // Solve real equations
                             }
                             else
                             {
-                                x = matrixH[i, i + 1];
+                                x = matrixH[ip1O + i];
-                                y = matrixH[i + 1, i];
+                                y = matrixH[iO + ip1];
                                 q = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]);
                                 t = ((x*s) - (z*r))/q;
-                                matrixH[i, n] = t;
+                                matrixH[nO + i] = t;
                                 if (Math.Abs(x) > Math.Abs(z))
                                 {
-                                    matrixH[i + 1, n] = (-r - (w * t)) / x;
+                                    matrixH[nO + ip1] = (-r - (w*t))/x;
                                 }
                                 else
                                 {
-                                    matrixH[i + 1, n] = (-s - (y * t)) / z;
+                                    matrixH[nO + ip1] = (-s - (y*t))/z;
                                 }
                             }
 
                             // Overflow control
-                            t = Math.Abs(matrixH[i, n]);
+                            t = Math.Abs(matrixH[nO + i]);
                             if ((eps*t)*t > 1)
                             {
                                 for (var j = i; j <= n; j++)
                                 {
-                                    matrixH[j, n] = matrixH[j, n] / t;
+                                    matrixH[nO + j] /= t;
                                 }
                             }
                         }
@@ -957,33 +975,38 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                     var l = n - 1;
 
                     // Last vector component imaginary so matrix is triangular
-                    if (Math.Abs(matrixH[n, n - 1]) > Math.Abs(matrixH[n - 1, n]))
+                    if (Math.Abs(matrixH[nm1O + n]) > Math.Abs(matrixH[nO + nm1]))
                     {
-                        matrixH[n - 1, n - 1] = q / matrixH[n, n - 1];
+                        matrixH[nm1O + nm1] = q/matrixH[nm1O + n];
-                        matrixH[n - 1, n] = (-(matrixH[n, n] - p)) / matrixH[n, n - 1];
+                        matrixH[nO + nm1] = (-(matrixH[nO + n] - p))/matrixH[nm1O + n];
                     }
                     else
                     {
-                        var res = Cdiv(0.0f, -matrixH[n - 1, n], matrixH[n - 1, n - 1] - p, q);
+                        var res = Cdiv(0.0f, -matrixH[nO + nm1], matrixH[nm1O + nm1] - p, q);
-                        matrixH[n - 1, n - 1] = res.Real;
+                        matrixH[nm1O + nm1] = res.Real;
-                        matrixH[n - 1, n] = res.Imaginary;
+                        matrixH[nO + nm1] = res.Imaginary;
                     }
 
-                    matrixH[n, n - 1] = 0.0f;
+                    matrixH[nm1O + n] = 0.0f;
-                    matrixH[n, n] = 1.0f;
+                    matrixH[nO + n] = 1.0f;
                     for (var i = n - 2; i >= 0; i--)
                     {
-                        float ra = 0.0f;
+                        var ip1 = i + 1;
-                        float sa = 0.0f;
+                        var iO = i*order;
+                        var ip1O = ip1*order;
+                        var ra = 0.0f;
+                        var sa = 0.0f;
                         for (var j = l; j <= n; j++)
                         {
-                            ra = ra + (matrixH[i, j] * matrixH[j, n - 1]);
+                            var jO = j*order;
-                            sa = sa + (matrixH[i, j] * matrixH[j, n]);
+                            var jOi = jO + i;
+                            ra = ra + (matrixH[jOi]*matrixH[nm1O + j]);
+                            sa = sa + (matrixH[jOi]*matrixH[nO + j]);
                         }
 
-                        w = matrixH[i, i] - p;
+                        w = matrixH[iO + i] - p;
 
-                        if (e[i] < 0.0f)
+                        if (e[i] < 0.0)
                         {
                             z = w;
                             r = ra;
@@ -992,49 +1015,49 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                         else
                         {
                             l = i;
-                            if (e[i] == 0.0f)
+                            if (e[i] == 0.0)
                             {
                                 var res = Cdiv(-ra, -sa, w, q);
-                                matrixH[i, n - 1] = res.Real;
+                                matrixH[nm1O + i] = res.Real;
-                                matrixH[i, n] = res.Imaginary;
+                                matrixH[nO + i] = res.Imaginary;
                             }
                             else
                             {
                                 // Solve complex equations
-                                x = matrixH[i, i + 1];
+                                x = matrixH[ip1O + i];
-                                y = matrixH[i + 1, i];
+                                y = matrixH[iO + ip1];
 
-                                float vr = ((d[i] - p) * (d[i] - p)) + (e[i] * e[i]) - (q * q);
+                                var vr = ((d[i] - p)*(d[i] - p)) + (e[i]*e[i]) - (q*q);
-                                float vi = (d[i] - p) * 2.0f * q;
+                                var vi = (d[i] - p)*2.0f*q;
                                 if ((vr == 0.0f) && (vi == 0.0f))
                                 {
-                                    vr = eps * norm * (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z));
+                                    vr = eps*norm*(float) (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z));
                                 }
 
                                 var res = Cdiv((x*r) - (z*ra) + (q*sa), (x*s) - (z*sa) - (q*ra), vr, vi);
-                                matrixH[i, n - 1] = res.Real;
+                                matrixH[nm1O + i] = res.Real;
-                                matrixH[i, n] = res.Imaginary;
+                                matrixH[nO + i] = res.Imaginary;
                                 if (Math.Abs(x) > (Math.Abs(z) + Math.Abs(q)))
                                 {
-                                    matrixH[i + 1, n - 1] = (-ra - (w * matrixH[i, n - 1]) + (q * matrixH[i, n])) / x;
+                                    matrixH[nm1O + ip1] = (-ra - (w*matrixH[nm1O + i]) + (q*matrixH[nO + i]))/x;
-                                    matrixH[i + 1, n] = (-sa - (w * matrixH[i, n]) - (q * matrixH[i, n - 1])) / x;
+                                    matrixH[nO + ip1] = (-sa - (w*matrixH[nO + i]) - (q*matrixH[nm1O + i]))/x;
                                 }
                                 else
                                 {
-                                    res = Cdiv(-r - (y * matrixH[i, n - 1]), -s - (y * matrixH[i, n]), z, q);
+                                    res = Cdiv(-r - (y*matrixH[nm1O + i]), -s - (y*matrixH[nO + i]), z, q);
-                                    matrixH[i + 1, n - 1] = res.Real;
+                                    matrixH[nm1O + ip1] = res.Real;
-                                    matrixH[i + 1, n] = res.Imaginary;
+                                    matrixH[nO + ip1] = res.Imaginary;
                                 }
                             }
 
                             // Overflow control
-                            t = Math.Max(Math.Abs(matrixH[i, n - 1]), Math.Abs(matrixH[i, n]));
+                            t = Math.Max(Math.Abs(matrixH[nm1O + i]), Math.Abs(matrixH[nO + i]));
                             if ((eps*t)*t > 1)
                             {
                                 for (var j = i; j <= n; j++)
                                 {
-                                    matrixH[j, n - 1] = matrixH[j, n - 1] / t;
+                                    matrixH[nm1O + j] /= t;
-                                    matrixH[j, n] = matrixH[j, n] / t;
+                                    matrixH[nO + j] /= t;
                                 }
                             }
                         }
@@ -1045,15 +1068,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             // Back transformation to get eigenvectors of original matrix
             for (var j = order - 1; j >= 0; j--)
             {
+                var jO = j*order;
                 for (var i = 0; i < order; i++)
                 {
                     z = 0.0f;
                     for (var k = 0; k <= j; k++)
                     {
-                        z = z + (a[(k * order) + i] * matrixH[k, j]);
+                        z = z + (a[k*order + i]*matrixH[jO + k]);
                     }
 
-                    a[(j * order) + i] = z;
+                    a[jO + i] = z;
                 }
             }
         }
@@ -1066,14 +1090,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// <param name="yreal">Real part of Y</param>
         /// <param name="yimag">Imaginary part of Y</param>
         /// <returns>Division result as a <see cref="Complex"/> number.</returns>
-        private static Complex32 Cdiv(float xreal, float ximag, float yreal, float yimag)
+        static Numerics.Complex32 Cdiv(float xreal, float ximag, float yreal, float yimag)
         {
             if (Math.Abs(yimag) < Math.Abs(yreal))
             {
-                return new Complex32((xreal + (ximag * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal))), (ximag - (xreal * (yimag / yreal))) / (yreal + (yimag * (yimag / yreal))));
+                return new Numerics.Complex32((xreal + (ximag*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))), (ximag - (xreal*(yimag/yreal)))/(yreal + (yimag*(yimag/yreal))));
             }
 
-            return new Complex32((ximag + (xreal * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag))), (-xreal + (ximag * (yreal / yimag))) / (yimag + (yreal * (yreal / yimag))));
+            return new Numerics.Complex32((ximag + (xreal*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))), (-xreal + (ximag*(yreal/yimag)))/(yimag + (yreal*(yreal/yimag))));
         }
 
         /// <summary>
@@ -1209,7 +1233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 for (var j = 0; j < order; j++)
                 {
                     value = 0;
-                    for (int i = 0; i < order; i++)
+                    for (var i = 0; i < order; i++)
                     {
                         value += ((DenseMatrix) MatrixEv).Values[(i*order) + j]*tmp[i];
                     }

src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs

@@ -80,9 +80,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
            
             IsSymmetric = true;
 
-            for (var i = 0; i < order & IsSymmetric; i++)
+            for (var i = 0; IsSymmetric && i < order; i++)
             {
-                for (var j = 0; j < order & IsSymmetric; j++)
+                for (var j = 0; IsSymmetric && j < order; j++)
                 {
                     IsSymmetric &= matrix.At(i, j) == matrix.At(j, i);
                 }

src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs

@@ -248,7 +248,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
             var matrix = _matrices["Square3x3"];
             var work = new float[18];
             var norm = Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, matrix.RowCount, matrix.ColumnCount, matrix.Values, work);
-            AssertHelpers.AlmostEqual(10.777754868246f, norm, 8);
+            AssertHelpers.AlmostEqual(10.777754868246f, norm, 6);
         }
 
         /// <summary>

src/UnitTests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs

@@ -30,7 +30,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
     using System.Numerics;
     using LinearAlgebra.Complex;
     using LinearAlgebra.Complex.Factorization;
-    using LinearAlgebra.Generic.Factorization;
     using NUnit.Framework;
 
     /// <summary>

src/UnitTests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs

@@ -30,7 +30,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
     using System.Numerics;
     using LinearAlgebra.Complex32;
     using LinearAlgebra.Complex32.Factorization;
-    using LinearAlgebra.Generic.Factorization;
     using NUnit.Framework;
     using Complex32 = Numerics.Complex32;
 
@@ -115,8 +114,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
         /// <summary>
         /// Can factorize a symmetric random square matrix.
         /// </summary> <param name="order">Matrix order.</param>
-        [Test, Ignore]
+        [Test]
-        public void CanFactorizeRandomSymmetricMatrix([Values(1, 2, 5, 10, 50, 100)] int order)
+        public void CanFactorizeRandomSymmetricMatrix([Values(1, 2, 5, 10)] int order)
         {
             var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order);
             MatrixHelpers.ForceConjugateSymmetric(matrixA);
@@ -201,13 +200,12 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
         /// Can solve a system of linear equations for a random vector and symmetric matrix (Ax=b).
         /// </summary>
         /// <param name="order">Matrix order.</param>
-        [Test, Ignore]
+        [Test]
         [TestCase(1)]
         [TestCase(2)]
         [TestCase(5)]
         [TestCase(10)]
         [TestCase(50)]
-        [TestCase(100)]
         public void CanSolveForRandomVectorAndSymmetricMatrix(int order)
         {
             var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order);
@@ -249,7 +247,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
         [TestCase(5)]
         [TestCase(10)]
         [TestCase(50)]
-        [TestCase(100)]
         public void CanSolveForRandomMatrixAndSymmetricMatrix(int order)
         {
             var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order);
@@ -292,13 +289,12 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
         /// Can solve a system of linear equations for a random vector and symmetric matrix (Ax=b) into a result vector.
         /// </summary>
         /// <param name="order">Matrix order.</param>
-        [Test, Ignore]
+        [Test]
         [TestCase(1)]
         [TestCase(2)]
         [TestCase(5)]
         [TestCase(10)]
         [TestCase(50)]
-        [TestCase(100)]
         public void CanSolveForRandomVectorAndSymmetricMatrixWhenResultVectorGiven(int order)
         {
             var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteHermitianDenseMatrix(order);

src/UnitTests/LinearAlgebraTests/Double/Factorization/EvdTests.cs

@@ -24,8 +24,6 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.LinearAlgebra.Generic;
-
 namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
 {
     using System;

src/UnitTests/LinearAlgebraTests/Single/Factorization/EvdTests.cs

@@ -77,7 +77,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization
         /// Can factorize a random square matrix.
         /// </summary>
         /// <param name="order">Matrix order.</param>
-        [Test, Ignore]
+        [Test]
         public void CanFactorizeRandomMatrix([Values(1, 2, 5, 10, 50, 100)] int order)
         {
             var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
@@ -108,8 +108,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization
         /// Can factorize a symmetric random square matrix.
         /// </summary>
         /// <param name="order">Matrix order.</param>
-        [Test, Ignore]
+        [Test]
-        public void CanFactorizeRandomSymmetricMatrix([Values(1, 2, 5, 10, 50, 100)] int order)
+        public void CanFactorizeRandomSymmetricMatrix([Values(1, 2, 5, 10, 50)] int order)
         {
             var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order);
             MatrixHelpers.ForceSymmetric(matrixA);
@@ -169,7 +169,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization
                 matrixA[i - 1, i] = 1;
                 matrixA[i + 1, i] = 1;
             }
-
             var factorEvd = matrixA.Evd();
 
             Assert.AreEqual(factorEvd.Determinant, 0);