我想使用来自英特尔MKL(2019 Update 2)的LAPACKE_dsyevd计算实对称矩阵的所有特征值和所有特征向量。
我在C#中使用以下方法:
public static class MKL
{
public static double[,] SymmetricEig(double[,] a, out double[] w)
{
int n1 = a.GetLength(0);
int n2 = a.GetLength(1);
if (n1 != n2) throw new System.Exception("Matrix must be square");
double[,] b = Copy(a);
int matrix_layout = 101; // row-major arrays
char jobz = 'V';
char uplo = 'U';
int n = n2;
int lda = n;
w = new double[n];
_mkl.LAPACKE_dsyevd(matrix_layout, jobz, uplo, n, b, lda, w);
return b;
}
}
同
class _mkl
{
[DllImport(DLLName, CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
internal static extern lapack_int LAPACKE_dsyevd(
int matrix_layout, char jobz, char uplo, lapack_int n, [In, Out] double[,] a, lapack_int lda, [In, Out] double[] w);
}
以及以下测试代码:
int n = 32766; // 32767 or greater --> Not enough memory to allocate work array in LAPACKE_dsyevd
double[,] A = CreateRandomSymmetricMatrix(n);
double[] w = new double[n];
double[,] B = MKL.SymmetricEig(A, out w);
同
static double[,] CreateRandomSymmetricMatrix(int n1)
{
double[,] m = new double[n1, n1];
for (int i1 = 0; i1 < n1; i1++)
{
for (int i2 = 0; i2 <= i1; i2++)
{
m[i1, i2] = r.NextDouble();
m[i2, i1] = m[i1, i2];
}
}
return m;
}
如果n
大于32766则失败并显示以下错误消息:
内存不足以在LAPACKE_dsyevd中分配工作数组
我的电脑有128 GB的RAM就足够了。我知道<gcAllowVeryLargeObjects enabled="true" />
,我把它设置为真。我很清楚C#中多维数组的65535 ^ 2限制,请参阅2d-Array with more than 65535^2 elements --> Array dimensions exceeded supported range。
顺便说一下,我可以使用MATLAB为n = 40000或更大的矩阵计算特征值分解。我认为MATLAB也在使用英特尔MKL来计算它。
那么如何使用英特尔MKL计算C#中非常大的矩阵(n> 40000)的特征值分解?
我不认为这是你的问题。以下定义表明w = new double[n];
在你的情况下是小的方式。
* LWORK (input) INTEGER
* The dimension of the array WORK.
* If N <= 1, LWORK must be at least 1.
* If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
* If JOBZ = 'V' and N > 1, LWORK must be at least
* 1 + 6*N + 2*N**2.
你真的应该总是做一个工作区查询。我知道,文档会向用户公开,但它非常方便,有助于避免出现这种情况。所以你知道,如何进行工作区查询?如果没有真正快速击中我。
这似乎是LAPACKE_dsyevd的一个错误。使用LAPACKE_dsyevr也可以使用更大的矩阵。
我在MKL
类中添加了以下行:
public static double[,] SymmetricEigRelativelyRobustRepresentations(double[,] a, out double[] w)
{
int n1 = a.GetLength(0);
int n2 = a.GetLength(1);
if (n1 != n2) throw new System.Exception("Matrix must be square");
double[,] b = Copy(a);
int matrix_layout = 101; // row-major arrays
char jobz = 'V'; // eigenvalues and eigenvectors are computed
char range = 'A'; // the routine computes all eigenvalues
char uplo = 'U'; // a stores the upper triangular part of A
int n = n2;
int lda = n;
int vl = 0;
int vu = 0;
int il = 0;
int iu = 0;
double abstol = 0;
int m = n;
w = new double[n];
double[,] z = new double[n, n];
int ldz = n;
int[] isuppz = new int[2 * n];
int info = _mkl.LAPACKE_dsyevr(matrix_layout, jobz, range, uplo, n, b, lda, vl, vu, il, iu, abstol, ref m, w, z, ldz, isuppz);
return z;
}
和以下几行到_mkl
类:
[DllImport(DLLName, CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
internal static extern lapack_int LAPACKE_dsyevr(
int matrix_layout, char jobz, char range, char uplo, lapack_int n, [In, Out] double[,] a, lapack_int lda,
double vl, double vu, lapack_int il, lapack_int iu, double abstol, [In, Out] ref lapack_int m, [In, Out] double[] w,
[In, Out] double[,] z, lapack_int ldz, [In, Out] lapack_int[] isuppz);