dgemm segfaulting con matrici di ordine F di grandi dimensioni in scipy

Question

Sto tentando di calcolare A * AT in Python usando dgemm di SciPy, ma ottenendo un segfault quando A ha una dimensione di riga grande (~ 50.000) e ho passato le matrici in F -ordine. Naturalmente, la matrice risultante è molto grande, ma entrambi sgemm e passando per dgemm nelle opere C-order,

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

Qualcuno ha sperimentato questo bug prima o hanno alcuna idea che cosa è la causa? Sto usando Python 2.7.3, numpy 1.8.0 e scipy 0.13.2.

MODIFICA: FWIW, questo è l'unico ordine che produce l'errore.

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

Entrambe le operazioni hanno esito positivo.

EDIT: informazioni 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'] 

Answer 1

Non sei autorizzato ad argomenti alias quando chiamata Fortran. Non sono sicuro che questo sia il tuo problema, ma potrebbe esserlo.

Le prime due chiamate BLAS non alias gli argomenti perché gli array temporanei vengono creati prima di chiamare fortran. Cioè, a causa della mancata corrispondenza del dtype e dell'ordine del C, rispettivamente.

La terza chiamata BLAS alias gli argomenti. Prova con b = A.copy(). T invece.