第九章：高性能计算

// 加法：c = a + b
__m256 c = _mm256_add_ps(a, b);

// 乘加：d = a * b + c
__m256 d = _mm256_fmadd_ps(a, b, c);

// 广播：将标量复制到向量所有位置
__m256 vec = _mm256_set1_ps(scalar);

struct BadLayout {
    int thread1_data;  // 线程1写
    int thread2_data;  // 线程2写，但在同一缓存行！
};

struct GoodLayout {
    alignas(64) int thread1_data;  // 独占缓存行
    alignas(64) int thread2_data;  // 独占缓存行
};

// 绑定线程到NUMA节点
numa_run_on_node(node_id);

// 在指定节点分配内存
void* ptr = numa_alloc_onnode(size, node_id);

// 内存访问模式优化
#pragma omp parallel for schedule(static)
for (int i = 0; i < n; i++) {
    // 确保线程访问本地数据
    process_local_data(i);
}

# 使用perf统计缓存失效
perf stat -e cache-misses,cache-references ./program

# 监控分支预测
perf stat -e branch-misses,branch-instructions ./program

# 内存带宽测量
perf stat -e uncore_imc/data_reads/,uncore_imc/data_writes/ ./program

# 记录性能数据
perf record -g ./simulation

# 生成火焰图
perf script | flamegraph.pl > flame.svg

# 查看热点函数
perf report --stdio

#pragma omp parallel for
for (int i = 0; i < n; i++) {
    process(data[i]);  // 自动分配给多个线程
}

double sum = 0.0;
#pragma omp parallel for reduction(+:sum)
for (int i = 0; i < n; i++) {
    sum += array[i];
}

// 点对点通信
MPI_Send(buffer, count, MPI_DOUBLE, dest, tag, MPI_COMM_WORLD);
MPI_Recv(buffer, count, MPI_DOUBLE, source, tag, MPI_COMM_WORLD, &status);

// 集合通信
MPI_Bcast(buffer, count, MPI_DOUBLE, root, MPI_COMM_WORLD);
MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

MPI_Request request;
MPI_Isend(buffer, count, MPI_DOUBLE, dest, tag, MPI_COMM_WORLD, &request);
// 计算与通信重叠
do_computation();
MPI_Wait(&request, &status);

#pragma omp parallel sections
{
    #pragma omp section
    { compute_forces(); }     // 任务1

    #pragma omp section  
    { update_positions(); }   // 任务2
}

#pragma omp parallel for     // 数据并行
for (int i = 0; i < n; i++) {
    apply_constraints(particles[i]);
}

while (true) {
    Task* task = get_next_task();  // 原子获取
    if (!task) break;
    process_task(task);
}

// Kernel启动配置
dim3 blockDim(256);  // 每个block 256线程
dim3 gridDim((n + 255) / 256);  // 足够的blocks
kernel<<<gridDim, blockDim>>>(data, n);

__syncthreads();  // Block内同步
__syncwarp();     // Warp内同步（隐式）
// Grid级同步需要kernel边界

// Good: 合并访问
data[threadIdx.x + blockIdx.x * blockDim.x] = value;

// Bad: 跨步访问
data[threadIdx.x * stride] = value;

// Bad: 随机访问  
data[indices[threadIdx.x]] = value;

// 使用launch bounds限制寄存器
__global__ __launch_bounds__(256, 4)
void kernel() {
    // 编译器将优化到指定的占用率
}

