Cython execution speed vs. MSVC and GCC versions
Intro
I have a fairly simple Cython module - which I simplified even further for the specific tests I have carried on.
This module has only one class, which has only one interesting method (named run): this method accepts as an input a Fortran-ordered 2D NumPy array and two 1D NumPy arrays, and does some very, very simple things on those (see below for code).
For the sake of benchmarking, I have compiled the exact same module with MSVC, GCC 8, GCC 11, GCC 12 and GCC 13 on Windows 10 64 bit, Python 3.9.10 64 bit, Cython 0.29.32. All the GCC compilers I have obtained from the excellent Brecht Sanders GitHub page (https://github.com/brechtsanders).
Main Question
The overarching question of this very long post is: I am just curious to know if anyone has any explanation regarding why GCC12 and GCC13 are so much slower than GCC11 (which is the fastest of all). Looks like performances are going down at each release of GCC, rather than getting better...
Benchmarks
In the benchmarking, I simply vary the array dimensions of the 2D and 1D arrays (m and n) and the number on nonzero entries in the 2D and 1D arrays. I repeat the run method 20 times per compiler version, per set of m and n and nonzero entries.
Optimization settings I am using:
MVSC
MSVC_EXTRA_COMPILE_ARGS = ['/O2', '/GS-', '/fp:fast', '/Ob2', '/nologo', '/arch:AVX512', '/Ot', '/GL']
GCC
GCC_EXTRA_COMPILE_ARGS = ['-Ofast', '-funroll-loops', '-flto', '-ftree-vectorize', '-march=native', '-fno-asynchronous-unwind-tables']
GCC_EXTRA_LINK_ARGS = ['-flto'] + GCC_EXTRA_COMPILE_ARGS
What I am observing is the following:
- MSVC is by far the slowest at executing the benchmark (why would that be on Windows?)
- The progression GCC8 -> GCC11 is promising, as GCC11 is faster than GCC8
- GCC12 and GCC13 are both significantly slower than GCC11, with GCC13 being the worst (twice as slow as GCC11 and much worse than GCC12)
Table of Results:
Runtimes are in milliseconds (ms)
Graph (NOTE: Logarithmic Y axis!!):
Runtimes are in milliseconds (ms)
Code
Cython file:
###############################################################################
import numpy as np
cimport numpy as np
import cython
from cython.view cimport array as cvarray
from libc.float cimport FLT_MAX
DTYPE_float = np.float32
ctypedef np.float32_t DTYPE_float_t
cdef float SMALL = 1e-10
cdef int MAXSIZE = 1000000
cdef extern from "math.h" nogil:
cdef float fabsf(float x)
###############################################################################
@cython.final
cdef class CythonLoops:
cdef int m, n
cdef int [:] col_starts
cdef int [:] row_indices
cdef double [:] x
def __cinit__(self):
self.m = 0
self.n = 0
self.col_starts = cvarray(shape=(MAXSIZE,), itemsize=sizeof(int), format='i')
self.row_indices = cvarray(shape=(MAXSIZE,), itemsize=sizeof(int), format='i')
self.x = cvarray(shape=(MAXSIZE,), itemsize=sizeof(double), format='d')
@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False) # turn off negative index wrapping for entire function
@cython.nonecheck(False)
@cython.initializedcheck(False)
@cython.cdivision(True)
cpdef run(self, DTYPE_float_t[::1, :] matrix,
DTYPE_float_t[:] ub_values,
DTYPE_float_t[:] priority):
cdef Py_ssize_t i, j, m, n
cdef int nza, collen
cdef double too_large, ok, obj
cdef float ub, element
cdef int [:] col_starts = self.col_starts
cdef int [:] row_indices = self.row_indices
cdef double [:] x = self.x
m = matrix.shape[0]
n = matrix.shape[1]
self.m = m
self.n = n
nza = 0
collen = 0
for i in range(n):
for j in range(m+1):
if j == 0:
element = priority[i]
else:
element = matrix[j-1, i]
if fabsf(element) < SMALL:
continue
if j == 0:
obj = <double>element
# Do action 1 with external library
else:
collen = nza + 1
col_starts[collen] = i+1
row_indices[collen] = j
x[collen] = <double>element
nza += 1
ub = ub_values[i]
if ub > FLT_MAX:
too_large = 0.0
# Do action 2 with external library
elif ub > SMALL:
ok = <double>ub
# Do action 3 with external library
# Use x, row_indices and col_starts in the external library
Setup file:
I use the following to compile it:
python setup.py build_ext --inplace --compiler=mingw32 gcc13
Where the last argument is the compiler I want to test
#!/usr/bin/env python
from setuptools import setup
from setuptools import Extension
from Cython.Build import cythonize
from Cython.Distutils import build_ext
import numpy as np
import os
import shutil
import sys
import getpass
MODULE = 'loop_cython_%s'
GCC_EXTRA_COMPILE_ARGS = ['-Ofast', '-funroll-loops', '-flto', '-ftree-vectorize', '-march=native', '-fno-asynchronous-unwind-tables']
GCC_EXTRA_LINK_ARGS = ['-flto'] + GCC_EXTRA_COMPILE_ARGS
MSVC_EXTRA_COMPILE_ARGS = ['/O2', '/GS-', '/fp:fast', '/Ob2', '/nologo', '/arch:AVX512', '/Ot', '/GL']
MSVC_EXTRA_LINK_ARGS = MSVC_EXTRA_COMPILE_ARGS
def remove_builds(kind):
for folder in ['build', 'bin']:
if os.path.isdir(folder):
if folder == 'bin':
continue
shutil.rmtree(folder, ignore_errors=True)
if os.path.isfile(MODULE + '_%s.c'%kind):
os.remove(MODULE + '_%s.c'%kind)
def setenv(extra_args, doset=True, path=None, kind='gcc8'):
flags = ''
if doset:
flags = ' '.join(extra_args)
for key in ['CFLAGS', 'FFLAGS', 'CPPFLAGS']:
os.environ[key] = flags
user = getpass.getuser()
if doset:
path = os.environ['PATH']
if kind == 'gcc8':
os.environ['PATH'] = r'C:\Users\%s\Tools\MinGW64_8.0\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\Users\%s\WinPython39\WPy64-39100\python-3.9.10.amd64;'%(user, user)
elif kind == 'gcc11':
os.environ['PATH'] = r'C:\Users\%s\Tools\MinGW64\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\Users\%s\WinPython39\WPy64-39100\python-3.9.10.amd64;'%(user, user)
elif kind == 'gcc12':
os.environ['PATH'] = r'C:\Users\%s\Tools\MinGW64_12.2.0\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\Users\%s\WinPython39\WPy64-39100\python-3.9.10.amd64;'%(user, user)
elif kind == 'gcc13':
os.environ['PATH'] = r'C:\Users\%s\Tools\MinGW64_13.2.0\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\Users\%s\WinPython39\WPy64-39100\python-3.9.10.amd64;'%(user, user)
elif kind == 'msvc':
os.environ['PATH'] = r'C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.35.32215\bin\Hostx64\x64;C:\WINDOWS\system32;C:\WINDOWS;C:\Users\J0514162\WinPython39\WPy64-39100\python-3.9.10.amd64;C:\Program Files (x86)\Windows Kits\10\bin\10.0.22000.0\x64'
os.environ['LIB'] = r'C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.35.32215\lib\x64;C:\Program Files (x86)\Windows Kits\10\Lib\10.0.22000.0\um\x64;C:\Program Files (x86)\Windows Kits\10\Lib\10.0.22000.0\ucrt\x64'
os.environ["DISTUTILS_USE_SDK"] = '1'
os.environ["MSSdk"] = '1'
else:
os.environ['PATH'] = path
return path
class CustomBuildExt(build_ext):
def build_extensions(self):
# Override the compiler executables. Importantly, this
# removes the "default" compiler flags that would
# otherwise get passed on to to the compiler, i.e.,
# distutils.sysconfig.get_var("CFLAGS").
self.compiler.set_executable("compiler_so", "gcc -mdll -O -Wall -DMS_WIN64")
self.compiler.set_executable("compiler_cxx", "g++ -O -Wall -DMS_WIN64")
self.compiler.set_executable("linker_so", "gcc -shared -static")
self.compiler.dll_libraries = []
build_ext.build_extensions(self)
if __name__ == '__main__':
os.system('cls')
kind = None
for arg in sys.argv:
if arg.strip() in ['gcc8', 'gcc11', 'gcc12', 'gcc13', 'msvc']:
kind = arg
sys.argv.remove(arg)
break
base_file = os.path.join(os.getcwd(), MODULE[0:-3])
source = base_file + '.pyx'
target = base_file + '_%s.pyx'%kind
shutil.copyfile(source, target)
if kind == 'msvc':
extra_compile_args = MSVC_EXTRA_COMPILE_ARGS[:]
extra_link_args = MSVC_EXTRA_LINK_ARGS[:] + ['/MANIFEST']
else:
extra_compile_args = GCC_EXTRA_COMPILE_ARGS[:]
extra_link_args = GCC_EXTRA_LINK_ARGS[:]
path = setenv(extra_compile_args, kind=kind)
remove_builds(kind)
define_macros = [('WIN32', 1)]
nname = MODULE%kind
include_dirs = [np.get_include()]
if kind == 'msvc':
include_dirs += [r'C:\Program Files (x86)\Windows Kits\10\Include\10.0.22000.0\ucrt',
r'C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.35.32215\include',
r'C:\Program Files (x86)\Windows Kits\10\Include\10.0.22000.0\shared']
extensions = [
Extension(nname, [nname + '.pyx'],
extra_compile_args=extra_compile_args,
extra_link_args=extra_link_args,
include_dirs=include_dirs,
define_macros=define_macros)]
# build the core extension(s)
setup_kwargs = {'ext_modules': cythonize(extensions,
compiler_directives={'embedsignature' : False,
'boundscheck' : False,
'wraparound' : False,
'initializedcheck': False,
'cdivision' : True,
'language_level' : '3str',
'nonecheck' : False},
force=True,
cache=False,
quiet=False)}
if kind != 'msvc':
setup_kwargs['cmdclass'] = {'build_ext': CustomBuildExt}
setup(**setup_kwargs)
setenv([], False, path)
remove_builds(kind)
Test code:
import os
import numpy
import time
import loop_cython_msvc as msvc
import loop_cython_gcc8 as gcc8
import loop_cython_gcc11 as gcc11
import loop_cython_gcc12 as gcc12
import loop_cython_gcc13 as gcc13
# M N NNZ(matrix) NNZ(priority) NNZ(ub)
DIMENSIONS = [(1661 , 2608 , 3560 , 375 , 2488 ),
(2828 , 3512 , 4333 , 413 , 2973 ),
(780 , 985 , 646 , 23 , 984 ),
(799 , 1558 , 1883 , 301 , 1116 ),
(399 , 540 , 388 , 44 , 517 ),
(10545, 10486, 14799 , 1053 , 10041),
(3369 , 3684 , 3684 , 256 , 3242 ),
(2052 , 5513 , 4772 , 1269 , 3319 ),
(224 , 628 , 1345 , 396 , 594 ),
(553 , 1475 , 1315 , 231 , 705 )]
def RunTest():
print('M N NNZ MSVC GCC 8 GCC 11 GCC 12 GCC 13')
for m, n, nnz_mat, nnz_priority, nnz_ub in DIMENSIONS:
print('%-6d %-6d %-8d'%(m, n, nnz_mat), end='')
for solver, label in zip([msvc, gcc8, gcc11, gcc12, gcc13], ['MSVC', 'GCC 8', 'GCC 11', 'GCC 12', 'GCC 13']):
numpy.random.seed(123456)
size = m*n
idxes = numpy.arange(size)
matrix = numpy.zeros((size, ), dtype=numpy.float32)
idx_mat = numpy.random.choice(idxes, nnz_mat)
matrix[idx_mat] = numpy.random.uniform(0, 1000, size=(nnz_mat, ))
matrix = numpy.asfortranarray(matrix.reshape((m, n)))
idxes = numpy.arange(m)
priority = numpy.zeros((m, ), dtype=numpy.float32)
idx_pri = numpy.random.choice(idxes, nnz_priority)
priority[idx_pri] = numpy.random.uniform(0, 1000, size=(nnz_priority, ))
idxes = numpy.arange(n)
ub_values = numpy.inf*numpy.ones((n, ), dtype=numpy.float32)
idx_ub = numpy.random.choice(idxes, nnz_ub)
ub_values[idx_ub] = numpy.random.uniform(0, 1000, size=(nnz_ub, ))
solver = solver.CythonLoops()
time_tot = []
for i in range(20):
start = time.perf_counter()
solver.run(matrix, ub_values, priority)
elapsed = time.perf_counter() - start
time_tot.append(elapsed*1e3)
print('%-8.4g'%numpy.mean(time_tot), end=' ')
print()
if __name__ == '__main__':
os.system('cls')
RunTest()
EDIT
After @PeterCordes comments, I have changed the optimization flags to this:
MSVC_EXTRA_COMPILE_ARGS = ['/O2', '/GS-', '/fp:fast', '/Ob2', '/nologo', '/arch:AVX512', '/Ot', '/GL', '/QIntel-jcc-erratum']
GCC_EXTRA_COMPILE_ARGS = ['-Ofast', '-funroll-loops', '-flto', '-ftree-vectorize', '-march=native', '-fno-asynchronous-unwind-tables', '-Wa,-mbranches-within-32B-boundaries']
MSVC appears to be marginally faster than before (between 5% and 10%), but GCC12 and GCC13 are slower than before (between 3% and 20%). Below a graph with the results on the largest 2D matrix:
Note: "Current" means with the latest optimization flags suggested by @PeterCordes, "Previous" is the original set of flags.
This is a partial answer providing the generated C code produced by Cython once it has been simplified a bit to be shorter, more human-readable and easy to compile without any Cython, Python, or NumPy dependencies (the transformations are not expected to drastically impact the timings). It also show the generated assembly code and few comments.
Here is the C code:
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <float.h>
/* Remove external dependencies */
typedef int64_t Py_ssize_t;
typedef void memoryview_obj;
typedef void* VTable;
typedef struct {
struct memoryview_obj *memview;
char *data;
Py_ssize_t shape[8];
Py_ssize_t strides[8];
Py_ssize_t suboffsets[8];
} memviewslice;
struct Self
{
/* PyObject_HEAD */
struct VTable* tab;
int m;
int n;
memviewslice col_starts;
memviewslice row_indices;
memviewslice x;
};
static float test_SMALL = 1e-10;
extern void INC_MEMVIEW(memviewslice* slice, int count);
extern void XDEC_MEMVIEW(memviewslice* slice, int count);
void run(struct Self* self, memviewslice matrix, memviewslice ub_values, memviewslice priority)
{
Py_ssize_t i;
Py_ssize_t j;
Py_ssize_t m;
Py_ssize_t n;
int nza;
int collen;
double too_large;
double ok;
double obj;
float ub;
float element;
memviewslice col_starts = { 0, 0, { 0 }, { 0 }, { 0 } };
memviewslice row_indices = { 0, 0, { 0 }, { 0 }, { 0 } };
memviewslice x = { 0, 0, { 0 }, { 0 }, { 0 } };
memviewslice t_1 = { 0, 0, { 0 }, { 0 }, { 0 } };
memviewslice t_2 = { 0, 0, { 0 }, { 0 }, { 0 } };
Py_ssize_t t_3;
Py_ssize_t t_4;
Py_ssize_t t_6;
Py_ssize_t t_7;
int t_9;
Py_ssize_t t_10;
Py_ssize_t t_11;
t_1 = self->col_starts;
INC_MEMVIEW(&t_1, 1);
col_starts = t_1;
t_1.memview = NULL;
t_1.data = NULL;
t_1 = self->row_indices;
INC_MEMVIEW(&t_1, 1);
row_indices = t_1;
t_1.memview = NULL;
t_1.data = NULL;
t_2 = self->x;
INC_MEMVIEW(&t_2, 1);
x = t_2;
t_2.memview = NULL;
t_2.data = NULL;
m = matrix.shape[0];
n = matrix.shape[1];
self->m = m;
self->n = n;
nza = 0;
collen = 0;
t_3 = n;
t_4 = t_3;
for (i = 0; i < t_4; i++)
{
t_6 = m + 1;
t_7 = t_6;
for (j = 0; j < t_7; j++)
{
t_9 = (j == 0) != 0;
if (t_9)
{
t_10 = i;
element = (*((float*) ( /* dim=0 */ (priority.data + t_10 * priority.strides[0]) )));
goto L7;
}
{
t_10 = j - 1;
t_11 = i;
element = (*((float*) ( /* dim=1 */ (( /* dim=0 */ ((char *) (((float*) matrix.data) + t_10)) ) + t_11 * matrix.strides[1]) )));
}
L7:;
t_9 = (fabsf(element) < test_SMALL) != 0;
if (t_9)
{
goto L5_continue;
}
t_9 = ((j == 0) != 0);
if (t_9)
{
obj = (double)element;
goto L9;
}
{
collen = nza + 1;
t_11 = collen;
*((int *) ( /* dim=0 */ (col_starts.data + t_11 * col_starts.strides[0]) )) = i + 1;
t_11 = collen;
*((int *) ( /* dim=0 */ (row_indices.data + t_11 * row_indices.strides[0]) )) = j;
t_11 = collen;
*((double *) ( /* dim=0 */ (x.data + t_11 * x.strides[0]) )) = ((double)element);
nza = nza + 1;
}
L9:;
L5_continue:;
}
t_11 = i;
ub = (*((float*) ( /* dim=0 */ (ub_values.data + t_11 * ub_values.strides[0]) )));
t_9 = (ub > FLT_MAX) != 0;
if (t_9)
{
too_large = 0.0;
goto L10;
}
t_9 = (ub > test_SMALL) != 0;
if (t_9)
{
ok = (double)ub;
}
L10:;
}
XDEC_MEMVIEW(&col_starts, 1);
XDEC_MEMVIEW(&row_indices, 1);
XDEC_MEMVIEW(&x, 1);
}
Here is the resulting main loop compiled in assembly using GCC 11.4 (with the compiler option -Ofast -funroll-loops -ftree-vectorize -march=skylake -fno-asynchronous-unwind-tables):
.L4:
lea rdi, [r12-1]
xor eax, eax
and edi, 3
je .L5
vmovss xmm4, DWORD PTR [rcx]
mov eax, 1
vandps xmm5, xmm4, xmm1
vcomiss xmm9, xmm5
ja .L25
inc edx
movsx r14, edx
mov r15, r9
imul r15, r14
vcvtss2sd xmm6, xmm4, xmm4
mov DWORD PTR [r8+r15], esi
mov r15, r11
imul r15, r14
imul r14, rbp
mov DWORD PTR [r10+r15], 1
vmovsd QWORD PTR [rbx+r14], xmm6
.L25:
cmp rdi, 1
je .L5
cmp rdi, 2
je .L26
inc rax
vmovss xmm7, DWORD PTR [rcx-4+rax*4]
vandps xmm2, xmm7, xmm1
vcomiss xmm9, xmm2
ja .L26
inc edx
movsx rdi, edx
mov r14, r9
mov r15, r11
imul r14, rdi
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r8+r14], esi
vcvtss2sd xmm3, xmm7, xmm7
mov DWORD PTR [r10+r15], eax
vmovsd QWORD PTR [rbx+rdi], xmm3
.L26:
inc rax
vmovss xmm6, DWORD PTR [rcx-4+rax*4]
vandps xmm8, xmm6, xmm1
vcomiss xmm9, xmm8
jbe .L36
jmp .L5
.L7:
vmovss xmm10, DWORD PTR [rcx-4+rax*4]
vandps xmm11, xmm10, xmm1
vcomiss xmm9, xmm11
ja .L24
inc edx
movsx rdi, edx
mov r14, r9
mov r15, r11
imul r14, rdi
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r8+r14], esi
vcvtss2sd xmm12, xmm10, xmm10
mov DWORD PTR [r10+r15], eax
vmovsd QWORD PTR [rbx+rdi], xmm12
.L24:
lea r14, [rax+1]
vmovss xmm13, DWORD PTR [rcx-4+r14*4]
vandps xmm14, xmm13, xmm1
vcomiss xmm9, xmm14
ja .L27
inc edx
movsx rdi, edx
mov r15, r9
imul r15, rdi
vcvtss2sd xmm15, xmm13, xmm13
mov DWORD PTR [r8+r15], esi
mov r15, r11
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r10+r15], r14d
vmovsd QWORD PTR [rbx+rdi], xmm15
.L27:
lea r14, [rax+2]
vmovss xmm0, DWORD PTR [rcx-4+r14*4]
vandps xmm4, xmm0, xmm1
vcomiss xmm9, xmm4
ja .L28
inc edx
movsx rdi, edx
mov r15, r9
imul r15, rdi
vcvtss2sd xmm5, xmm0, xmm0
mov DWORD PTR [r8+r15], esi
mov r15, r11
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r10+r15], r14d
vmovsd QWORD PTR [rbx+rdi], xmm5
.L28:
add rax, 3
vmovss xmm6, DWORD PTR [rcx-4+rax*4]
vandps xmm7, xmm6, xmm1
vcomiss xmm9, xmm7
ja .L5
.L36:
inc edx
movsx rdi, edx
mov r14, r9
mov r15, r11
imul r14, rdi
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r8+r14], esi
vcvtss2sd xmm2, xmm6, xmm6
mov DWORD PTR [r10+r15], eax
vmovsd QWORD PTR [rbx+rdi], xmm2
.L5:
inc rax
cmp r12, rax
jne .L7
lea rax, [rsi+1]
add rcx, QWORD PTR [rsp+8]
cmp r13, rsi
je .L2
mov rsi, rax
jmp .L4
.L2:
Same with GCC 12.2:
.L4:
lea rdi, [r12-1]
xor eax, eax
and edi, 3
je .L5
vmovss xmm15, DWORD PTR [rcx]
mov eax, 1
vandps xmm4, xmm15, xmm1
vcomiss xmm0, xmm4
ja .L27
inc edx
movsx r14, edx
mov r15, r9
imul r15, r14
vcvtss2sd xmm5, xmm15, xmm15
mov DWORD PTR [r8+r15], esi
mov r15, r11
imul r15, r14
imul r14, rbp
mov DWORD PTR [r10+r15], 1
vmovsd QWORD PTR [rbx+r14], xmm5
.L27:
cmp rdi, 1
je .L5
cmp rdi, 2
je .L28
inc rax
vmovss xmm6, DWORD PTR [rcx-4+rax*4]
vandps xmm7, xmm6, xmm1
vcomiss xmm0, xmm7
ja .L28
inc edx
movsx rdi, edx
mov r14, r9
mov r15, r11
imul r14, rdi
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r8+r14], esi
vcvtss2sd xmm2, xmm6, xmm6
mov DWORD PTR [r10+r15], eax
vmovsd QWORD PTR [rbx+rdi], xmm2
.L28:
inc rax
vmovss xmm5, DWORD PTR [rcx-4+rax*4]
vandps xmm3, xmm5, xmm1
vcomiss xmm0, xmm3
jbe .L39
jmp .L5
.L41:
vmovss xmm8, DWORD PTR [rcx-4+rax*4]
vandps xmm9, xmm8, xmm1
vcomiss xmm0, xmm9
ja .L26
inc edx
movsx r15, edx
mov r14, r9
mov rdi, r11
imul r14, r15
imul rdi, r15
imul r15, rbp
mov DWORD PTR [r8+r14], esi
vcvtss2sd xmm10, xmm8, xmm8
mov DWORD PTR [r10+rdi], eax
vmovsd QWORD PTR [rbx+r15], xmm10
.L26:
lea r14, [rax+1]
vmovss xmm11, DWORD PTR [rcx-4+r14*4]
vandps xmm12, xmm11, xmm1
vcomiss xmm0, xmm12
ja .L29
inc edx
movsx rdi, edx
mov r15, r9
imul r15, rdi
vcvtss2sd xmm13, xmm11, xmm11
mov DWORD PTR [r8+r15], esi
mov r15, r11
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r10+r15], r14d
vmovsd QWORD PTR [rbx+rdi], xmm13
.L29:
lea r14, [rax+2]
vmovss xmm14, DWORD PTR [rcx-4+r14*4]
vandps xmm15, xmm14, xmm1
vcomiss xmm0, xmm15
ja .L30
inc edx
movsx rdi, edx
mov r15, r9
imul r15, rdi
vcvtss2sd xmm4, xmm14, xmm14
mov DWORD PTR [r8+r15], esi
mov r15, r11
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r10+r15], r14d
vmovsd QWORD PTR [rbx+rdi], xmm4
.L30:
add rax, 3
vmovss xmm5, DWORD PTR [rcx-4+rax*4]
vandps xmm6, xmm5, xmm1
vcomiss xmm0, xmm6
ja .L5
.L39:
inc edx
movsx rdi, edx
mov r14, r9
mov r15, r11
imul r14, rdi
imul r15, rdi
imul rdi, rbp
mov DWORD PTR [r8+r14], esi
vcvtss2sd xmm7, xmm5, xmm5
mov DWORD PTR [r10+r15], eax
vmovsd QWORD PTR [rbx+rdi], xmm7
.L5:
inc rax
cmp r12, rax
jne .L41
mov rdi, QWORD PTR [rsp+8]
lea rax, [rsi+1]
add rcx, rdi
cmp rsi, r13
je .L2
mov rsi, rax
jmp .L4
.L2:
Here is the same with GCC 13.1 (GCC 13.2 is not available on Godbolt):
.L4:
mov QWORD PTR [rsp+16], r14
xor eax, eax
.L20:
mov r8, rax
sub r8, r9
not r8
and r8d, 3
je .L6
inc rax
vmovss xmm15, DWORD PTR [rcx-4+rax*4]
vandps xmm0, xmm15, xmm2
vcomiss xmm1, xmm0
ja .L20
inc edx
mov r15, r10
mov rdi, QWORD PTR [rsp+24]
vcvtss2sd xmm4, xmm15, xmm15
movsx r14, edx
imul r15, r14
mov DWORD PTR [rdi+r15], esi
mov r15, rbx
imul r15, r14
imul r14, r13
mov DWORD PTR [r11+r15], eax
vmovsd QWORD PTR [r12+r14], xmm4
cmp r8, 1
je .L6
cmp r8, 2
je .L21
inc rax
vmovss xmm5, DWORD PTR [rcx-4+rax*4]
vandps xmm6, xmm5, xmm2
vcomiss xmm1, xmm6
ja .L20
inc edx
mov r14, r10
vcvtss2sd xmm7, xmm5, xmm5
movsx r8, edx
imul r14, r8
mov DWORD PTR [rdi+r14], esi
mov rdi, rbx
imul rdi, r8
imul r8, r13
mov DWORD PTR [r11+rdi], eax
vmovsd QWORD PTR [r12+r8], xmm7
.L21:
inc rax
vmovss xmm5, DWORD PTR [rcx-4+rax*4]
vandps xmm3, xmm5, xmm2
vcomiss xmm1, xmm3
ja .L20
inc edx
mov r14, r10
mov rdi, QWORD PTR [rsp+24]
movsx r8, edx
imul r14, r8
.L32:
mov r15, rbx
mov DWORD PTR [rdi+r14], esi
vcvtss2sd xmm8, xmm5, xmm5
imul r15, r8
imul r8, r13
mov DWORD PTR [r11+r15], eax
vmovsd QWORD PTR [r12+r8], xmm8
.L6:
lea r8, [rax+1]
mov rax, r8
cmp r9, r8
je .L7
vmovss xmm9, DWORD PTR [rcx-4+r8*4]
vandps xmm10, xmm9, xmm2
vcomiss xmm1, xmm10
ja .L20
inc edx
mov r15, r10
mov rdi, QWORD PTR [rsp+24]
vcvtss2sd xmm11, xmm9, xmm9
movsx r14, edx
mov DWORD PTR [rsp+4], edx
imul r15, r14
mov DWORD PTR [rdi+r15], esi
mov r15, rbx
imul r15, r14
imul r14, r13
mov DWORD PTR [r11+r15], eax
inc rax
vmovss xmm12, DWORD PTR [rcx-4+rax*4]
vmovsd QWORD PTR [r12+r14], xmm11
vandps xmm13, xmm12, xmm2
vcomiss xmm1, xmm13
ja .L20
inc edx
mov r15, r10
vcvtss2sd xmm14, xmm12, xmm12
movsx r14, edx
imul r15, r14
mov DWORD PTR [rdi+r15], esi
mov r15, rbx
imul r15, r14
imul r14, r13
mov DWORD PTR [r11+r15], eax
lea rax, [r8+2]
vmovss xmm15, DWORD PTR [rcx-4+rax*4]
vmovsd QWORD PTR [r12+r14], xmm14
vandps xmm0, xmm15, xmm2
vcomiss xmm1, xmm0
ja .L20
mov edx, DWORD PTR [rsp+4]
mov r15, r10
vcvtss2sd xmm4, xmm15, xmm15
add edx, 2
movsx r14, edx
imul r15, r14
mov DWORD PTR [rdi+r15], esi
mov r15, rbx
imul r15, r14
imul r14, r13
mov DWORD PTR [r11+r15], eax
lea rax, [r8+3]
vmovss xmm5, DWORD PTR [rcx-4+rax*4]
vmovsd QWORD PTR [r12+r14], xmm4
vandps xmm6, xmm5, xmm2
vcomiss xmm1, xmm6
ja .L20
mov edx, DWORD PTR [rsp+4]
mov r14, r10
add edx, 3
movsx r8, edx
imul r14, r8
jmp .L32
.L7:
mov rdi, QWORD PTR [rsp+8]
mov r14, QWORD PTR [rsp+16]
lea rax, [rsi+1]
add rcx, rdi
cmp r14, rsi
je .L2
mov rsi, rax
jmp .L4
.L2:
Discussion
First of all, we can see that the main loop is full of conditional branches mixed with goto. Such kind of code are generally not written by humans. Compiler can expect them to be rare and so optimize them less efficiently than ones without goto (since assertions are easier to compute). The thing is the j == 0 condition is not really useful here: it can be done only once before the j-based loop. GCC does not seem to see that (especially due to the complex control-flow, amplified by Cython goto instructions). On to of that, the compiler do not know if fabsf(element) < SMALL is often true or often false, or even just hard to predict at runtime. If the loop is easy to predict (the first two cases), then a good compiler can write a more efficient code than the one produced by GCC (possibly using SIMD instructions).
Actually, I am not even sure -funroll-loops is useful here because of the complex control-flow. It makes the hot-loop significantly bigger (i.e., more space in caches) and the processor migh take more time to predict the probability of the taken/not-taken branches. This is dependent of the target processor. Having smaller code is also good for profiling and optimizing your code: simpler code are easier to optimize.
To optimize the code, the compiler make assumptions on the cost of the loops and the probability of the branches to be taken. If can assume for example that the inner loop if far more expensive than the outer loop. Heuristics are used to optimize the code and their parameters are carefully tuned so to maximize the performance of benchmark suites (and avoid regressions in code like the Linux kernel for GCC). Thus, a side effect is that your code can be slower while most code can be faster using a newer version of GCC, or most compilers. This is also one reason the performance of a given compiler is not very stable from one version to another. For example, if a matrix has a lot of zeros, then fabsf(element) < SMALL will be often true, causing the outer-loop to be more expensive than expected. The thing is the outer loop might less well optimized in GCC 13.2 than GCC 12.2, resulting in a slower execution time with GCC 13.2. For example, GCC 13.2 starts the outer loop with:
mov r8, rax
sub r8, r9
not r8
and r8d, 3
je .L6
Note that modern processors, like Skylake, can execute multiple instructions in parallel per cycle. This sequence of instructions can only be executed sequentially because of the dependency chain on the r8 register. Meanwhile, GCC 12.2 generate this in the beginning of the outer loop:
lea rdi, [r12-1]
xor eax, eax
and edi, 3
je .L5
The dependency chain (on rdi/edi) is apparently significantly smaller in this case. In practice, the assembly code is pretty big, so I do not expect this to be the only issue. One need to track the actual path taken for the input matrices in all versions of GCC (very tedious), or just run the code and profile it with tools like VTune (the recommended solution in this case).
We can see that there is a lot of imul instructions. Such instructions are typically due to the strides of the array that are not known at compile-time. GCC does not generate code for the strides that are 1. You should really mention this information in the memory views using ::1 if possible (as mentioned by ead in the comments) since imul instructions tends to be a bit expensive in hot loops like this (not to mention they can also increase the size of the code which is already bigger due to the partial unrolling).