优化C++矩阵乘法性能
我将设计一个性能优化的矩阵乘法实现，包括基础实现和多种优化技术。让我们从基础开始，逐步添加优化。
基础实现
，我们实现一个基础的矩阵乘法，然后逐步优化：
#include <iostream>
#include <vector>
#include <chrono>
#include <cstring>
#include <immintrin.h> // AVX指令集
#include <omp.h> // OpenMP并行
// 基础矩阵乘法
void matrix_multiply_basic(const std::vector<float>& A,
const std::vector<float>& B,
std::vector<float>& C,
int M, int N, int K) {
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; ++j) {
float sum = 0.0f;
for (int k = 0; k < K; ++k) {
sum += A[i * K + k] * B[k * N + j];
}
C[i * N + j] = sum;
}
}
}
// 优化1：循环重排序（提高缓存局部性）
void matrix_multiply_loop_reorder(const std::vector<float>& A,
const std::vector<float>& B,
std::vector<float>& C,
int M, int N, int K) {
// 将B矩阵转置以提高缓存命中率
std::vector<float> B_transpose(N * K);
for (int i = 0; i < K; ++i) {
for (int j = 0; j < N; ++j) {
B_transpose[j * K + i] = B[i * N + j];
}
}
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; ++j) {
float sum = 0.0f;
const float* A_row = &A[i * K];
const float* B_col = &B_transpose[j * K];
for (int k = 0; k < K; ++k) {
sum += A_row[k] * B_col[k];
}
C[i * N + j] = sum;
}
}
}
// 优化2：循环分块（Blocking）
void matrix_multiply_blocking(const std::vector<float>& A,
const std::vector<float>& B,
std::vector<float>& C,
int M, int N, int K,
int block_size = 64) {
// 清零C矩阵
std::memset(C.data(), 0, M * N * sizeof(float));
for (int i = 0; i < M; i += block_size) {
for (int j = 0; j < N; j += block_size) {
for (int k = 0; k < K; k += block_size) {
// 计算当前块的大小
int i_end = std::min(i + block_size, M);
int j_end = std::min(j + block_size, N);
int k_end = std::min(k + block_size, K);
// 处理当前块
for (int ii = i; ii < i_end; ++ii) {
for (int kk = k; kk < k_end; ++kk) {
float a_val = A[ii * K + kk];
for (int jj = j; jj < j_end; ++jj) {
C[ii * N + jj] += a_val * B[kk * N + jj];
}
}
}
}
}
}
}
// 优化3：使用AVX指令集（SIMD）
#ifdef __AVX__
void matrix_multiply_avx(const std::vector<float>& A,
const std::vector<float>& B,
std::vector<float>& C,
int M, int N, int K) {
// 将B矩阵转置
std::vector<float> B_transpose(N * K);
for (int i = 0; i < K; ++i) {
for (int j = 0; j < N; ++j) {
B_transpose[j * K + i] = B[i * N + j];
}
}
// 清零C矩阵
std::memset(C.data(), 0, M * N * sizeof(float));
// AVX一次处理8个float
const int simd_width = 8;
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; j += simd_width) {
// 加载C矩阵的值
__m256 c_vec = _mm256_loadu_ps(&C[i * N + j]);
for (int k = 0; k < K; ++k) {
// 广播A[i][k]
__m256 a_vec = _mm256_set1_ps(A[i * K + k]);
// 加载B的一行
__m256 b_vec = _mm256_loadu_ps(&B_transpose[j * K + k]);
// 乘积累加
c_vec = _mm256_fmadd_ps(a_vec, b_vec, c_vec);
}
// 存储结果
_mm256_storeu_ps(&C[i * N + j], c_vec);
}
}
}
#endif
// 优化4：OpenMP多线程并行
void matrix_multiply_parallel(const std::vector<float>& A,
const std::vector<float>& B,
std::vector<float>& C,
int M, int N, int K) {
#pragma omp parallel for
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; ++j) {
float sum = 0.0f;
for (int k = 0; k < K; ++k) {
sum += A[i * K + k] * B[k * N + j];
}
C[i * N + j] = sum;
}
}
}
// 优化5：综合优化（分块+AVX+OpenMP）
void matrix_multiply_optimized(const std::vector<float>& A,
const std::vector<float>& B,
std::vector<float>& C,
int M, int N, int K,
int block_size = 64) {
// 清零C矩阵
std::memset(C.data(), 0, M * N * sizeof(float));
#pragma omp parallel for
for (int i = 0; i < M; i += block_size) {
for (int j = 0; j < N; j += block_size) {
for (int k = 0; k < K; k += block_size) {
// 计算当前块的大小
int i_end = std::min(i + block_size, M);
int j_end = std::min(j + block_size, N);
int k_end = std::min(k + block_size, K);
// 处理当前块
for (int ii = i; ii < i_end; ++ii) {
for (int kk = k; kk < k_end; ++kk) {
float a_val = A[ii * K + kk];
// 使用AVX优化内层循环
#ifdef __AVX__
int jj = j;
__m256 a_vec = _mm256_set1_ps(a_val);
// AVX处理（每次8个元素）
for (; jj + 7 < j_end; jj += 8) {
__m256 b_vec = _mm256_loadu_ps(&B[kk * N + jj]);
__m256 c_vec = _mm256_loadu_ps(&C[ii * N + jj]);
c_vec = _mm256_fmadd_ps(a_vec, b_vec, c_vec);
_mm256_storeu_ps(&C[ii * N + jj], c_vec);
}
// 处理剩余元素
for (; jj < j_end; ++jj) {
C[ii * N + jj] += a_val * B[kk * N + jj];
}
#else
// 非AVX版本
for (int jj = j; jj < j_end; ++jj) {
C[ii * N + jj] += a_val * B[kk * N + jj];
}
#endif
}
}
}
}
}
}
// 性能测试函数
template<typename Func>
double benchmark(Func func, const std::vector<float>& A,
const std::vector<float>& B, std::vector<float>& C,
int M, int N, int K, int iterations = 5) {
double total_time = 0.0;
for (int iter = 0; iter < iterations; ++iter) {
auto start = std::chrono::high_resolution_clock::now();
func(A, B, C, M, N, K);
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> duration = end - start;
total_time += duration.count();
}
return total_time / iterations;
}
// 验证结果正确性
bool verify_results(const std::vector<float>& C1,
const std::vector<float>& C2,
int size, float epsilon = 1e-5) {
for (int i = 0; i < size; ++i) {
if (std::abs(C1[i] - C2[i]) > epsilon) {
std::cout << "Mismatch at index " << i
<< ": " << C1[i] << " vs " << C2[i] << std::endl;
return false;
}
}
return true;
}
int main() {
// 设置矩阵大小
const int M = 512;
const int N = 512;
const int K = 512;
// 分配内存
std::vector<float> A(M * K);
std::vector<float> B(K * N);
std::vector<float> C1(M * N);
std::vector<float> C2(M * N);
std::vector<float> C3(M * N);
std::vector<float> C4(M * N);
std::vector<float> C5(M * N);
std::vector<float> C6(M * N);
// 初始化矩阵
for (int i = 0; i < M * K; ++i) A[i] = static_cast<float>(rand()) / RAND_MAX;
for (int i = 0; i < K * N; ++i) B[i] = static_cast<float>(rand()) / RAND_MAX;
std::cout << "矩阵大小: " << M << "x" << K << " * " << K << "x" << N << std::endl;
std::cout << "==========================================" << std::endl;
// 测试各种实现
double time_basic = benchmark(matrix_multiply_basic, A, B, C1, M, N, K);
std::cout << "基础实现: " << time_basic << " 秒" << std::endl;
double time_reorder = benchmark(matrix_multiply_loop_reorder, A, B, C2, M, N, K);
std::cout << "循环重排序: " << time_reorder << " 秒 (加速比: "
<< time_basic / time_reorder << "x)" << std::endl;
double time_blocking = benchmark(matrix_multiply_blocking, A, B, C3, M, N, K);
std::cout << "循环分块: " << time_blocking << " 秒 (加速比: "
<< time_basic / time_blocking << "x)" << std::endl;
#ifdef __AVX__
double time_avx = benchmark(matrix_multiply_avx, A, B, C4, M, N, K);
std::cout << "AVX指令集: " << time_avx << " 秒 (加速比: "
<< time_basic / time_avx << "x)" << std::endl;
#endif
double time_parallel = benchmark(matrix_multiply_parallel, A, B, C5, M, N, K);
std::cout << "OpenMP并行: " << time_parallel << " 秒 (加速比: "
<< time_basic / time_parallel << "x)" << std::endl;
double time_optimized = benchmark(matrix_multiply_optimized, A, B, C6, M, N, K);
std::cout << "综合优化: " << time_optimized << " 秒 (加速比: "
<< time_basic / time_optimized << "x)" << std::endl;
// 验证结果正确性
std::cout << "\n验证结果正确性:" << std::endl;
std::cout << "基础 vs 循环重排序: "
<< (verify_results(C1, C2, M * N) ? "通过" : "失败") << std::endl;
std::cout << "基础 vs 循环分块: "
<< (verify_results(C1, C3, M * N) ? "通过" : "失败") << std::endl;
#ifdef __AVX__
std::cout << "基础 vs AVX: "
<< (verify_results(C1, C4, M * N) ? "通过" : "失败") << std::endl;
#endif
std::cout << "基础 vs 综合优化: "
<< (verify_results(C1, C6, M * N) ? "通过" : "失败") << std::endl;
// 计算FLOPS
double flops = 2.0 * M * N * K;
std::cout << "\n性能指标 (GFLOPS):" << std::endl;
std::cout << "基础实现: " << (flops / time_basic) / 1e9 << " GFLOPS" << std::endl;
std::cout << "综合优化: " << (flops / time_optimized) / 1e9 << " GFLOPS" << std::endl;
return 0;
}编译和运行说明
编译选项
# 基础编译（无优化）
g++ -O2 -o matrix_multiply matrix_multiply.cpp
# 启用AVX指令集
g++ -O2 -mavx -mfma -o matrix_multiply_avx matrix_multiply.cpp
# 启用OpenMP和AVX
g++ -O2 -mavx -mfma -fopenmp -o matrix_multiply_opt matrix_multiply.cpp
# 最高优化级别
g++ -O3 -march=native -fopenmp -o matrix_multiply_fast matrix_multiply.cpp运行程序
./matrix_multiply_opt优化技术总结
-
循环重排序：通过转置B矩阵，将内存访问模式从列优先改为行优先，提高缓存命中率。
-
循环分块（Blocking）：将大矩阵分解为小块，使数据在缓存中停留更久，减少缓存失效。
-
SIMD指令集（AVX）：使用单指令多数据技术，一次处理多个浮点数运算。
-
多线程并行（OpenMP）：利用多核CPU并行计算，提高吞吐量。
-
综合优化：结合分块、SIMD和并行化，实现性能。
预期性能提升
在典型硬件上（如Intel Core i7/i9处理器）：
- 基础实现：~1-2 GFLOPS
- 循环重排序：~3-5 GFLOPS（2-3倍加速）
- 循环分块：~10-20 GFLOPS（5-10倍加速）
- AVX指令集：~30-50 GFLOPS（15-25倍加速）
- OpenMP并行：~100-200 GFLOPS（50-100倍加速）
- 综合优化：~200-400 GFLOPS（100-200倍加速）
注意事项
- 矩阵大小应为块大小的整数倍，以获得性能
- 实际性能取决于CPU架构、缓存大小和内存带宽
- 对于非常大的矩阵，可能需要考虑内存分页和NUMA优化
- 某些优化可能需要特定的CPU指令集支持
这个实现展示了从基础到高级的矩阵乘法优化技术，可以根据具体硬件和需求选择合适的优化组合。
作者声明:本文包含人工智能生成内容。
