tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

added native Complex EVD support

pull/290/head
Marcus Cuda 11 years ago
parent
129702d49e
commit
296ff72b4b
3 changed files with 122 additions and 14 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 27
      src/NativeProviders/MKL/lapack.cpp
  2. 54
      src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs
  3. 55
      src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs

27
src/NativeProviders/MKL/lapack.cpp
View File

@ -255,8 +255,8 @@ inline MKL_INT complex_svd_factor(bool compute_vectors, MKL_INT m, MKL_INT n, T
return info;
}
template<typename T, typename GEES, typename TREVC>
inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], GEES gees, TREVC trevc)
template<typename T, typename R, typename GEES, typename TREVC>
inline MKL_INT eigen_factor(MKL_INT n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc)
{
T* clone_a = Clone(n, n, a);
T* wr = new T[n];
@ -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<typename T, typename R, typename SYEV>
template<typename R, typename T, typename SYEV>
inline MKL_INT sym_eigen_factor(MKL_INT n, T a[], T vectors[], MKL_Complex16 values[], T d[], SYEV syev)
{
T* clone_a = Clone(n, n, a);
R* w = new R[n];
R* w = new R[n];
MKL_INT info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w);
if (info != 0)
@ -668,7 +668,7 @@ extern "C" {
{
if (isSymmetric)
{
return sym_eigen_factor<float, float>(n, a, vectors, values, d, LAPACKE_ssyev);
return sym_eigen_factor<float>(n, a, vectors, values, d, LAPACKE_ssyev);
}
else
{
@ -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<double, double>(n, a, vectors, values, d, LAPACKE_dsyev);
return sym_eigen_factor<double>(n, a, vectors, values, d, LAPACKE_dsyev);
}
else
{
return eigen_factor(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<MKL_Complex8, float>(n, a, vectors, values, d, LAPACKE_cheev);
return sym_eigen_factor<float>(n, a, vectors, values, d, LAPACKE_cheev);
}
else
{
return -1;
//return eigen_factor<MKL_Complex16, LAPACK_Z_SELECT1>(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc);
return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc);
}
}
DLLEXPORT MKL_INT z_eigen(bool isSymmetric, MKL_INT n, MKL_Complex16 a[], MKL_Complex16 vectors[], MKL_Complex16 values[], MKL_Complex16 d[])
{
if (isSymmetric)
{
return sym_eigen_factor<MKL_Complex16, double>(n, a, vectors, values, d, LAPACKE_zheev);
return sym_eigen_factor<double>(n, a, vectors, values, d, LAPACKE_zheev);
}
else
{

54
src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs
View File

@ -1314,6 +1314,60 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl
SafeNativeMethods.z_vector_divide(x.Length, x, y, result);
}
/// <summary>
/// Computes the eigenvalues and eigenvectors of a matrix.
/// </summary>
/// <param name="isSymmetric">Wether the matrix is symmetric or not.</param>
/// <param name="order">The order of the matrix.</param>
/// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
/// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
/// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
/// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
public override void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD)
{
if (matrix == null)
{
throw new ArgumentNullException("matrix");
}
if (matrix.Length != order * order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix");
}
if (matrixEv == null)
{
throw new ArgumentNullException("matrixEv");
}
if (matrixEv.Length != order * order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv");
}
if (vectorEv == null)
{
throw new ArgumentNullException("vectorEv");
}
if (vectorEv.Length != order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv");
}
if (matrixD == null)
{
throw new ArgumentNullException("matrixD");
}
if (matrixD.Length != order * order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD");
}
SafeNativeMethods.z_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD);
}
}
}

55
src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs
View File

@ -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);
}
/// <summary>
/// Computes the eigenvalues and eigenvectors of a matrix.
/// </summary>
/// <param name="isSymmetric">Wether the matrix is symmetric or not.</param>
/// <param name="order">The order of the matrix.</param>
/// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
/// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
/// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
/// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
public override void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD)
{
if (matrix == null)
{
throw new ArgumentNullException("matrix");
}
if (matrix.Length != order * order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix");
}
if (matrixEv == null)
{
throw new ArgumentNullException("matrixEv");
}
if (matrixEv.Length != order * order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv");
}
if (vectorEv == null)
{
throw new ArgumentNullException("vectorEv");
}
if (vectorEv.Length != order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv");
}
if (matrixD == null)
{
throw new ArgumentNullException("matrixD");
}
if (matrixD.Length != order * order)
{
throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD");
}
SafeNativeMethods.c_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD);
}
}
}

Write Preview
Loading…
Cancel
Save