def calc_occupancy(regs_per_thread, shared_per_block, threads_per_block):
    max_warps_per_sm = 64  # 架构相关
    max_regs_per_sm = 65536
    max_shared_per_sm = 49152

    warps_limited_by_regs = max_regs_per_sm // (regs_per_thread * 32)
    blocks_limited_by_shared = max_shared_per_sm // shared_per_block
    warps_limited_by_shared = blocks_limited_by_shared * (threads_per_block // 32)

    active_warps = min(warps_limited_by_regs, warps_limited_by_shared)
    return active_warps / max_warps_per_sm

template<int dim>
class SPGrid {
    static constexpr int page_size = 4096;
    static constexpr int block_size = (dim == 2) ? 4 : 8;

    uint64_t* page_table;  // 页表跟踪
    T* data;               // 虚拟内存池

    size_t offset(const Vec<dim>& idx) {
        // Morton编码保持空间局部性
        return morton_encode(idx);
    }
};

Root (dynamic hash map)
  └── Internal Node (32³ or 16³)
       └── Internal Node (16³ or 8³)  
            └── Leaf Node (8³)

struct HierarchicalGrid {
    struct Node {
        union {
            Node* children[8];  // 内部节点
            T data[8*8*8];      // 叶子节点
        };
        uint8_t mask;           // 子节点存在性

        bool is_leaf() const { return mask & 0x80; }
    };

    void refine(Node* node, const Vec3& pos) {
        if (need_refinement(node, pos)) {
            split_node(node);
        }
    }
};

struct SparseMask {
    uint64_t mask[BLOCKS/64];  // 每bit表示一个block

    bool is_active(int idx) {
        return mask[idx/64] & (1ULL << (idx%64));
    }

    void activate(int idx) {
        mask[idx/64] |= (1ULL << (idx%64));
    }

    int popcount() {  // 活跃block数
        int count = 0;
        for (auto m : mask) {
            count += __builtin_popcountll(m);
        }
        return count;
    }
};

class ConcurrentSparseGrid {
    struct Block {
        std::atomic<uint32_t> version;
        T data[BLOCK_SIZE];
    };

    void update(const Vec3& pos, T value) {
        Block* block = get_or_create_block(pos);
        uint32_t ver = block->version.load();

        // 乐观并发控制
        block->data[local_idx(pos)] = value;

        if (!block->version.compare_exchange_weak(ver, ver+1)) {
            // 重试或使用其他策略
        }
    }
};

// 原始循环
for (int i = 0; i < n; i++) {
    sum += a[i] * b[i];
}

// 4倍展开
for (int i = 0; i < n-3; i += 4) {
    sum += a[i] * b[i] + a[i+1] * b[i+1] + 
           a[i+2] * b[i+2] + a[i+3] * b[i+3];
}
// 处理剩余元素
for (int i = n - n%4; i < n; i++) {
    sum += a[i] * b[i];
}

// 矩阵乘法的缓存优化版本
const int TILE = 64;  // 选择适合L1缓存的块大小

for (int ii = 0; ii < n; ii += TILE) {
    for (int jj = 0; jj < n; jj += TILE) {
        for (int kk = 0; kk < n; kk += TILE) {
            // 处理一个TILE×TILE的块
            for (int i = ii; i < min(ii+TILE, n); i++) {
                for (int j = jj; j < min(jj+TILE, n); j++) {
                    double sum = C[i][j];
                    for (int k = kk; k < min(kk+TILE, n); k++) {
                        sum += A[i][k] * B[k][j];
                    }
                    C[i][j] = sum;
                }
            }
        }
    }
}

// 分离的循环
for (int i = 0; i < n; i++) {
    a[i] = b[i] + c[i];
}
for (int i = 0; i < n; i++) {
    d[i] = a[i] * e[i];
}

// 融合后：减少内存流量
for (int i = 0; i < n; i++) {
    a[i] = b[i] + c[i];
    d[i] = a[i] * e[i];  // a[i]还在缓存中
}

// 原始：列优先访问（缓存不友好）
for (int j = 0; j < n; j++) {
    for (int i = 0; i < n; i++) {
        process(A[i][j]);
    }
}

// 交换后：行优先访问（缓存友好）
for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
        process(A[i][j]);
    }
}

// Array of Structures (AoS) - 缓存不友好
struct Particle_AoS {
    float x, y, z;
    float vx, vy, vz;
};
Particle_AoS particles[N];

// Structure of Arrays (SoA) - SIMD友好
struct Particles_SoA {
    float x[N], y[N], z[N];
    float vx[N], vy[N], vz[N];
};

