从 Eigen 调用 Lapack 函数时出现分段错误答案

【问题标题】：Segmentation Fault while Calling Lapack function from Eigen从 Eigen 调用 Lapack 函数时出现分段错误
【发布时间】：2014-08-11 18:19:14
【问题描述】：

我的程序是用 C++ 编写的，我使用 Eigen 库进行其中的矩阵运算。里面有一个巨大的矩阵积，尺寸分别为50000*1000和1000*50000。所以我尝试从 MKL 库中调用 BLAS 函数来提高性能。然后调用 dgemm 函数时出现分段错误。

这里是代码

        double alpha = 1, beta = 0;
    double *s1;
    MKL_INT mm1 = q, nn1 = q, kk1 = ncol1;
    s1 = (double *)malloc(q*q*sizeof(double));
    cout << 14 << endl;

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,mm1, nn1, kk1, alpha, V.data(), mm1, A01.data(), kk1, beta, s1, mm1);

代码在小尺寸下工作正常。我用以下代码编译代码：

icpc lapack.cpp generators.cpp SimpleRNG.cpp example.cpp -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm  -DMKL_ILP64 -o new_example.o

或

icpc lapack.cpp generators.cpp SimpleRNG.cpp example.cpp -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -o new_example.o

即：我尝试了 LP64 接口和 ILP64 接口，但它们都不起作用，有人可以帮我吗？我在Linux服务器上运行程序，内存很大。

非常感谢！

【问题讨论】：

50000 * 50000 * sizeof(double) 约为 20GB（假设为 8 字节 doubles）。你检查你的malloc的返回值了吗？
这个问题标记为C++，你为什么用malloc而不是std::vector？此外，鉴于您使用的是malloc，那么对free 的调用在哪里？另外：你为什么使用 C 风格的演员表？
@T.C.，谢谢提醒，我把内存分配部分改成'double s1 = (double)mkl_malloc(sizeof(double)*nn1*nn1, 64)'，并使用 ilp 接口编译，现在代码可以工作了。你介意回答这个问题吗，我会马上接受的。
@RobertAllanHenniganLeahy，感谢您的评论，这里我使用 CBLAS 作为接口来调用用 fortran 编写的 BLAS 函数。 CBLAS是C接口，所以我使用malloc()，函数末尾有free()。
我使用bpeek查看服务器的实时输出，有时输出很长，我需要使用什么命令将输出重定向到服务器上的文件，所以我可以检查?谢谢！

标签： c++ segmentation-fault eigen lapack blas

【解决方案1】：

以下讨论假设：

sizeof(double) == 8
MKL_INT 是 int 和 sizeof(int) == 4
sizeof(std::size_t) == 8
CHAR_BIT == 8

这些在典型的 64 位系统中应该是正确的。

这一行发生了一些非常有趣的事情：

s1 = (double *)malloc(q*q*sizeof(double));

如果q 是50000，那么q*q 是2500000000。如果q 是int，那么这会导致有符号整数溢出，从而导致未定义的行为。在这种特殊情况下，编译器可能会简单地回绕（有效地减去 2³²），从而产生-1794967296。

但是，当您将-1794967296 与sizeof(double) 相乘时，std::size_t 是一个无符号整数类型，有趣的事情就会发生。如果 size_t 是 64 位，则编译器需要将 -1794967296 转换为无符号的 64 位数字，从概念上讲，这种转换是通过将 2⁶⁴ 加到数字上来实现的，从而得到 @987654341 @。当你将它乘以sizeof(double) 时，它再次溢出，但无符号溢出是明确定义的，并且对于 64 位操作数返回结果模 2⁶⁴，因此最终结果是 18446744059349813248。（有关计算，请参见此处的demo）。

因此，您的原始代码最终会尝试分配 18446744059349813248 字节的内存（这几乎是 16 exabytes）。哎哟。显然分配将失败并返回一个空指针。由于您没有检查返回值，因此您稍后会遇到分段错误。

当你把它改写为

s1 = (double *)malloc(sizeof(double) * q * q);

然后首先评估sizeof(double) * q。此乘法会将q 转换为std::size_t，但由于q 是正数，因此转换不会影响其值。因此结果是明确定义的，并且是一个std::size_t，其值为400000。第二个乘法同样定义良好——q 再次转换为std::size_t，结果乘法产生20000000000，它不会溢出std::size_t，因此您的malloc 调用实际上请求了20GB 内存.

【讨论】：