2337.sagews

%cython
cimport cython
import numpy as np 
cimport numpy as np 
ctypedef double dtype_t

cdef extern from "cblas.h": 
    double ddot "cblas_ddot"(int N, double *X, int incX, double *Y, int incY)
                                                   
@cython.boundscheck(False) 
@cython.wraparound(False)
def matmul(np.ndarray[dtype_t, ndim=2] A,
           np.ndarray[dtype_t, ndim=2] B,
           np.ndarray[dtype_t, ndim=2] C=None):
    for i in range(A.shape[0]):
        for j in range(B.shape[1]): 
            s=0
            for k in range(A.shape[1]): 
                s += A[i, k] * B[k, j]
        C[i,j] = s 
    return C

@cython.boundscheck(False) 
@cython.wraparound(False)    
def matmul_gedp(np.ndarray[dtype_t, ndim=2] A,
            np.ndarray[dtype_t, ndim=2] B,
            np.ndarray[dtype_t, ndim=2] out): 
    cdef Py_ssize_t i, j, k
    cdef dtype_t s 
    for i in range(A.shape[0]):
        for j in range(B.shape[1]):
            out[i,j] = ddot(A.shape[1], 
                            <dtype_t*>(A.data + i*A.strides[0]),
                            A.strides[1] // sizeof(dtype_t), 
                            <dtype_t*>(B.data + j*B.strides[1]), 
                            B.strides[0] // sizeof(dtype_t))
    return out

from numpy import array, random
M = N = K = 100
ALPHA = 1.0
BETA = 0.0
LDA = LDB = LDC = 1
A = array([[random.random() for k in xrange(K)] for m in xrange(M)])
B = array([[random.random() for n in xrange(N)] for k in xrange(K)])
C = array([[random.random() for n in xrange(N)] for m in xrange(M)])

timeit('matmul(A, B, C)', number=1)
timeit('matmulgedp(A, B, C)', number=1)

1 loops, best of 3: 966 ms per loop 1 loops, best of 3: 1.66 ms per loop