我正在尝试使用以下代码求解线性系统:
#include <stdio.h>
#include <lapacke.h>
int main () {
lapack_complex_double mat[4];
lapack_complex_double vec[2];
lapack_int p[2];
mat[0] = lapack_make_complex_double(1,0);
mat[1] = lapack_make_complex_double(1,0);
mat[2] = lapack_make_complex_double(1,0);
mat[3] = lapack_make_complex_double(-1,0);
vec[0] = lapack_make_complex_double(1,0);
vec[1] = lapack_make_complex_double(1,0);
LAPACKE_zgetrf(LAPACK_ROW_MAJOR, 2, 2, mat, 2, p);
LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N', 2, 1, mat, 2, p, vec, 2);
printf("%g %g\n", lapack_complex_double_real(vec[0]),
lapack_complex_double_imag(vec[0]));
return 0;
}
由于某些原因,这会导致 LAPACKE_zgetrs
中的非法内存访问(由 valgrind
和我的大程序在 zgetrs
中检测到,因为“glibc 检测到损坏或双重释放”)。为简洁起见,我没有将其包含在我的 SSCCE 中,但所有返回的 LAPACKE
例程都返回 0。
具有 LAPACK_COL_MAJOR
的相同代码运行且 valgrind 完美无缺。
我的lapacke、lapack等是为Ubuntu 12.04自建的。我在 lapack CMake 文件中使用了以下设置:
BUILD_COMPLEX ON
BUILD_COMPLEX16 ON
BUILD_DOUBLE ON
BUILD_SHARED_LIBS ON
BUILD_SINGLE ON
BUILD_STATIC_LIBS ON
BUILD_TESTING ON
CMAKE_BUILD_TYPE Release
LAPACKE ON
LAPACKE_WITH_TMG ON
其余的(优化的 blas/lapack 和 xblas)关闭。构建过程中没有错误,所有测试均成功。
我哪里搞砸了?
编辑: 我刚刚用 Fedora21 和打包的 lapacke 尝试了这个。它没有重现错误。
编辑 2: 虽然它不会重现内存失败,但它会产生错误的解决方案,即 (1 + 0I, 1 + 0I)
上面的输入 (应该是 (1,0)
)
最佳答案
经过更多的研究和思考,我找到了罪魁祸首:
使用 LAPACK_ROW_MAJOR
切换 ld*
前导维度参数的含义。虽然普通 Fortran 数组的主要维度是行的数量,但切换到 ROW_MAJOR
会将其含义切换为列的数量。所以正确的调用(给出正确的结果)应该是:
LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N', 2, 1, mat, 2, p, vec, 1);
其中第二个 2
是 mat
的列(不是行!)的数量,最后一个参数必须等于右侧 nrhs
(不是变量的数量!)。我将这个调用隔离开来,因为我项目中的所有其他调用都处理方矩阵,因此“错误”调用不会因对称性而产生任何负面影响。
像往常一样,如果您在末尾跳过列,前导尺寸会相应变大,就像在正常设置中跳过行一样。
显然,Fortran 文档中没有提到这一点。不幸的是,我确实在 Lapacke 文档中没有看到这样的评论,这本来可以节省我几个小时的时间。 :)
关于c - 将 LAPACKE_zgetrs 与 LAPACK_ROW_MAJOR 一起使用会导致非法内存访问,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/30212743/