# [SciPy-User] fast small matrix multiplication with cython?

Skipper Seabold jsseabold@gmail....
Mon Dec 6 16:34:19 CST 2010

```I'm wondering if anyone might have a look at my cython code that does
matrix multiplication and see where I can speed it up or offer some
pointers/reading.  I'm new to Cython and my knowledge of C is pretty
basic based on trial and (mostly) error, so I am sure the code is
still very naive.

import numpy as np
from matmult import dotAB, multAB

A = np.array([[ 1.,  3.,  4.],
[ 5.,  6.,  3.]])
B = A.T.copy()

timeit dotAB(A,B)
# 1 loops, best of 3: 826 ms per loop

timeit multAB(A,B)
# 1 loops, best of 3: 1.16 s per loop

As you can see my multAB results in a negative speedup of about .75.

I compile the cython code with

cython -a matmult.pyx
gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing
-I/usr/include/python2.6
-I/usr/local/lib/python2.6/dist-packages/numpy/core/include/ -o
matmult.so matmult.c

Cython code is attached and inlined below.

Profile is here (some of which I don't understand why there are
bottlenecks) http://eagle1.american.edu/~js2796a/matmult/matmult.html
-----------------------------------------------------------

from numpy cimport float64_t, ndarray, NPY_DOUBLE, npy_intp
cimport cython
from numpy import dot

ctypedef float64_t DOUBLE

cdef extern from "numpy/arrayobject.h":
cdef void import_array()
cdef object PyArray_SimpleNew(int nd, npy_intp *dims, int typenum)

import_array()

@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline object matmult(ndarray[DOUBLE, ndim=2, mode='c'] A,
ndarray[DOUBLE, ndim=2, mode='c'] B):
cdef int lda = A.shape[0]
cdef int n = B.shape[1]
cdef npy_intp *dims = [lda, n]
cdef ndarray[DOUBLE, ndim=2] out = PyArray_SimpleNew(2, dims, NPY_DOUBLE)
cdef int i,j,k
cdef double s
for i in xrange(lda):
for j in xrange(n):
s = 0
for k in xrange(A.shape[1]):
s += A[i,k] * B[k,j]
out[i,j] = s
return out

def multAB(ndarray[DOUBLE, ndim=2] A, ndarray[DOUBLE, ndim=2] B):
for i in xrange(1000000):
C = matmult(A,B)
return C

def dotAB(ndarray[DOUBLE, ndim=2] A, ndarray[DOUBLE, ndim=2] B):
for i in xrange(1000000):
C = dot(A,B)
return C

Skipper
-------------- next part --------------
A non-text attachment was scrubbed...
Name: matmult.pyx
Type: application/octet-stream
Size: 1248 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20101206/9074fb90/attachment.obj
```