我希望标题中问题的答案是我在做一些愚蠢的事情！
这是问题所在。我想计算一个实数对称矩阵的所有特征值和特征向量。我已经使用GNU Scientific Library在 MATLAB（实际上，我使用 Octave 运行它）和 C++ 中实现了代码。我在下面为这两种实现提供了我的完整代码。
据我所知，GSL 附带了它自己的 BLAS API 实现（以下我将其称为 GSLCBLAS）并使用我编译的这个库：
g++ -O3 -lgsl -lgslcblas
GSL在此建议使用替代 BLAS 库，例如自优化ATLAS库，以提高性能。我正在运行 Ubuntu 12.04，并且已经从Ubuntu 存储库安装了 ATLAS 软件包。在这种情况下，我编译使用：
g++ -O3 -lgsl -lcblas -latlas -lm
对于这三种情况，我以 100 的步长对大小为 100 到 1000 的随机生成的矩阵进行了实验。对于每个大小，我使用不同的矩阵执行 10 次特征分解，并平均花费的时间。结果如下：
性能上的差异是荒谬的。对于大小为 1000 的矩阵，Octave 在不到一秒的时间内执行分解；GSLCBLAS 和 ATLAS 大约需要 25 秒。
我怀疑我可能错误地使用了 ATLAS 库。欢迎任何解释；提前致谢。
关于代码的一些注释：
在 C++ 实现中，不需要使矩阵对称，因为该函数只使用它的下三角部分。
在 Octave 中，线
triu(A) + triu(A, 1)'强制矩阵是对称的。如果你想在你自己的 Linux 机器上编译 C++ 代码，你还需要添加 flag
-lrt，因为这个clock_gettime函数。不幸的是，我认为不会
clock_gettime在其他平台上退出。考虑将其更改为gettimeofday.
八度码
K = 10;
fileID = fopen('octave_out.txt','w');
for N = 100:100:1000
AverageTime = 0.0;
for k = 1:K
A = randn(N, N);
A = triu(A) + triu(A, 1)';
tic;
eig(A);
AverageTime = AverageTime + toc/K;
end
disp([num2str(N), " ", num2str(AverageTime), "\n"]);
fprintf(fileID, '%d %f\n', N, AverageTime);
end
fclose(fileID);
C++ 代码
#include <iostream>
#include <fstream>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
int main()
{
const int K = 10;
gsl_rng * RandomNumberGenerator = gsl_rng_alloc(gsl_rng_default);
gsl_rng_set(RandomNumberGenerator, 0);
std::ofstream OutputFile("atlas.txt", std::ios::trunc);
for (int N = 100; N <= 1000; N += 100)
{
gsl_matrix* A = gsl_matrix_alloc(N, N);
gsl_eigen_symmv_workspace* EigendecompositionWorkspace = gsl_eigen_symmv_alloc(N);
gsl_vector* Eigenvalues = gsl_vector_alloc(N);
gsl_matrix* Eigenvectors = gsl_matrix_alloc(N, N);
double AverageTime = 0.0;
for (int k = 0; k < K; k++)
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
gsl_matrix_set(A, i, j, gsl_ran_gaussian(RandomNumberGenerator, 1.0));
}
}
timespec start, end;
clock_gettime(CLOCK_MONOTONIC_RAW, &start);
gsl_eigen_symmv(A, Eigenvalues, Eigenvectors, EigendecompositionWorkspace);
clock_gettime(CLOCK_MONOTONIC_RAW, &end);
double TimeElapsed = (double) ((1e9*end.tv_sec + end.tv_nsec) - (1e9*start.tv_sec + start.tv_nsec))/1.0e9;
AverageTime += TimeElapsed/K;
std::cout << "N = " << N << ", k = " << k << ", Time = " << TimeElapsed << std::endl;
}
OutputFile << N << " " << AverageTime << std::endl;
gsl_matrix_free(A);
gsl_eigen_symmv_free(EigendecompositionWorkspace);
gsl_vector_free(Eigenvalues);
gsl_matrix_free(Eigenvectors);
}
return 0;
}
图 比较 Matlab 和 C。“+ vec”表示代码包括特征向量的计算。CPU% 是在 N=1000 时对 CPU 使用率的粗略观察，上限为 800%，尽管它们应该完全使用所有 8 个内核。Matlab 和 C 之间的差距小于 8 倍。
图 比较 Mathematica 中不同的矩阵类型。程序自动选择的算法。