GPUKernelContest/S1/3/reduce_sum_algorithm.maca

598 lines
21 KiB
Plaintext
Executable File
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#include "test_utils.h"
#include "performance_utils.h"
#include "yaml_reporter.h"
#include <iostream>
#include <vector>
#include <iomanip>
#include <stdint.h>
#include <cmath>
// ============================================================================
// 实现标记宏 - 参赛者修改实现时请将此宏设为0
// ============================================================================
#ifndef USE_DEFAULT_REF_IMPL
#define USE_DEFAULT_REF_IMPL 0 // 1=默认实现, 0=参赛者自定义实现
#endif
#if USE_DEFAULT_REF_IMPL
#include <thrust/reduce.h>
#include <thrust/device_vector.h>
#include <thrust/execution_policy.h>
#include <thrust/functional.h>
#else
// ========================================
// 参赛者自定义实现区域
// ========================================
/**
* @brief 使用warp shuffle指令进行warp内归约
*/
__device__ __forceinline__ double warpShuffle(double value) {
#pragma unroll
for (int delta = 16; delta > 0; delta >>= 1) {
value += __shfl_down_sync(0xffffffff, value, delta);
}
return value;
}
/**
* @brief 优化的块级归约基于warp shuffle
*/
__device__ __forceinline__ double blockWiseReduce(double localSum, double* sharedMem) {
const int tid = threadIdx.x;
const int laneId = tid & 0x1F; // tid % 32
const int warpId = tid >> 5; // tid / 32
// Phase 1: warp内归约
localSum = warpShuffle(localSum);
// Phase 2: warp leader写入共享内存
if (laneId == 0) {
sharedMem[warpId] = localSum;
}
__syncthreads();
// Phase 3: 第一个warp处理所有warp的结果
if (warpId == 0) {
localSum = (tid < (blockDim.x >> 5)) ? sharedMem[tid] : 0.0;
localSum = warpShuffle(localSum);
}
return localSum;
}
// ============================================================================
// 策略1: 小数据量优化 (num_items <= 1000000)
// 使用grid-stride循环 + 高展开因子
// ============================================================================
/**
* @brief 第一阶段归约 - 小数据量专用
*/
__global__ void firstStageReduce_Small(const float* __restrict__ input,
double* __restrict__ partialSums,
int totalElements) {
double threadSum = 0.0;
const int tid = threadIdx.x;
const int blockId = blockIdx.x;
const int threadsPerBlock = blockDim.x;
const int totalBlocks = gridDim.x;
const int totalThreads = totalBlocks * threadsPerBlock;
extern __shared__ double shmem[];
// 8路循环展开
const int UNROLL = 8;
int globalIdx = blockId * threadsPerBlock + tid;
const int stride = totalThreads * UNROLL;
// 主循环:处理对齐的数据
for (; globalIdx + (totalThreads * (UNROLL - 1)) < totalElements; globalIdx += stride) {
threadSum += static_cast<double>(input[globalIdx]);
threadSum += static_cast<double>(input[globalIdx + totalThreads]);
threadSum += static_cast<double>(input[globalIdx + totalThreads * 2]);
threadSum += static_cast<double>(input[globalIdx + totalThreads * 3]);
threadSum += static_cast<double>(input[globalIdx + totalThreads * 4]);
threadSum += static_cast<double>(input[globalIdx + totalThreads * 5]);
threadSum += static_cast<double>(input[globalIdx + totalThreads * 6]);
threadSum += static_cast<double>(input[globalIdx + totalThreads * 7]);
}
// 尾部处理
for (; globalIdx < totalElements; globalIdx += totalThreads) {
threadSum += static_cast<double>(input[globalIdx]);
}
// 块内归约
threadSum = blockWiseReduce(threadSum, shmem);
// 块leader写出结果
if (tid == 0) {
partialSums[blockId] = threadSum;
}
}
/**
* @brief 第二阶段归约 - 小数据量专用
*/
__global__ void secondStageReduce_Small(const double* __restrict__ partialSums,
float* __restrict__ finalResult,
int numPartials,
float initVal) {
double threadSum = 0.0;
const int tid = threadIdx.x;
const int blockSize = blockDim.x;
extern __shared__ double shmem[];
// Grid-stride模式收集所有部分和
for (int idx = tid; idx < numPartials; idx += blockSize) {
threadSum += partialSums[idx];
}
// 块内归约
threadSum = blockWiseReduce(threadSum, shmem);
// 写出最终结果
if (tid == 0) {
threadSum += static_cast<double>(initVal);
*finalResult = static_cast<float>(threadSum);
}
}
// ============================================================================
// 策略2: 大数据量优化 (num_items > 1000000)
// 每个线程处理多个元素,使用更激进的展开
// ============================================================================
/**
* @brief 第一阶段归约 - 大数据量专用4路展开
*/
template <typename T>
__global__ void firstStageReduce_Large(const T* __restrict__ input,
T* __restrict__ partialSums,
int totalElements) {
extern __shared__ T localCache[];
const int tid = threadIdx.x;
const int blockSize = blockDim.x;
const int baseIdx = blockIdx.x * blockSize * 4 + tid; // 每个线程处理4个元素
T accumulator = static_cast<T>(0);
// 4路展开加载和累加
if (baseIdx < totalElements) {
accumulator += input[baseIdx];
}
if (baseIdx + blockSize < totalElements) {
accumulator += input[baseIdx + blockSize];
}
if (baseIdx + blockSize * 2 < totalElements) {
accumulator += input[baseIdx + blockSize * 2];
}
if (baseIdx + blockSize * 3 < totalElements) {
accumulator += input[baseIdx + blockSize * 3];
}
localCache[tid] = accumulator;
__syncthreads();
// 树状归约 - 分两段处理
// 第一段普通归约到warp级别
for (int offset = blockSize / 2; offset > 32; offset >>= 1) {
if (tid < offset) {
localCache[tid] += localCache[tid + offset];
}
__syncthreads();
}
// 第二段warp内无同步归约
if (tid < 32) {
volatile T* volatileCache = localCache;
if (blockSize >= 64) volatileCache[tid] += volatileCache[tid + 32];
if (blockSize >= 32) volatileCache[tid] += volatileCache[tid + 16];
if (blockSize >= 16) volatileCache[tid] += volatileCache[tid + 8];
if (blockSize >= 8) volatileCache[tid] += volatileCache[tid + 4];
if (blockSize >= 4) volatileCache[tid] += volatileCache[tid + 2];
if (blockSize >= 2) volatileCache[tid] += volatileCache[tid + 1];
}
// 块leader写出部分和
if (tid == 0) {
partialSums[blockIdx.x] = localCache[0];
}
}
/**
* @brief 第二阶段归约 - 大数据量专用
*/
template <typename T>
__global__ void secondStageReduce_Large(const T* __restrict__ partialSums,
T* __restrict__ finalResult,
int numPartials,
T initVal) {
extern __shared__ T localCache[];
const int tid = threadIdx.x;
const int blockSize = blockDim.x;
// 每个线程收集一个或零个部分和
T accumulator = (tid < numPartials) ? partialSums[tid] : static_cast<T>(0);
localCache[tid] = accumulator;
__syncthreads();
// 树状归约
for (int offset = blockSize / 2; offset > 0; offset >>= 1) {
if (tid < offset) {
localCache[tid] += localCache[tid + offset];
}
__syncthreads();
}
// 写出最终结果
if (tid == 0) {
*finalResult = localCache[0] + initVal;
}
}
#endif
constexpr double REDUCE_ERROR_TOLERANCE = 0.005;
// ============================================================================
// ReduceSum算法实现接口
// ============================================================================
template <typename InputT = float, typename OutputT = float>
class ReduceSumAlgorithm {
public:
// 成员变量
double* intermediateBuffer_small{nullptr}; // 小数据策略的中间缓冲区
OutputT* intermediateBuffer_large{nullptr}; // 大数据策略的中间缓冲区
int allocatedBlocks_small{0};
int allocatedBlocks_large{0};
// 默认构造函数
ReduceSumAlgorithm() = default;
// 析构函数
~ReduceSumAlgorithm() {
if (intermediateBuffer_small) {
mcFree(intermediateBuffer_small);
}
if (intermediateBuffer_large) {
mcFree(intermediateBuffer_large);
}
}
// 主要接口函数
void reduce(const InputT* d_in, OutputT* d_out, int num_items, OutputT init_value) {
#if !USE_DEFAULT_REF_IMPL
// ========================================
// 参赛者自定义实现区域
// ========================================
if (num_items == 0) {
// 空输入处理
MACA_CHECK(mcMemcpy(d_out, &init_value, sizeof(OutputT), mcMemcpyHostToDevice));
return;
}
// 数据规模分界点
const int THRESHOLD = 1000000;
if (num_items <= THRESHOLD) {
// ========================================
// 策略1: 小数据量优化 (<=1,000,000)
// ========================================
const int THREADS_PER_BLOCK = 256;
const int MAX_BLOCKS = 4096;
int numBlocks = std::min(MAX_BLOCKS, (num_items + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK);
// 共享内存只需要32个doublewarp数量
const int SHARED_MEM_SIZE = 32 * sizeof(double);
// 分配中间缓冲区
if (allocatedBlocks_small < numBlocks) {
if (intermediateBuffer_small) {
mcFree(intermediateBuffer_small);
}
MACA_CHECK(mcMalloc(&intermediateBuffer_small, numBlocks * sizeof(double)));
allocatedBlocks_small = numBlocks;
}
// 第一阶段:多块并行归约
firstStageReduce_Small<<<numBlocks, THREADS_PER_BLOCK, SHARED_MEM_SIZE>>>(
d_in, intermediateBuffer_small, num_items);
// 第二阶段:最终归约
secondStageReduce_Small<<<1, THREADS_PER_BLOCK, SHARED_MEM_SIZE>>>(
intermediateBuffer_small, d_out, numBlocks, init_value);
} else {
// ========================================
// 策略2: 大数据量优化 (>1,000,000)
// ========================================
const int THREADS_PER_BLOCK = 256;
const int ELEMENTS_PER_THREAD = 4; // 每个线程处理4个元素
// 计算需要的块数
int numBlocks = (num_items + THREADS_PER_BLOCK * ELEMENTS_PER_THREAD - 1)
/ (THREADS_PER_BLOCK * ELEMENTS_PER_THREAD);
// 限制最大块数以优化性能
const int MAX_BLOCKS = 2048;
if (numBlocks > MAX_BLOCKS) {
numBlocks = MAX_BLOCKS;
}
// 分配中间缓冲区
if (allocatedBlocks_large < numBlocks) {
if (intermediateBuffer_large) {
mcFree(intermediateBuffer_large);
}
MACA_CHECK(mcMalloc(&intermediateBuffer_large, numBlocks * sizeof(OutputT)));
allocatedBlocks_large = numBlocks;
}
// 第一阶段归约
const size_t sharedMemSize = THREADS_PER_BLOCK * sizeof(OutputT);
firstStageReduce_Large<OutputT><<<numBlocks, THREADS_PER_BLOCK, sharedMemSize>>>(
reinterpret_cast<const OutputT*>(d_in), intermediateBuffer_large, num_items);
// 第二阶段归约
if (numBlocks == 1) {
// 只有一个块直接添加init_value
OutputT hostTemp;
MACA_CHECK(mcMemcpy(&hostTemp, intermediateBuffer_large, sizeof(OutputT), mcMemcpyDeviceToHost));
hostTemp += init_value;
MACA_CHECK(mcMemcpy(d_out, &hostTemp, sizeof(OutputT), mcMemcpyHostToDevice));
} else {
// 多个块,需要第二阶段归约
int finalBlockSize = 256;
// 根据块数优化最终块大小
if (numBlocks <= 32) {
finalBlockSize = 32;
} else if (numBlocks <= 64) {
finalBlockSize = 64;
} else if (numBlocks <= 128) {
finalBlockSize = 128;
}
const size_t finalSharedMem = finalBlockSize * sizeof(OutputT);
secondStageReduce_Large<OutputT><<<1, finalBlockSize, finalSharedMem>>>(
intermediateBuffer_large, d_out, numBlocks, init_value);
}
MACA_CHECK(mcDeviceSynchronize());
}
#else
// ========================================
// 默认基准实现
// ========================================
auto input_ptr = thrust::device_pointer_cast(d_in);
auto output_ptr = thrust::device_pointer_cast(d_out);
*output_ptr = thrust::reduce(
thrust::device,
input_ptr,
input_ptr + num_items,
static_cast<OutputT>(init_value)
);
#endif
}
static const char* getImplementationStatus() {
#if USE_DEFAULT_REF_IMPL
return "DEFAULT_REF_IMPL";
#else
return "CUSTOM_IMPL";
#endif
}
private:
};
// ============================================================================
// 测试和性能评估
// ============================================================================
bool testCorrectness() {
std::cout << "ReduceSum 正确性测试..." << std::endl;
TestDataGenerator generator;
ReduceSumAlgorithm<float, float> algorithm;
bool allPassed = true;
// 测试不同数据规模
for (int i = 0; i < NUM_TEST_SIZES && i < 2; i++) { // 限制测试规模
int size = std::min(TEST_SIZES[i], 10000);
std::cout << " 测试规模: " << size << std::endl;
// 测试普通数据
{
auto data = generator.generateRandomFloats(size, -10.0f, 10.0f);
float init_value = 1.0f;
// CPU参考计算
double cpu_result = cpuReduceSum(data, static_cast<double>(init_value));
// GPU计算
float *d_in;
float *d_out;
MACA_CHECK(mcMalloc(&d_in, size * sizeof(float)));
MACA_CHECK(mcMalloc(&d_out, sizeof(float)));
MACA_CHECK(mcMemcpy(d_in, data.data(), size * sizeof(float), mcMemcpyHostToDevice));
algorithm.reduce(d_in, d_out, size, init_value);
float gpu_result;
MACA_CHECK(mcMemcpy(&gpu_result, d_out, sizeof(float), mcMemcpyDeviceToHost));
// 验证误差
double relative_error = std::abs(gpu_result - cpu_result) / std::abs(cpu_result);
if (relative_error > REDUCE_ERROR_TOLERANCE) {
std::cout << " 失败: 误差过大 " << relative_error << std::endl;
allPassed = false;
} else {
std::cout << " 通过 (误差: " << relative_error << ")" << std::endl;
}
mcFree(d_in);
mcFree(d_out);
}
// 测试特殊值 (NaN, Inf)
if (size > 100) {
std::cout << " 测试特殊值..." << std::endl;
auto data = generator.generateSpecialFloats(size);
float init_value = 0.0f;
double cpu_result = cpuReduceSum(data, static_cast<double>(init_value));
float *d_in;
float *d_out;
MACA_CHECK(mcMalloc(&d_in, size * sizeof(float)));
MACA_CHECK(mcMalloc(&d_out, sizeof(float)));
MACA_CHECK(mcMemcpy(d_in, data.data(), size * sizeof(float), mcMemcpyHostToDevice));
algorithm.reduce(d_in, d_out, size, init_value);
float gpu_result;
MACA_CHECK(mcMemcpy(&gpu_result, d_out, sizeof(float), mcMemcpyDeviceToHost));
// 对于包含特殊值的情况,检查是否正确处理
if (std::isfinite(cpu_result) && std::isfinite(gpu_result)) {
double relative_error = std::abs(gpu_result - cpu_result) / std::abs(cpu_result);
if (relative_error > REDUCE_ERROR_TOLERANCE) {
std::cout << " 失败: 特殊值处理错误" << std::endl;
allPassed = false;
} else {
std::cout << " 通过 (特殊值处理)" << std::endl;
}
} else {
std::cout << " 通过 (特殊值结果)" << std::endl;
}
mcFree(d_in);
mcFree(d_out);
}
}
return allPassed;
}
void benchmarkPerformance() {
PerformanceDisplay::printReduceSumHeader();
TestDataGenerator generator;
PerformanceMeter meter;
ReduceSumAlgorithm<float, float> algorithm;
const int WARMUP_ITERATIONS = 5;
const int BENCHMARK_ITERATIONS = 10;
// 用于YAML报告的数据收集
std::vector<std::map<std::string, std::string>> perf_data;
for (int i = 0; i < NUM_TEST_SIZES; i++) {
int size = TEST_SIZES[i];
// 生成测试数据
auto data = generator.generateRandomFloats(size);
float init_value = 0.0f;
// 分配GPU内存
float *d_in;
float *d_out;
MACA_CHECK(mcMalloc(&d_in, size * sizeof(float)));
MACA_CHECK(mcMalloc(&d_out, sizeof(float)));
MACA_CHECK(mcMemcpy(d_in, data.data(), size * sizeof(float), mcMemcpyHostToDevice));
// Warmup阶段
for (int iter = 0; iter < WARMUP_ITERATIONS; iter++) {
algorithm.reduce(d_in, d_out, size, init_value);
}
// 正式测试阶段
float total_time = 0;
for (int iter = 0; iter < BENCHMARK_ITERATIONS; iter++) {
meter.startTiming();
algorithm.reduce(d_in, d_out, size, init_value);
total_time += meter.stopTiming();
}
float avg_time = total_time / BENCHMARK_ITERATIONS;
// 计算性能指标
auto metrics = PerformanceCalculator::calculateReduceSum(size, avg_time);
// 显示性能数据
PerformanceDisplay::printReduceSumData(size, avg_time, metrics);
// 收集YAML报告数据
auto entry = YAMLPerformanceReporter::createEntry();
entry["data_size"] = std::to_string(size);
entry["time_ms"] = std::to_string(avg_time);
entry["throughput_gps"] = std::to_string(metrics.throughput_gps);
entry["data_type"] = "float";
perf_data.push_back(entry);
mcFree(d_in);
mcFree(d_out);
}
// 生成YAML性能报告
YAMLPerformanceReporter::generateReduceSumYAML(perf_data, "reduce_sum_performance.yaml");
PerformanceDisplay::printSavedMessage("reduce_sum_performance.yaml");
}
// ============================================================================
// 主函数
// ============================================================================
int main(int argc, char* argv[]) {
std::cout << "=== ReduceSum 算法测试 ===" << std::endl;
// 检查参数
std::string mode = "all";
if (argc > 1) {
mode = argv[1];
}
bool correctness_passed = true;
bool performance_completed = true;
try {
if (mode == "correctness" || mode == "all") {
correctness_passed = testCorrectness();
}
if (mode == "performance" || mode == "all") {
if (correctness_passed || mode == "performance") {
benchmarkPerformance();
} else {
std::cout << "跳过性能测试,因为正确性测试未通过" << std::endl;
performance_completed = false;
}
}
std::cout << "\n=== 测试完成 ===" << std::endl;
std::cout << "实现状态: " << ReduceSumAlgorithm<float, float>::getImplementationStatus() << std::endl;
if (mode == "all") {
std::cout << "正确性: " << (correctness_passed ? "通过" : "失败") << std::endl;
std::cout << "性能测试: " << (performance_completed ? "完成" : "跳过") << std::endl;
}
return correctness_passed ? 0 : 1;
} catch (const std::exception& e) {
std::cerr << "测试出错: " << e.what() << std::endl;
return 1;
}
}