void process_with_prefetch(float* data, int n) {
    const int PREFETCH_DISTANCE = 8;

    for (int i = 0; i < n; i++) {
        // 预取未来的数据
        if (i + PREFETCH_DISTANCE < n) {
            __builtin_prefetch(&data[i + PREFETCH_DISTANCE], 0, 3);
        }

        // 处理当前数据
        expensive_computation(data[i]);
    }
}

// 压缩存储减少内存带宽
struct CompressedParticle {
    // 位置使用相对坐标（16位整数）
    int16_t dx, dy, dz;
    // 速度使用半精度浮点
    half vx, vy, vz;
};

// 解压时转换回全精度
void decompress(CompressedParticle& cp, Particle& p, Vec3 base) {
    p.x = base.x + cp.dx * SCALE;
    p.y = base.y + cp.dy * SCALE;
    p.z = base.z + cp.dz * SCALE;
    p.vx = float(cp.vx);
    p.vy = float(cp.vy);
    p.vz = float(cp.vz);
}

// 标量版本
void add_scalar(float* a, float* b, float* c, int n) {
    for (int i = 0; i < n; i++) {
        c[i] = a[i] + b[i];
    }
}

// AVX向量化版本
void add_avx(float* a, float* b, float* c, int n) {
    int i = 0;
    // 主循环：8个元素一组
    for (; i < n - 7; i += 8) {
        __m256 va = _mm256_load_ps(&a[i]);
        __m256 vb = _mm256_load_ps(&b[i]);
        __m256 vc = _mm256_add_ps(va, vb);
        _mm256_store_ps(&c[i], vc);
    }
    // 处理剩余元素
    for (; i < n; i++) {
        c[i] = a[i] + b[i];
    }
}

// 确保32字节对齐（AVX需要）
float* aligned_alloc_float(size_t n) {
    void* ptr = nullptr;
    if (posix_memalign(&ptr, 32, n * sizeof(float)) != 0) {
        return nullptr;
    }
    return static_cast<float*>(ptr);
}

// 使用aligned属性
struct alignas(32) AlignedData {
    float values[8];
};

// Bad: gather操作性能差
__m256 gather_bad(float* data, int* indices) {
    return _mm256_i32gather_ps(data, _mm256_load_si256((__m256i*)indices), 4);
}

// Good: 重组数据避免gather
void transpose_block_8x8(float* src, float* dst, int stride) {
    __m256 row0 = _mm256_load_ps(src + 0*stride);
    __m256 row1 = _mm256_load_ps(src + 1*stride);
    // ... 加载所有行

    // 使用shuffle和blend进行转置
    // ... 转置操作

    _mm256_store_ps(dst + 0*8, col0);
    _mm256_store_ps(dst + 1*8, col1);
    // ... 存储所有列
}

// MPI非阻塞通信示例
void overlap_compute_comm(float* local_data, float* halo_data, int n) {
    MPI_Request requests[4];

    // 启动非阻塞接收
    MPI_Irecv(halo_left, size, MPI_FLOAT, left_rank, 0, 
              MPI_COMM_WORLD, &requests[0]);
    MPI_Irecv(halo_right, size, MPI_FLOAT, right_rank, 1, 
              MPI_COMM_WORLD, &requests[1]);

    // 启动非阻塞发送
    MPI_Isend(send_left, size, MPI_FLOAT, left_rank, 1, 
              MPI_COMM_WORLD, &requests[2]);
    MPI_Isend(send_right, size, MPI_FLOAT, right_rank, 0, 
              MPI_COMM_WORLD, &requests[3]);

    // 计算内部区域（不依赖halo）
    compute_interior(local_data, n);

    // 等待通信完成
    MPI_Waitall(4, requests, MPI_STATUSES_IGNORE);

    // 计算边界区域（使用halo数据）
    compute_boundary(local_data, halo_left, halo_right, n);
}

template<typename T>
class DoubleBuffer {
    T* buffers[2];
    int current = 0;

public:
    void swap() { current = 1 - current; }
    T* read_buffer() { return buffers[current]; }
    T* write_buffer() { return buffers[1 - current]; }
};

