BLAS

Question

ile daha hızlı python iç ürünü uygulama Python'daki standart numpy doğrusal cebir rutinleri üzerinde büyük hız geliştirmeleri elde etmek için düşük düzey BLAS işlevlerini (Cython'da uygulandı) kullanarak this useful tutorial numaralı telefonu buldum. Şimdi, iyi çalışıyor vektör ürün başarıyla aldım. Önce linalg.pyx olarak aşağıdaki tasarruf:

import cython 
import numpy as np 
cimport numpy as np 

from libc.math cimport exp 
from libc.string cimport memset 

from scipy.linalg.blas import fblas 

REAL = np.float64 
ctypedef np.float64_t REAL_t 

cdef extern from "/home/jlorince/flda/voidptr.h": 
    void* PyCObject_AsVoidPtr(object obj) 

ctypedef double (*ddot_ptr) (const int *N, const double *X, const int *incX, const double *Y, const int *incY) nogil 
cdef ddot_ptr ddot=<ddot_ptr>PyCObject_AsVoidPtr(fblas.ddot._cpointer) # vector-vector multiplication 

cdef int ONE = 1 
def vec_vec(syn0, syn1, size): 
    cdef int lSize = size 
    f = <REAL_t>ddot(&lSize, <REAL_t *>(np.PyArray_DATA(syn0)), &ONE, <REAL_t *>(np.PyArray_DATA(syn1)), &ONE) 
    return f 

(mevcut voidptr.h için kaynak kodu here)

bunu derlemek sonra, gayet iyi çalışıyor ve kesinlikle daha hızlı np.inner daha:

In [1]: import linalg 
In [2]: import numpy as np 
In [3]: x = np.random.random(100) 
In [4]: %timeit np.inner(x,x) 
1000000 loops, best of 3: 1.61 µs per loop 
In [5]: %timeit linalg.vec_vec(x,x,100) 
1000000 loops, best of 3: 483 ns per loop 
In [8]: np.all(np.inner(x,x)==linalg.vec_vec(x,x,100)) 
Out[8]: True 

Şimdi, bu harika, ancak sadece iki vektörü'un nokta/iç çarpımının (bu durumda eşdeğer) hesaplanması durumunda işe yarar. Şimdi yapmam gereken şey, vektör matrisi iç ürünlerini yapmak için benzer işlevleri (benzer hızlar sunacağını umduğum) uygular. Bu 1D dizisi ve bir 2 boyutlu matris geçirilen zaman np.inner işlevselliğini çoğaltmak istediğiniz gibidir:

In [4]: x = np.random.random(5) 
In [5]: y = np.random.random((5,5)) 
In [6]: np.inner(x,y) 
Out[6]: array([ 1.42116225, 1.13242989, 1.95690196, 1.87691992, 0.93967486]) 

Bu 1D dizinin (yine 1D diziler için eşdeğer) iç/nokta ürünleri hesaplama eşdeğerdir ve matris her satır: Ben BLAS belgelerin gördüklerime

In [32]: [np.inner(x,row) for row in y] 
Out[32]: 
[1.4211622497461549, 1.1324298918119025, 1.9569019618096966,1.8769199192990056, 0.93967485730285505] 

, ben (dgemv kullanarak) böyle bir şey ile başlamak gerekir düşünüyorum:

ctypedef double (*dgemv_ptr) (const str *TRANS, const int *M, const int *N, const double *ALPHA, const double *A, const int *LDA, const double *X, const int *incX, const double *BETA, const double *Y, const int *incY) 
cdef dgemv_ptr dgemv=<dgemv>PyCObject_AsVoidPtr(fblas.dgemv._cpointer) # matrix vector multiplication 

Ancak, (a) Python'da kullanabileceğim gerçek işlevi tanımlamak için yardıma ihtiyacım var (örn. Yukarıdaki vec_vec benzeri bir vec-matrix fonksiyonu ve (b) uygulamam için modele ihtiyacım olan np.inner davranışını düzgün şekilde çoğaltmak için nasıl kullanılacağını bilmek.

Düzenleme:Here Ben kullanılarak çözülmesi gereken burada teyit edildiği o, dgemv için ilgili BLAS belgelerine link:

In [13]: np.allclose(scipy.linalg.blas.fblas.dgemv(1.0,y,x), np.inner(x,y)) 
Out[13]: True 

Ama böyle kutudan kullanmaya aslında yavaştır Saf np.inner'dan daha fazla, bu yüzden Cython uygulamasında hala yardıma ihtiyacım var.

cdef int ONE = 1 
cdef char tr = 'n' 
cdef REAL_t ZEROF = <REAL_t>0.0 
cdef REAL_t ONEF = <REAL_t>1.0 
def mat_vec(mat,vec,mat_rows,mat_cols): 
    cdef int m = mat_rows 
    cdef int n = mat_cols 
    out = <REAL_t>dgemv(&tr, &m, &n, &ONEF, <REAL_t *>(np.PyArray_DATA(mat)), &m, <REAL_t *>(np.PyArray_DATA(vec)), &ONE, &ZEROF, NULL, &ONE) 
    return out 

derleme sonra kullanıyorum (, çalıştırmayı denediğinizde:

Edit2 Burada iyi derler ama bunu çalıştırmayı denediğinizde her parçalama arızası ile piton çöküyor bu benim son girişim, var Aynı x ve y yukarıdaki gibi) ama bu sadece çöküyor. Ben yakınım, ama ... değiştirmeye başka ne

Answer 1

Başına @Pietro Saccardi'nin bilmiyorum:

int dgemv_(char *trans, integer *m, integer *n, doublereal * 
      alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, 
      doublereal *beta, doublereal *y, integer *incy) 

... 

Y  - DOUBLE PRECISION array of DIMENSION at least 
     (1 + (m - 1)*abs(INCY)) when TRANS = 'N' or 'n' 
     and at least 
     (1 + (n - 1)*abs(INCY)) otherwise. 
     Before entry with BETA non-zero, the incremented array Y 
     must contain the vector y. On exit, Y is overwritten by the 
     updated vector y. 

Sana çağrısında Y için NULL kullanabilirsiniz şüpheliyim.