当前位置：首页 > news >正文

4-GGML:看一下加法算子!

news 2025/11/9 12:08:31

GGML

- CPU加法
- GPU加法

矩阵加法作为矩阵运算中比较简单的计算

M_d = M_0 + M_1

只需要把对应位置的元素加起来就好了,对于GGML的矩阵加法,其实有一点扩展的意味,他在一定程度上允许了两个矩阵的维度不是很一致也就是进行repeat的操作,如下定义的,这里注意一下顺序,其实在用的过程中他又反了一下.也就是

M_0

的维度是

M_1

的整数倍

// 要求 t1的每个维度都必须是t0的对应维度的整数倍
bool ggml_can_repeat(const struct ggml_tensor * t0, const struct ggml_tensor * t1) {static_assert(GGML_MAX_DIMS == 4, "GGML_MAX_DIMS is not 4 - update this function");return ggml_is_empty(t0) ? ggml_is_empty(t1) :(t1->ne[0]%t0->ne[0] == 0) &&(t1->ne[1]%t0->ne[1] == 0) &&(t1->ne[2]%t0->ne[2] == 0) &&(t1->ne[3]%t0->ne[3] == 0);
}

CPU加法

多线程的矩阵加法代码如下所示,对于关键的信息一步一步来看.

索引这个线程计算范围,也就是该线程计算多少个数从哪开始计算,矩阵计算中最小的单位是行,比如一共10个线程,矩阵一共100行,每个线程算10行就好了;如果矩阵只有1行,那么第一个线程算1行,其他线程休息.

const int ith = params->ith;// 线程编号
const int nth = params->nth;// 多线程总数
const int nr  = ggml_nrows(src0);// 行中变量的个数(tensor->ne[1]*tensor->ne[2]*tensor->ne[3])
// rows per thread 每个线程分配的行的个数
const int dr = (nr + nth - 1)/nth;
// row range for this thread 行的范围编号
const int ir0 = dr*ith;// 开始
const int ir1 = MIN(ir0 + dr, nr);// 结束,保证别超范围

对于src是连续的情况,连续的意思是,在0维度比如2*2的矩阵{{1,2},{3,4}}存储的时候连续存的,暂时还没遇到非连续的,遇到后会补充一下例子.这里只分析连续的情况吧,因为不连续的情况计算索引的原理是一致的.
计算在 $M_0$ 矩阵下的这一维度的索引,将所有的行号索引还原回各个维度上的索引,也就是这个1,2,3维度的索引值,同时你可以看到他是int64的单位

// 每个i03块包含ne02*ne01行
const int64_t i03 = ir/(ne02*ne01);                     // 这一行在i3上的编号
// 剩下的子小块再锁定ne01
const int64_t i02 = (ir - i03*ne02*ne01)/ne01;          // 这一行在i2上的编号
const int64_t i01 = (ir - i03*ne02*ne01 - i02*ne01);    // 这一行在i1上的编号

处理repeat的情况.
可以看i03是在0矩阵上的索引,对于1矩阵的维度可以是0矩阵整数倍缩小的.以下是计算矩阵1的行数据的索引.同时nr0代表的是如果行也是倍数关系是需要把行也复制多遍的

const int64_t i13 = i03 % ne13;
const int64_t i12 = i02 % ne12;
const int64_t i11 = i01 % ne11;
const int64_t nr0 = ne00 / ne10;// 这个是0维度的整数倍率

找到行地址.上面我们确定了计算的这一行位于矩阵各个维度的索引,用索引可以获得存储在内存空间的位置.注意这里用的是nb,关于nb不清楚的可以看解释tensor的那一个博文.

// 目标地址的行号 ptr
float * dst_ptr  = (float *) ((char *) dst->data  + i03*nb3  + i02*nb2  + i01*nb1 );
float * src0_ptr = (float *) ((char *) src0->data + i03*nb03 + i02*nb02 + i01*nb01);
float * src1_ptr = (float *) ((char *) src1->data + i13*nb13 + i12*nb12 + i11*nb11);

对每一行进行批量拷贝

for (int64_t r = 0; r < nr0; ++r) {// 拷贝倍率次// 个数 输出地址 输入0地址 输入1地址ggml_vec_add_f32(ne10, dst_ptr + r*ne10, src0_ptr + r*ne10, src1_ptr);
}

整体的代码可以过一遍了,梳理下流程:

分配任务:以行为单位,计算需要算哪行
将行号翻译成每一行在src0,src1上的矩阵维度
按唯独标号计算内存空间位置
批量加

// ggml前向推理中的add 32函数
static void ggml_compute_forward_add_f32(const struct ggml_compute_params * params,// 线程相关信息struct ggml_tensor * dst) {const struct ggml_tensor * src0 = dst->src[0];// add = 1+2 中的1const struct ggml_tensor * src1 = dst->src[1];// add = 1+2 中的2// src0与dst的维度要一样 && src0是src1的整数倍GGML_ASSERT(ggml_can_repeat(src1, src0) && ggml_are_same_shape(src0, dst));const int ith = params->ith;const int nth = params->nth;const int nr  = ggml_nrows(src0);// 行中变量的个数(tensor->ne[1]*tensor->ne[2]*tensor->ne[3])GGML_TENSOR_BINARY_OP_LOCALSGGML_ASSERT( nb0 == sizeof(float));GGML_ASSERT(nb00 == sizeof(float));// rows per thread 每个线程分配的行的个数const int dr = (nr + nth - 1)/nth;// row range for this thread 行的范围编号const int ir0 = dr*ith;const int ir1 = MIN(ir0 + dr, nr);if (nb10 == sizeof(float)) {// nb10是单个float大小 也就是src1是连续的 加一的最小单位是float 4个字节// 在所需要计算的行上编号上进行计算for (int ir = ir0; ir < ir1; ++ir) {// 对总行数进行区分// src1 is broadcastable across src0 and dst in i1, i2, i3// 每个i03块包含ne02*ne01行const int64_t i03 = ir/(ne02*ne01);                     // 这一行在i3上的编号// 剩下的子小块再锁定ne01const int64_t i02 = (ir - i03*ne02*ne01)/ne01;          // 这一行在i2上的编号const int64_t i01 = (ir - i03*ne02*ne01 - i02*ne01);    // 这一行在i1上的编号// 这一行在i3上的编号(考虑倍率) 如果ne13是等于ne03的，则肯定就按i03来计算// 其实这是为了处理src1和src0维度不一样的情况const int64_t i13 = i03 % ne13;const int64_t i12 = i02 % ne12;const int64_t i11 = i01 % ne11;const int64_t nr0 = ne00 / ne10;// 这个是0维度的整数倍率// 目标地址的行号 ptrfloat * dst_ptr  = (float *) ((char *) dst->data  + i03*nb3  + i02*nb2  + i01*nb1 );float * src0_ptr = (float *) ((char *) src0->data + i03*nb03 + i02*nb02 + i01*nb01);float * src1_ptr = (float *) ((char *) src1->data + i13*nb13 + i12*nb12 + i11*nb11);for (int64_t r = 0; r < nr0; ++r) {// 拷贝倍率次
#ifdef GGML_USE_ACCELERATEvDSP_vadd(src0_ptr + r*ne10, 1, src1_ptr, 1, dst_ptr + r*ne10, 1, ne10);
#else// 个数 输出地址 输入0地址 输入1地址ggml_vec_add_f32(ne10, dst_ptr + r*ne10, src0_ptr + r*ne10, src1_ptr);
#endif}}} else {// src1 is not contiguous 不连续的情况for (int ir = ir0; ir < ir1; ++ir) {// src1 is broadcastable across src0 and dst in i1, i2, i3const int64_t i03 = ir/(ne02*ne01);const int64_t i02 = (ir - i03*ne02*ne01)/ne01;const int64_t i01 = (ir - i03*ne02*ne01 - i02*ne01);const int64_t i13 = i03 % ne13;const int64_t i12 = i02 % ne12;const int64_t i11 = i01 % ne11;float * dst_ptr  = (float *) ((char *) dst->data  + i03*nb3  + i02*nb2  + i01*nb1 );float * src0_ptr = (float *) ((char *) src0->data + i03*nb03 + i02*nb02 + i01*nb01);for (int64_t i0 = 0; i0 < ne0; ++i0) {// 行的元素逐个进行计算const int64_t i10 = i0 % ne10;// 每次单独计算src1的地址float * src1_ptr = (float *) ((char *) src1->data + i13*nb13 + i12*nb12 + i11*nb11 + i10*nb10);dst_ptr[i0] = src0_ptr[i0] + *src1_ptr;}}}
}