Manipulating matrices in Cython

asked 2022-07-12 14:18:14 +0100

Tilpo gravatar image

I need to do low-level arithmetic with (dense, integer) matrices in Cython. Right now I'm converting back and forward between typed memoryviews and sage Matrix_integer_dense objects, and it seems that these conversions are actually somewhat of a bottleneck.

I don't really understand how to work with matrices on a low level however. For example, consider the following simple function:

from sage.rings.integer cimport Integer
from sage.rings.integer_ring import ZZ
from sage.matrix.matrix_integer_dense cimport Matrix_integer_dense
from sage.matrix.matrix_space import MatrixSpace

def func(Matrix_integer_dense M):
    cdef Py_ssize_t i,j
    cdef Integer k 
    mat_space = MatrixSpace(ZZ, M.nrows(), M.ncols(), False, None)
    cdef Matrix_integer_dense new_mat = Matrix_integer_dense(mat_space)
    for i in range(M.nrows()):
        for j in range(M.ncols()):
            k = M[i,j]*M[j,i]
            new_mat[i,j] = k
    return new_mat

Looking annotated code Cython produces when compiling, it's clear that the k = M[i,j]*M[j,i] requires a coercion which is handled by Python, and very slow.

It is really unclear to me what the underlying type of the entries of an integer matrix are. It would be really nice if I can cheaply convert a Matrix_integer_dense object to a C-array of ints of some type and work with it directly. It seems that under the hood an integer matrix is a fancy container for a FLINT matrix, but from reading the source code it's not really clear to me how to extract this.

edit retag flag offensive close merge delete