使用f2py封装Fortran矩阵乘法时，增大维度至50触发matmul_r8_avx2断言失败崩溃

看起来你遇到的是GCC libgfortran库中AVX2优化版matmul函数的断言bug，这个问题在矩阵维度超过一定阈值（这里是50）时触发——因为此时libgfortran会自动切换到AVX2优化路径，但该路径对数组的类型标记检查存在问题，即便你已经用np.asfortranarray确保了Fortran内存布局，还是会触发断言失败导致崩溃。下面给你几个可行的解决方案：

方案1：禁用AVX2优化编译 #

直接在f2py编译时添加-mno-avx2编译选项，让libgfortran使用非AVX2的matmul实现，避开有问题的优化路径。修改后的编译命令如下：

python -m numpy.f2py -c --f90flags='-O3 -mno-avx2' -m test test.f90

方案2：升级GCC编译器 #

这个断言错误是旧版本GCC（比如5.x、6.x系列）的libgfortran的已知bug，升级到GCC 7及以上版本后，该问题已经被修复。如果你的环境允许升级编译器，这是一劳永逸的解决方式。

方案3：手动实现矩阵乘法替代内置matmul #

如果不想修改编译选项或升级编译器，可以在Fortran代码中手动实现矩阵乘法循环，绕过libgfortran的内置matmul函数。修改后的test.f90代码如下：

!f2py dp selected_real_kind(15, 307)
subroutine wrap(A, B, C)
    real(kind=8), intent(in) :: A(:, :), B(:, :)
    real(kind=8), intent(inout) :: C(:, :)
    integer :: m, n, p, i, j, k

    m = size(A, 1)
    n = size(A, 2)
    p = size(B, 2)

    if (n /= size(B, 1)) then
        print *, "Error: Incompatible dimensions for matrix multiplication"
        return
    end if

    ! 手动实现矩阵乘法，替代内置matmul
    do i = 1, m
        do j = 1, p
            C(i,j) = 0.0d0
            do k = 1, n
                C(i,j) = C(i,j) + A(i,k) * B(k,j)
            end do
        end do
    end do

end subroutine wrap

之后用原来的编译命令重新编译即可。

备注：内容来源于stack exchange，提问作者Chungji