void pipeline_processing() {
    DoubleBuffer<Data> buffer;

    // 启动第一次读取
    async_read(buffer.write_buffer());

    for (int i = 0; i < iterations; i++) {
        buffer.swap();

        // 并行：处理当前数据 + 读取下一批
        auto future = std::async(std::launch::async, [&]() {
            async_read(buffer.write_buffer());
        });

        process(buffer.read_buffer());

        future.wait();
    }
}

void multi_stream_processing(float* h_data, float* d_data, int n, int chunks) {
    const int chunk_size = n / chunks;
    cudaStream_t streams[chunks];

    // 创建流
    for (int i = 0; i < chunks; i++) {
        cudaStreamCreate(&streams[i]);
    }

    // 异步处理
    for (int i = 0; i < chunks; i++) {
        int offset = i * chunk_size;

        // 异步拷贝到GPU
        cudaMemcpyAsync(d_data + offset, h_data + offset, 
                       chunk_size * sizeof(float), 
                       cudaMemcpyHostToDevice, streams[i]);

        // 在流中启动kernel
        process_kernel<<<blocks, threads, 0, streams[i]>>>
                      (d_data + offset, chunk_size);

        // 异步拷贝回CPU
        cudaMemcpyAsync(h_data + offset, d_data + offset, 
                       chunk_size * sizeof(float), 
                       cudaMemcpyDeviceToHost, streams[i]);
    }

    // 同步所有流
    for (int i = 0; i < chunks; i++) {
        cudaStreamSynchronize(streams[i]);
        cudaStreamDestroy(streams[i]);
    }
}

struct Domain1D {
    int global_size;
    int num_gpus;
    int gpu_id;

    int local_size() const {
        return (global_size + num_gpus - 1) / num_gpus;
    }

    int local_start() const {
        return gpu_id * local_size();
    }

    int local_end() const {
        return min(local_start() + local_size(), global_size);
    }
};

struct Domain3D {
    int3 global_size;
    int3 gpu_grid;    // GPU的3D排列
    int3 gpu_coords;  // 当前GPU的坐标

    int3 local_size() const {
        return make_int3(
            (global_size.x + gpu_grid.x - 1) / gpu_grid.x,
            (global_size.y + gpu_grid.y - 1) / gpu_grid.y,
            (global_size.z + gpu_grid.z - 1) / gpu_grid.z
        );
    }

    // 计算表面积/体积比
    float surface_volume_ratio() const {
        int3 ls = local_size();
        int surface = 2 * (ls.x*ls.y + ls.y*ls.z + ls.z*ls.x);
        int volume = ls.x * ls.y * ls.z;
        return float(surface) / volume;
    }
};

class HaloExchange {
    int3 local_size;
    int3 halo_width;
    int3 total_size;  // local + 2*halo

    float* data;

    // 6个方向的邻居GPU
    int neighbors[6];  // -x, +x, -y, +y, -z, +z

    void exchange_halos() {
        // X方向交换
        exchange_x_direction();

        // Y方向交换（包含X方向的halo）
        exchange_y_direction();

        // Z方向交换（包含X,Y方向的halo）
        exchange_z_direction();
    }

    void exchange_x_direction() {
        size_t yz_slice_size = halo_width.x * local_size.y * local_size.z;

        // 准备发送缓冲区
        pack_x_minus(send_buffer_xm);
        pack_x_plus(send_buffer_xp);

        // 非阻塞P2P传输
        cudaMemcpyPeerAsync(recv_xm, gpu_xm, send_xp, gpu_id, 
                           yz_slice_size * sizeof(float), stream_xm);
        cudaMemcpyPeerAsync(recv_xp, gpu_xp, send_xm, gpu_id, 
                           yz_slice_size * sizeof(float), stream_xp);

        // 同步
        cudaStreamSynchronize(stream_xm);
        cudaStreamSynchronize(stream_xp);

        // 解包到halo区域
        unpack_x_minus(recv_buffer_xm);
        unpack_x_plus(recv_buffer_xp);
    }
};

