diff --git a/src/NativeProviders/MKL/lapack.cpp b/src/NativeProviders/MKL/lapack.cpp index d51976fe..2f466777 100644 --- a/src/NativeProviders/MKL/lapack.cpp +++ b/src/NativeProviders/MKL/lapack.cpp @@ -255,8 +255,8 @@ inline MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T return info; } -template -inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], GEES gees, TREVC trevc) +template +inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc) { T* clone_a = Clone(n, n, a); T* wr = new T[n]; @@ -284,7 +284,7 @@ inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[ for (MKL_INT index = 0; index < n; ++index) { - values[index] = MKL_Complex16(wr[index], wi[index]); + values[index] = R(wr[index], wi[index]); } for (MKL_INT i = 0; i < n; ++i) @@ -343,11 +343,11 @@ inline MKL_INT eigen_complex_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 return info; } -template +template inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], SYEV syev) { T* clone_a = Clone(n, n, a); - R* w = new R[n]; + R* w = new R[n]; MKL_INT info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w); if (info != 0) @@ -668,7 +668,7 @@ extern "C" { { if (isSymmetric) { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_ssyev); + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_ssyev); } else { @@ -676,15 +676,15 @@ extern "C" { } } - DLLEXPORT MKL_INT d_eigen(bool isSymmetric, MKL_INT n, double a[], double vectors[], MKL_Complex16 values[], double d[]) + DLLEXPORT MKL_INT d_eigen(bool isSymmetric, MKL_INT n, double a[], double vectors[], MKL_Complex16 values[], double d[]) { if (isSymmetric) { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_dsyev); + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_dsyev); } else { - return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc); + return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc); } } @@ -692,20 +692,19 @@ extern "C" { { if (isSymmetric) { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_cheev); + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_cheev); } else { - return -1; - //return eigen_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); + return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc); } } - + DLLEXPORT MKL_INT z_eigen(bool isSymmetric, MKL_INT n, MKL_Complex16 a[], MKL_Complex16 vectors[], MKL_Complex16 values[], MKL_Complex16 d[]) { if (isSymmetric) { - return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_zheev); + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_zheev); } else { diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs index 7a4b8e71..4bed4a5d 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs @@ -1314,6 +1314,60 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl SafeNativeMethods.z_vector_divide(x.Length, x, y, result); } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Wether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + SafeNativeMethods.z_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs index 01322ae4..0f6e3f62 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs @@ -31,6 +31,7 @@ #if NATIVEMKL using System; +using System.Numerics; using System.Security; using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; @@ -1313,6 +1314,60 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl SafeNativeMethods.c_vector_divide(x.Length, x, y, result); } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Wether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + SafeNativeMethods.c_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + } } }