分块板子((N+Q)*(根号N))

本文介绍了一种使用分块思想优化区间更新和查询操作的方法。通过预处理数据和维护块内和块间的数据结构,实现对大规模数据集的有效操作。文章详细解释了算法的实现过程,包括初始化数据、处理区间更新和查询请求。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
int data[6000000],l[500000],r[500000],pos[500000],add[500000],n,m,t;
long long sum[200000];
int change(int l1,int r1,int d){
	int pl=pos[l1],pr=pos[r1];
	if(pl==pr){
		for(int i=l1;i<=r1;i++)data[i]++;
		sum[pl]+=d*(r1-l1+1);
	}else{
		for(int i=pl+1;i<=pr-1;i++)add[i]+=d;
		for(int i=l1;i<=r[pl];i++)data[i]+=d;
		sum[pl]+=(r[pl]-l1+1)*d;
		for(int i=l[pr];i<=r1;i++)data[i]+=d;
		sum[pr]+=(r1-l[pr]+1)*d;
	}//是小块朴素,大段维护的思想 
}
int ask(int l1,int r1){
	int pl=pos[l1],pr=pos[r1];
	long long res=0;
	if(pl==pr){
		for(int i=l1;i<=r1;i++)res+=(data[i]+add[pl]);		
	}else{
		for(int i=pl+1;i<=pr-1;i++)res+=(sum[i]+add[i]*(r[i]-l[i]+1));
		for(int i=l1;i<=r[pl];i++)res+=(data[i]+add[i]);
		for(int i=l[pr];i<=r1;i++)res+=(data[i]+add[i]);
	}
	return res;
}
int main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++)scanf("%d",&data[i]);
	//预处理左右断点 
	t=sqrt(n);
	for(int i=1;i<=t;i++){
		l[i]=(i-1)*sqrt(n)+1;
		r[i]=i*sqrt(n);
	}
	if(r[t]<n){
		t++;
		l[t]=r[t-1]+1;r[t]=n;
	}
	//预处理每个点在哪一段,并且求出每一段的和 
	for(int i=1;i<=t;i++){
		for(int j=l[i];j<=r[i];j++){
			pos[j]=i;
			sum[i]+=data[j];
		}
	} 
	char c;
	int x,y,d;
	while(m--){
		cin>>c;
		if(c=='C'){
			cin>>x>>y>>d;
			if(x>y)swap(x,y);
			change(x,y,d);
		}else{
			cin>>x>>y;
			cout<<ask(x,y);
		}
	}	
}

 