void enable_peer_access(int num_gpus) {
    for (int i = 0; i < num_gpus; i++) {
        cudaSetDevice(i);
        for (int j = 0; j < num_gpus; j++) {
            if (i != j) {
                int can_access;
                cudaDeviceCanAccessPeer(&can_access, i, j);
                if (can_access) {
                    cudaDeviceEnablePeerAccess(j, 0);
                }
            }
        }
    }
}

void timestep_with_overlap() {
    // 1. 启动halo交换（非阻塞）
    start_halo_exchange();

    // 2. 计算内部区域（不需要halo）
    dim3 interior_blocks(
        (local_size.x - 2*halo_width.x + 255) / 256,
        local_size.y - 2*halo_width.y,
        local_size.z - 2*halo_width.z
    );
    compute_interior<<<interior_blocks, 256, 0, compute_stream>>>(
        data, local_size, halo_width
    );

    // 3. 等待halo交换完成
    finish_halo_exchange();

    // 4. 计算边界区域
    compute_boundary_x<<<...>>>();
    compute_boundary_y<<<...>>>();
    compute_boundary_z<<<...>>>();

    // 5. 同步所有计算
    cudaDeviceSynchronize();
}

class MultiGPUManager {
    struct GPUContext {
        int device_id;
        cudaStream_t compute_stream;
        cudaStream_t transfer_stream[2];  // 双缓冲传输

        float* d_data;      // 设备内存
        float* h_buffer[2]; // 钉扎内存缓冲区
    };

    std::vector<GPUContext> contexts;

    void process_multi_gpu() {
        // 为每个GPU创建上下文
        for (int i = 0; i < num_gpus; i++) {
            cudaSetDevice(i);
            contexts[i].compute_stream = create_stream();
            contexts[i].transfer_stream[0] = create_stream();
            contexts[i].transfer_stream[1] = create_stream();

            // 分配钉扎内存加速传输
            cudaHostAlloc(&contexts[i].h_buffer[0], buffer_size, 
                         cudaHostAllocPortable);
            cudaHostAlloc(&contexts[i].h_buffer[1], buffer_size, 
                         cudaHostAllocPortable);
        }

        // Pipeline处理
        for (int iter = 0; iter < num_iterations; iter++) {
            for (int gpu = 0; gpu < num_gpus; gpu++) {
                cudaSetDevice(gpu);
                auto& ctx = contexts[gpu];

                int buf_idx = iter % 2;

                // 异步传输下一批数据
                if (iter + 1 < num_iterations) {
                    prepare_next_data(ctx.h_buffer[1-buf_idx], iter+1);
                    cudaMemcpyAsync(ctx.d_data, ctx.h_buffer[1-buf_idx],
                                   data_size, cudaMemcpyHostToDevice,
                                   ctx.transfer_stream[1-buf_idx]);
                }

                // 在当前数据上计算
                compute_kernel<<<grid, block, 0, ctx.compute_stream>>>
                              (ctx.d_data, params);

                // 等待上一次传输完成
                cudaStreamSynchronize(ctx.transfer_stream[buf_idx]);
            }
        }
    }
};

class EventSync {
    cudaEvent_t events[MAX_GPUS][MAX_EVENTS];

    void cross_gpu_dependency() {
        // GPU 0完成计算
        cudaSetDevice(0);
        kernel_gpu0<<<...>>>();
        cudaEventRecord(events[0][0]);

        // GPU 1等待GPU 0
        cudaSetDevice(1);
        cudaStreamWaitEvent(stream_gpu1, events[0][0]);
        kernel_gpu1<<<..., stream_gpu1>>>();

        // 可以形成复杂的依赖DAG
    }
};

class WorkStealingScheduler {
    struct WorkQueue {
        std::atomic<int> head;
        std::atomic<int> tail;
        Task* tasks[QUEUE_SIZE];

        bool try_push(Task* task) {
            int h = head.load(std::memory_order_acquire);
            int t = tail.load(std::memory_order_relaxed);

            if (t - h >= QUEUE_SIZE) return false;

            tasks[t % QUEUE_SIZE] = task;
            tail.store(t + 1, std::memory_order_release);
            return true;
        }

