我试着找出做矩阵乘法的最快方法,试了3种不同的方法:

  • 纯python实现:这里没有什么意外.
  • 使用numpy.dot(a, b)实现Numpy
  • 使用Python中的ctypes模块与C进行接口.

这是转换为共享库的C代码:

#include <stdio.h>
#include <stdlib.h>

void matmult(float* a, float* b, float* c, int n) {
    int i = 0;
    int j = 0;
    int k = 0;

    /*float* c = malloc(nay * sizeof(float));*/

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            int sub = 0;
            for (k = 0; k < n; k++) {
                sub = sub + a[i * n + k] * b[k * n + j];
            }
            c[i * n + j] = sub;
        }
    }
    return ;
}

以及调用它的Python代码:

def C_mat_mult(a, b):
    libmatmult = ctypes.CDLL("./matmult.so")

    dima = len(a) * len(a)
    dimb = len(b) * len(b)

    array_a = ctypes.c_float * dima
    array_b = ctypes.c_float * dimb
    array_c = ctypes.c_float * dima

    suma = array_a()
    sumb = array_b()
    sumc = array_c()

    inda = 0
    for i in range(0, len(a)):
        for j in range(0, len(a[i])):
            suma[inda] = a[i][j]
            inda = inda + 1
        indb = 0
    for i in range(0, len(b)):
        for j in range(0, len(b[i])):
            sumb[indb] = b[i][j]
            indb = indb + 1

    libmatmult.matmult(ctypes.byref(suma), ctypes.byref(sumb), ctypes.byref(sumc), 2);

    res = numpy.zeros([len(a), len(a)])
    indc = 0
    for i in range(0, len(sumc)):
        res[indc][i % len(a)] = sumc[i]
        if i % len(a) == len(a) - 1:
            indc = indc + 1

    return res

我敢打赌,使用C语言的版本会更快……我就输了!下面是我的基准,它似乎显示我要么做错了,要么numpy太快了:

benchmark

我想理解为什么numpy版本比ctypes版本快,我甚至不谈论纯Python实现,因为它有点明显.

推荐答案

我对Numpy不太熟悉,但源代码在Github上.dot产品的一部分是在https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/arraytypes.c.src中实现的,我假设它被转换 for each 数据类型的特定C实现.例如:

/**begin repeat
 *
 * #name = BYTE, UBYTE, SHORT, USHORT, INT, UINT,
 * LONG, ULONG, LONGLONG, ULONGLONG,
 * FLOAT, DOUBLE, LONGDOUBLE,
 * DATETIME, TIMEDELTA#
 * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint,
 * npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 * npy_float, npy_double, npy_longdouble,
 * npy_datetime, npy_timedelta#
 * #out = npy_long, npy_ulong, npy_long, npy_ulong, npy_long, npy_ulong,
 * npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 * npy_float, npy_double, npy_longdouble,
 * npy_datetime, npy_timedelta#
 */
static void
@name@_dot(char *ip1, npy_intp is1, char *ip2, npy_intp is2, char *op, npy_intp n,
           void *NPY_UNUSED(ignore))
{
    @out@ tmp = (@out@)0;
    npy_intp i;

    for (i = 0; i < n; i++, ip1 += is1, ip2 += is2) {
        tmp += (@out@)(*((@type@ *)ip1)) *
               (@out@)(*((@type@ *)ip2));
    }
    *((@type@ *)op) = (@type@) tmp;
}
/**end repeat**/

这似乎可以计算一维点积,即向量上的点积.在我浏览Github的几分钟时间里,我无法找到矩阵的源代码,但它可能会对结果矩阵中的每个元素调用一次FLOAT_dot.这意味着该函数中的循环对应于最内部的循环.

它们之间的一个区别是"跨步"——输入中连续元素之间的差异——在调用函数之前显式计算一次.在您的情况下,没有步幅,每次都会计算每个输入的偏移量,例如a[i * n + k].我本以为一个好的编译器会将其优化为类似于Numpy步幅的东西,但它可能无法证明步幅是一个常数(或者它没有被优化).

Numpy也可能在调用该函数的更高级别代码中使用了缓存效果.一个常见的技巧是考虑每一行是连续的,还是每一列都是连续的——然后try 先迭代每个连续的部分.似乎很难做到完美的最优,对于每个点积,一个输入矩阵必须按行遍历,另一个按列遍历(除非它们碰巧以不同的主顺序存储).但它至少可以为结果元素做到这一点.

Numpy还包含从不同的基本实现中 Select 特定操作(包括"点")实现的代码.例如,它可以使用BLAS个库.从上面的讨论来看,似乎使用了CBLA.这是从Fortran翻译成C的.我想你测试中使用的实现应该是在这里找到的:http://www.netlib.org/clapack/cblas/sdot.c.

请注意,该程序是由一台机器编写的,供另一台机器读取.但你可以在底部看到,它使用一个展开的循环一次处理5个元素:

for (i = mp1; i <= *n; i += 5) {
stemp = stemp + SX(i) * SY(i) + SX(i + 1) * SY(i + 1) + SX(i + 2) * 
    SY(i + 2) + SX(i + 3) * SY(i + 3) + SX(i + 4) * SY(i + 4);
}

这个展开因子很可能是在分析了几次之后 Select 的.但它的一个理论优势是在每个分支点之间进行更多的算术运算,编译器和CPU在如何优化调度它们以获得尽可能多的指令流水线方面有更多的 Select .

C++相关问答推荐

使用sd-设备列举设备导致seg错误

如何用C(使用两个s补数算术的32位程序)计算

为什么我一直收到分段错误?

如何创建一个C程序来存储5种动物的名字,并在用户 Select 其中任何一种动物时打印内存地址?

为什么cudaFree不需要数据 struct 的地址?

为什么我不能只在内存地址中添加一个int来寻址任何数组?

如何将字符**传递给需要常量字符指针的常量数组的函数

整型文字后缀在左移中的用途

在Rust和C之间使用ffi时如何通过 struct 中的[U8;1]成员传递指针

1处的解析器错误:yacc语法的语法错误

有什么方法可以将字符串与我们 Select 的子字符串分开吗?喜欢:SIN(LOG(10))

为什么用非常数指针变量改变常量静态变量时会出现分段错误?

如何读取程序中嵌入的数据S自己的ELF?

*S=0;正在优化中.可能是GCC 13号虫?或者是一些不明确的行为?

通过char*访问指针的对象表示是未定义的行为吗?

SSE 向量与 Epsilon 的比较

使用邻接表创建图

在 C23 之前如何对空指针使用nullptr?

初始化动态分配的布尔二维数组的最佳方法是什么?

我们可以在不违反标准的情况下向标准函数声明添加属性吗?