#include <immintrin.h> #include <stdio.h> // 用于调试时的printf // M, N, K 维度定义 // A: M x K, 列主序 // B: K x N, 行主序 // C: M x N, 行主序 /* 矩阵元素访问宏 */ // A_T 是列主序 #define A_T(i,j,lda,aptr) (aptr)[(j)*(lda) + (i)] // B 是行主序 #define B_RM(i,j,ldb,bptr) (bptr)[(i)*(ldb) + (j)] // C 是行主序 #define C_RM(i,j,ldc,cptr) (cptr)[(i)*(ldc) + (j)] // min 函数宏 #define min( i, j ) ( (i)<(j) ? (i): (j) ) // 分块大小 #ifndef MC #define MC 240 // M 维度的块大小 (A的行,C的行) 8的倍数 #endif #ifndef NC #define NC 768 // N 维度的块大小 (B的列,C的列) 24的倍数 #endif #ifndef KC #define KC 128 // K 维度的块大小 (A的列,B的行) #endif #define M_TILE 8 #define N_TILE 24 // AddDot static void AddDot8x24( int k_iter, const double *a_tile, int lda_full, /* M_dim for A */ const double *b_tile, int ldb_full, /* N_dim for B */ double *c_tile, int ldc_full ); /* N_dim for C */ // Kernel_McNcKc 函数:计算 C 的一个 MC x NC 子块 static void Kernel_McNcKc( int current_mc, int current_nc, int current_kc, const double *a_block, int global_lda, const double *b_block, int global_ldb, double *c_block, int global_ldc) { int i_block_loop, j_block_loop; // 每次处理一个 M_TILE x N_TILE 的tile for (i_block_loop = 0; i_block_loop < current_mc; i_block_loop += M_TILE) { for (j_block_loop = 0; j_block_loop < current_nc; j_block_loop += N_TILE) { const double *a_tile_ptr = &A_T(i_block_loop, 0, global_lda, a_block); const double *b_tile_ptr = &B_RM(0, j_block_loop, global_ldb, b_block); double *c_tile_ptr = &C_RM(i_block_loop, j_block_loop, global_ldc, c_block); AddDot8x24(current_kc, /*即K维度子块的大小 */ a_tile_ptr, global_lda, b_tile_ptr, global_ldb, c_tile_ptr, global_ldc); } } } void MY_MMult( int m_dim, int n_dim, int k_dim, const double *global_a_ptr, int lda_full, /* A的真实lda (等于M) */ const double *global_b_rm_ptr, int ldb_full, /* B的真实列数 (等于N) */ double *global_c_rm_ptr, int ldc_full ) /* C的真实列数 (等于N) */ { int ic_loop, pc_loop, jc_loop; int current_mc_val, current_nc_val, current_kc_val; //先n for (jc_loop = 0; jc_loop < n_dim; jc_loop += NC) { current_nc_val = min(n_dim - jc_loop, NC); //相当于没有 for (pc_loop = 0; pc_loop < k_dim; pc_loop += KC) { current_kc_val = min(k_dim - pc_loop, KC); const double *a_panel_for_kc_block = &A_T(0, pc_loop, lda_full, global_a_ptr); const double *b_panel_for_kc_block = &B_RM(pc_loop, 0, ldb_full, global_b_rm_ptr); for (ic_loop = 0; ic_loop < m_dim; ic_loop += MC) { current_mc_val = min(m_dim - ic_loop, MC); // A: A(ic_loop : ic_loop+current_mc_val-1, pc_loop : pc_loop+current_kc_val-1) // B: B(pc_loop : pc_loop+current_kc_val-1, jc_loop : jc_loop+current_nc_val-1) // C: C(ic_loop : ic_loop+current_mc_val-1, jc_loop : jc_loop+current_nc_val-1) // Calculate start pointers for the panel/sub-block for Kernel_McNcKc // A panel starts at global A [ic_loop, pc_loop] const double *a_panel_ptr = &A_T(ic_loop, 0, lda_full, a_panel_for_kc_block); // B panel starts at global B [pc_loop, jc_loop] const double *b_panel_ptr = &B_RM(0, jc_loop, ldb_full, b_panel_for_kc_block); // C sub-block starts at global C [ic_loop, jc_loop] double *c_sub_tile_ptr = &C_RM(ic_loop, jc_loop, ldc_full, global_c_rm_ptr); Kernel_McNcKc(current_mc_val, current_nc_val, current_kc_val, a_panel_ptr, lda_full, b_panel_ptr, ldb_full, c_sub_tile_ptr, ldc_full); } } } } static void AddDot8x24( int k_iter, const double *a_tile, int lda_full, const double *b_tile, int ldb_full, double *c_tile, int ldc_full ) { __m512d c_reg_00, c_reg_01, c_reg_02; __m512d c_reg_10, c_reg_11, c_reg_12; __m512d c_reg_20, c_reg_21, c_reg_22; __m512d c_reg_30, c_reg_31, c_reg_32; __m512d c_reg_40, c_reg_41, c_reg_42; __m512d c_reg_50, c_reg_51, c_reg_52; __m512d c_reg_60, c_reg_61, c_reg_62; __m512d c_reg_70, c_reg_71, c_reg_72; c_reg_00 = _mm512_load_pd(c_tile + 0*ldc_full + 0*8); c_reg_01 = _mm512_load_pd(c_tile + 0*ldc_full + 1*8); c_reg_02 = _mm512_load_pd(c_tile + 0*ldc_full + 2*8); c_reg_10 = _mm512_load_pd(c_tile + 1*ldc_full + 0*8); c_reg_11 = _mm512_load_pd(c_tile + 1*ldc_full + 1*8); c_reg_12 = _mm512_load_pd(c_tile + 1*ldc_full + 2*8); c_reg_20 = _mm512_load_pd(c_tile + 2*ldc_full + 0*8); c_reg_21 = _mm512_load_pd(c_tile + 2*ldc_full + 1*8); c_reg_22 = _mm512_load_pd(c_tile + 2*ldc_full + 2*8); c_reg_30 = _mm512_load_pd(c_tile + 3*ldc_full + 0*8); c_reg_31 = _mm512_load_pd(c_tile + 3*ldc_full + 1*8); c_reg_32 = _mm512_load_pd(c_tile + 3*ldc_full + 2*8); c_reg_40 = _mm512_load_pd(c_tile + 4*ldc_full + 0*8); c_reg_41 = _mm512_load_pd(c_tile + 4*ldc_full + 1*8); c_reg_42 = _mm512_load_pd(c_tile + 4*ldc_full + 2*8); c_reg_50 = _mm512_load_pd(c_tile + 5*ldc_full + 0*8); c_reg_51 = _mm512_load_pd(c_tile + 5*ldc_full + 1*8); c_reg_52 = _mm512_load_pd(c_tile + 5*ldc_full + 2*8); c_reg_60 = _mm512_load_pd(c_tile + 6*ldc_full + 0*8); c_reg_61 = _mm512_load_pd(c_tile + 6*ldc_full + 1*8); c_reg_62 = _mm512_load_pd(c_tile + 6*ldc_full + 2*8); c_reg_70 = _mm512_load_pd(c_tile + 7*ldc_full + 0*8); c_reg_71 = _mm512_load_pd(c_tile + 7*ldc_full + 1*8); c_reg_72 = _mm512_load_pd(c_tile + 7*ldc_full + 2*8); __m512d a_bcast_0, a_bcast_1, a_bcast_2, a_bcast_3; __m512d b_vec_0, b_vec_1, b_vec_2; const double *a_k_ptr = a_tile; const double *b_k_ptr = b_tile; for (int p = 0; p < k_iter; ++p) { b_vec_0 = _mm512_load_pd(b_k_ptr + 0*8); b_vec_1 = _mm512_load_pd(b_k_ptr + 1*8); b_vec_2 = _mm512_load_pd(b_k_ptr + 2*8); a_bcast_0 = _mm512_set1_pd(*(a_k_ptr + 0)); a_bcast_1 = _mm512_set1_pd(*(a_k_ptr + 1)); a_bcast_2 = _mm512_set1_pd(*(a_k_ptr + 2)); a_bcast_3 = _mm512_set1_pd(*(a_k_ptr + 3)); c_reg_00 = _mm512_fmadd_pd(a_bcast_0, b_vec_0, c_reg_00); c_reg_01 = _mm512_fmadd_pd(a_bcast_0, b_vec_1, c_reg_01); c_reg_02 = _mm512_fmadd_pd(a_bcast_0, b_vec_2, c_reg_02); c_reg_10 = _mm512_fmadd_pd(a_bcast_1, b_vec_0, c_reg_10); c_reg_11 = _mm512_fmadd_pd(a_bcast_1, b_vec_1, c_reg_11); c_reg_12 = _mm512_fmadd_pd(a_bcast_1, b_vec_2, c_reg_12); c_reg_20 = _mm512_fmadd_pd(a_bcast_2, b_vec_0, c_reg_20); c_reg_21 = _mm512_fmadd_pd(a_bcast_2, b_vec_1, c_reg_21); c_reg_22 = _mm512_fmadd_pd(a_bcast_2, b_vec_2, c_reg_22); c_reg_30 = _mm512_fmadd_pd(a_bcast_3, b_vec_0, c_reg_30); c_reg_31 = _mm512_fmadd_pd(a_bcast_3, b_vec_1, c_reg_31); c_reg_32 = _mm512_fmadd_pd(a_bcast_3, b_vec_2, c_reg_32); a_bcast_0 = _mm512_set1_pd(*(a_k_ptr + 4)); a_bcast_1 = _mm512_set1_pd(*(a_k_ptr + 5)); a_bcast_2 = _mm512_set1_pd(*(a_k_ptr + 6)); a_bcast_3 = _mm512_set1_pd(*(a_k_ptr + 7)); c_reg_40 = _mm512_fmadd_pd(a_bcast_0, b_vec_0, c_reg_40); c_reg_41 = _mm512_fmadd_pd(a_bcast_0, b_vec_1, c_reg_41); c_reg_42 = _mm512_fmadd_pd(a_bcast_0, b_vec_2, c_reg_42); c_reg_50 = _mm512_fmadd_pd(a_bcast_1, b_vec_0, c_reg_50); c_reg_51 = _mm512_fmadd_pd(a_bcast_1, b_vec_1, c_reg_51); c_reg_52 = _mm512_fmadd_pd(a_bcast_1, b_vec_2, c_reg_52); c_reg_60 = _mm512_fmadd_pd(a_bcast_2, b_vec_0, c_reg_60); c_reg_61 = _mm512_fmadd_pd(a_bcast_2, b_vec_1, c_reg_61); c_reg_62 = _mm512_fmadd_pd(a_bcast_2, b_vec_2, c_reg_62); c_reg_70 = _mm512_fmadd_pd(a_bcast_3, b_vec_0, c_reg_70); c_reg_71 = _mm512_fmadd_pd(a_bcast_3, b_vec_1, c_reg_71); c_reg_72 = _mm512_fmadd_pd(a_bcast_3, b_vec_2, c_reg_72); a_k_ptr += lda_full; b_k_ptr += ldb_full; } _mm512_store_pd(c_tile + 0*ldc_full + 0*8, c_reg_00); _mm512_store_pd(c_tile + 0*ldc_full + 1*8, c_reg_01); _mm512_store_pd(c_tile + 0*ldc_full + 2*8, c_reg_02); _mm512_store_pd(c_tile + 1*ldc_full + 0*8, c_reg_10); _mm512_store_pd(c_tile + 1*ldc_full + 1*8, c_reg_11); _mm512_store_pd(c_tile + 1*ldc_full + 2*8, c_reg_12); _mm512_store_pd(c_tile + 2*ldc_full + 0*8, c_reg_20); _mm512_store_pd(c_tile + 2*ldc_full + 1*8, c_reg_21); _mm512_store_pd(c_tile + 2*ldc_full + 2*8, c_reg_22); _mm512_store_pd(c_tile + 3*ldc_full + 0*8, c_reg_30); _mm512_store_pd(c_tile + 3*ldc_full + 1*8, c_reg_31); _mm512_store_pd(c_tile + 3*ldc_full + 2*8, c_reg_32); _mm512_store_pd(c_tile + 4*ldc_full + 0*8, c_reg_40); _mm512_store_pd(c_tile + 4*ldc_full + 1*8, c_reg_41); _mm512_store_pd(c_tile + 4*ldc_full + 2*8, c_reg_42); _mm512_store_pd(c_tile + 5*ldc_full + 0*8, c_reg_50); _mm512_store_pd(c_tile + 5*ldc_full + 1*8, c_reg_51); _mm512_store_pd(c_tile + 5*ldc_full + 2*8, c_reg_52); _mm512_store_pd(c_tile + 6*ldc_full + 0*8, c_reg_60); _mm512_store_pd(c_tile + 6*ldc_full + 1*8, c_reg_61); _mm512_store_pd(c_tile + 6*ldc_full + 2*8, c_reg_62); _mm512_store_pd(c_tile + 7*ldc_full + 0*8, c_reg_70); _mm512_store_pd(c_tile + 7*ldc_full + 1*8, c_reg_71); _mm512_store_pd(c_tile + 7*ldc_full + 2*8, c_reg_72); } 请告诉我这段代码完成了什么优化工作,请详细说明
最新发布
05-29
<think>好的,我现在需要分析用户提供的这段C代码,了解它完成了什么优化工作。用户提到这是一个矩阵乘法的实现,使用了分块技术和AVX-512指令集。我需要详细解释代码中的优化策略,并指出各个部分的作用。 首先,我会通读代码的结构,了解各个函数的作用。主函数是MY_MMult,它负责处理矩阵的分块,然后调用Kernel_McNcKc来处理每个分块后的子矩阵。Kernel_McNcKc进一步将子矩阵划分为更小的块(8x24),并调用AddDot8x24进行计算。AddDot8x24函数使用了AVX-512指令进行向量化计算,这是关键的性能优化部分。接下来,我需要确定分块的具体参数。代码中定义了MC、NC、KC三个宏,分别代表M、N、K维度上的分块大小。MC为240,NC为768,KC为128。这些值的选择通常是为了适应CPU的缓存结构,减少缓存未命中,提高数据局部性。例如,较大的KC值可能有助于在计算时重复利用A和B的子块,减少内存访问次数。 然后,Kernel_McNcKc中的循环结构是将当前的分块(MCxNC)进一步划分为8x24的小块。这可能是为了适配AVX-512的寄存器大小,每个AVX-512寄存器可以处理8个双精度浮点数(512位)。而24这个数值可能涉及到循环展开的策略,使得每次迭代能处理更多的数据,减少循环开销。在AddDot8x24函数中,使用了大量的__m512d类型的变量,这些是AVX-512的向量寄存器,每个可以存储8个双精度浮点数。代码中加载和存储数据使用了_mm512_load_pd和_mm512_store_pd,这些指令能够高效地读写内存。计算部分使用了_mm512_fmadd_pd,这是融合乘加指令,能够在单个指令中完成乘法和加法操作,提高计算吞吐量。观察循环结构,AddDot8x24中的循环处理k_iter次,每次迭代从A和B中加载数据,进行广播和向量乘加。这里将A的8行数据分成两组,每组4行,分别与B的三列向量进行乘加。这样的设计可能充分利用了寄存器的数量,避免寄存器溢出,同时提高指令级并行度。 另外,代码中的矩阵存储顺序需要注意。A矩阵是列主序,而B和C是行主序。这可能影响分块时的数据访问模式,列主序的A在分块时需要按列访问,而行主序的B和C则按行访问。分块策略可能优化了不同存储顺序下的内存访问,减少缓存行的冲突。整个优化过程综合应用了多种技术:分块(Blocking)以减少缓存未命中,循环展开以提高指令流水线效率,向量化(AVX-512)利用SIMD指令并行计算,以及寄存器重用减少内存访问。这些优化手段共同提升了矩阵乘法的性能,尤其是在大规模矩阵运算时效果显著。需要进一步确认的是,分块大小是否适配目标CPU的缓存大小。例如,KC=128可能使得A的子块(MCxKC)和B的子块(KCxNC)能够放入L2或L3缓存中,减少重复访问主存的次数。此外,8x24的小块处理可能确保每次计算的数据都能保存在寄存器中,减少临时数据的存储开销。 总结来说,这段代码通过分块技术优化了缓存利用率,利用AVX-512指令进行向量化计算,提高了并行处理能力,同时通过循环展开和寄存器重用减少了指令开销和内存访问延迟,从而实现了高效的矩阵乘法运算。</think>这段代码实现了一个高度优化的矩阵乘法函数,主要采用了分块技术(Blocking)和AVX-512向量化指令集优化。以下是关键优化点的详细说明: --- ### 1. **三级分块策略** - **MC (240)**: M维度分块大小(A的行数/C的行数) - **NC (768)**: N维度分块大小(B的列数/C的列数) - **KC (128)**: K维度分块大小(A的列数/B的行数) | 分块类型 | 作用 | |---------|------| | 外层分块 | 将矩阵划分为适合L3缓存的块,减少内存带宽压力 | | 中层分块 | 优化L2缓存利用率,提升数据重用率 | | 内层分块 | 通过8x24的微内核适配AVX-512寄存器 | --- ### 2. **AVX-512向量化内核** 在`AddDot8x24`函数中使用了以下优化: ```c // 关键向量化代码示例 __m512d c_reg_00 = _mm512_load_pd(c_tile); // 加载8个double __m512d b_vec_0 = _mm512_load_pd(b_k_ptr); c_reg_00 = _mm512_fmadd_pd(a_bcast_0, b_vec_0, c_reg_00); // 融合乘加 ``` - **寄存器重用**:同时保持8x24个中间结果寄存器 - **数据广播**:`_mm512_set1_pd`实现标量到向量的广播 - **内存对齐**:`_mm512_load_pd`要求64字节对齐 --- ### 3. **内存访问优化** - **A矩阵列主序**:通过`A_T(i,j)`宏实现列优先访问 - **B/C矩阵行主序**:通过`B_RM/C_RM`宏实现行优先访问 - **缓存友好访问模式**:分块后数据呈现局部性特征 --- ### 4. **微内核设计** ```c #define M_TILE 8 // 对应AVX-512的8个double #define N_TILE 24 // 3组8列向量 ``` - **计算密度**:每次处理8行x24列的矩阵块 - **流水线优化**:通过循环展开隐藏指令延迟 --- ### 5. **性能提升效果** | 优化手段 | 理论加速比 | |---------|-----------| | AVX-512向量化 | 8x (双精度) | | 分块缓存优化 | 2-5x | | 指令级并行 | 1.5-3x | --- ### 6. **代码执行流程** 1. **外层循环**:沿N维度分块(`jc_loop`) 2. **中层循环**:沿K维度分块(`pc_loop`) 3. **内层循环**:沿M维度分块(`ic_loop`) 4. **微内核计算**:8x24块的高强度向量化计算 ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值