dgemm segfaulting 与 scipy 中的大 F 阶矩阵答案

【问题标题】：dgemm segfaulting with large F-order matrices in scipydgemm segfaulting 与 scipy 中的大 F 阶矩阵
【发布时间】：2013-12-28 09:27:44
【问题描述】：

我正在尝试使用 SciPy 的 dgemm 在 Python 中计算 A*A.T，但是当 A 具有较大的行维度（~50,000）并且我以 F 顺序传递矩阵时会出现段错误。当然，生成的矩阵非常大，但是 sgemm 和在 C 阶中传递给 dgemm 都有效，

>>> import numpy as np
>>> import scipy.linalg.blas
>>> A = np.ones((50000,100))
#sgemm works, A.T is in F-order
>>> C = scipy.linalg.blas.sgemm(alpha=1.0, a=A.T, b=A.T, trans_a=True);
#dgemm works, A is in C-order (slower)
>>> C = scipy.linalg.blas.dgemm(alpha=1.0, a=A, b=A, trans_b=True); 
#dgemm segfaults when both are in F order
>>> C = scipy.linalg.blas.dgemm(alpha=1.0, a=A.T, b=A.T, trans_a=True);
Segmentation fault (core dumped)

之前有没有人遇到过这个错误，或者知道是什么原因造成的？我正在使用 Python 2.7.3、numpy 1.8.0 和 scipy 0.13.2。

编辑：FWIW，这是唯一产生错误的顺序。

>>> C = scipy.linalg.blas.dgemm(alpha=1.0, a=A.T, b=A, trans_a=True, trans_b=True)
>>> C = scipy.linalg.blas.dgemm(alpha=1.0, a=A, b=A.T)

以上都成功了。

编辑：BLAS 信息

blas_opt_info:
libraries = ['ptf77blas', 'ptcblas', 'atlas']
library_dirs = ['/usr/lib/atlas-base']
define_macros = [('ATLAS_INFO', '"\\"3.8.4\\""')]
language = c
include_dirs = ['/usr/include/atlas']

【问题讨论】：

您使用的是哪个 blas 库？ np.__config__.blas_info.
编辑描述以包含信息

标签： python numpy matrix scipy blas

【解决方案1】：

在调用 Fortran 时不允许给参数起别名。我不确定这是否是您的问题，但可能是。

前两个 BLAS 调用没有为参数设置别名，因为临时数组是在调用 fortran 之前创建的。也就是说，分别是由于 dtype 不匹配和 C 排序。

第三个 BLAS 调用为参数设置别名。请尝试使用 b=A.copy().T。

【讨论】：