8

I'm attempting to compute A*A.T in Python using SciPy's dgemm, but getting a segfault when A has large row dimension (~50,000) and I pass the matrices in in F-order. Of course, the resulting matrix is very large, but both sgemm and passing to dgemm in C-order works,

>>> 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)

Has anyone experienced this bug before or have any idea whats causing it? I'm using Python 2.7.3, numpy 1.8.0 and scipy 0.13.2.

EDIT: FWIW, that is the only ordering that produces the error.

>>> 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)

Both of the above succeed.

EDIT: BLAS info

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']
Brielin Brown
  • 81
  • 3

1 Answers1

1

You are not allowed to alias arguments when calling Fortran. I am not sure if this is your problem, but it might be.

The first two BLAS calls do not alias the arguments because temporary arrays are made before calling fortran. That is, due to dtype mismatch and C-ordering, respectively.

The third BLAS call aliases the arguments. Try with b=A.copy().T instead.

Sturla Molden
  • 1,124
  • 7
  • 7