        Task* try_pop() {
            int h = head.load(std::memory_order_relaxed);
            int t = tail.load(std::memory_order_acquire);

            if (h >= t) return nullptr;

            Task* task = tasks[h % QUEUE_SIZE];
            head.store(h + 1, std::memory_order_release);
            return task;
        }

        Task* try_steal() {
            int h = head.load(std::memory_order_acquire);
            int t = tail.load(std::memory_order_acquire);

            if (h >= t) return nullptr;

            Task* task = tasks[h % QUEUE_SIZE];
            if (head.compare_exchange_strong(h, h+1)) {
                return task;
            }
            return nullptr;
        }
    };

    WorkQueue queues[MAX_GPUS];

    void gpu_worker(int gpu_id) {
        cudaSetDevice(gpu_id);
        WorkQueue& my_queue = queues[gpu_id];

        while (!done) {
            Task* task = my_queue.try_pop();

            if (!task) {
                // 尝试从其他GPU窃取
                for (int victim = 0; victim < num_gpus; victim++) {
                    if (victim != gpu_id) {
                        task = queues[victim].try_steal();
                        if (task) break;
                    }
                }
            }

            if (task) {
                execute_task(task, gpu_id);
            } else {
                std::this_thread::yield();
            }
        }
    }
};

class LoadMonitor {
    struct GPULoad {
        std::atomic<float> utilization;
        std::atomic<int> task_count;
        std::atomic<int64_t> total_time;
    };

    GPULoad loads[MAX_GPUS];

    void update_load(int gpu_id, float util, int64_t task_time) {
        loads[gpu_id].utilization.store(util);
        loads[gpu_id].task_count.fetch_add(1);
        loads[gpu_id].total_time.fetch_add(task_time);
    }

    int find_least_loaded_gpu() {
        float min_load = 1.0f;
        int best_gpu = 0;

        for (int i = 0; i < num_gpus; i++) {
            float load = loads[i].utilization.load();
            if (load < min_load) {
                min_load = load;
                best_gpu = i;
            }
        }

        return best_gpu;
    }

    void rebalance_work() {
        // 计算平均负载
        float avg_load = 0.0f;
        for (int i = 0; i < num_gpus; i++) {
            avg_load += loads[i].utilization.load();
        }
        avg_load /= num_gpus;

        // 识别过载和空闲GPU
        std::vector<int> overloaded, underloaded;
        for (int i = 0; i < num_gpus; i++) {
            float load = loads[i].utilization.load();
            if (load > avg_load * 1.2f) {
                overloaded.push_back(i);
            } else if (load < avg_load * 0.8f) {
                underloaded.push_back(i);
            }
        }

        // 迁移任务
        for (int src : overloaded) {
            for (int dst : underloaded) {
                migrate_tasks(src, dst);
            }
        }
    }
};

void nccl_allreduce_example() {
    ncclComm_t comms[num_gpus];

    // 初始化NCCL
    ncclCommInitAll(comms, num_gpus, device_ids);

    // 每个GPU执行allreduce
    #pragma omp parallel num_threads(num_gpus)
    {
        int gpu = omp_get_thread_num();
        cudaSetDevice(gpu);

        // 局部计算
        compute_local<<<...>>>(d_data[gpu]);

        // 全局归约
        ncclAllReduce(d_data[gpu], d_data[gpu], data_size,
                     ncclFloat, ncclSum, comms[gpu], 
                     streams[gpu]);

        // 等待完成
        cudaStreamSynchronize(streams[gpu]);
    }

    // 清理
    for (int i = 0; i < num_gpus; i++) {
        ncclCommDestroy(comms[i]);
    }
}

struct Particle {
    float x, y, z;
    float vx, vy, vz;
    float mass;
    int type;
};

void update_positions(Particle* p, int n, float dt) {
    for (int i = 0; i < n; i++) {
        p[i].x += p[i].vx * dt;
        p[i].y += p[i].vy * dt;
        p[i].z += p[i].vz * dt;
    }
}

__global__ void transpose(float* out, float* in, int n) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;

    if (x < n && y < n) {
        out[y * n + x] = in[x * n + y];
    }
}