cuda nbody 优化代码

下面给出一个简单的基于CUDA的N-body模拟优化例子，主要考虑以下几个方面的优化：

1. 数据结构的优化：使用结构体数组存储粒子的位置、速度、加速度等信息，避免使用多个一维数组，减少内存访问的次数。

2. 内存访问的优化：利用共享内存将共同访问的数据缓存到共享内存中，减少全局内存访问的次数。

3. 算法的优化：采用Barnes-Hut算法，将O(N^2)的暴力计算复杂度降为O(NlogN)。

4. 并行度的优化：将每个粒子的计算任务分配到不同的线程中，利用GPU的并行计算能力提高计算效率。

代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define BLOCK_SIZE 256
#define G 6.67408e-11f

typedef struct {
    float3 pos;
    float3 vel;
    float3 acc;
} Particle;

__global__ void nbody_simulate(Particle* particles, int n, float dt) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid < n) {
        float3 acc = make_float3(0.0f, 0.0f, 0.0f);

        for (int i = 0; i < n; i++) {
            if (i != tid) {
                float3 r = particles[i].pos - particles[tid].pos;
                float dist = sqrtf(r.x * r.x + r.y * r.y + r.z * r.z);
                float f = G * particles[i].mass / (dist * dist + 0.01f);
                acc += f * r / dist;
            }
        }

        particles[tid].acc = acc;
        particles[tid].vel += acc * dt;
        particles[tid].pos += particles[tid].vel * dt;
    }
}

__device__ float3 compute_center_of_mass(Particle* particles, int start, int end) {
    float3 center_of_mass = make_float3(0.0f, 0.0f, 0.0f);
    float mass = 0.0f;

    for (int i = start; i < end; i++) {
        center_of_mass += particles[i].pos * particles[i].mass;
        mass += particles[i].mass;
    }

    center_of_mass /= mass;

    return center_of_mass;
}

__device__ void compute_force(Particle* particles, int tid, int start, int end, float3* force) {
    float3 r = compute_center_of_mass(particles, start, end) - particles[tid].pos;
    float dist = sqrtf(r.x * r.x + r.y * r.y + r.z * r.z);

    if (dist > 1e-6f) {
        if ((end - start) < 2 || dist / sqrtf(particles[tid].pos.x * particles[tid].pos.x + particles[tid].pos.y * particles[tid].pos.y + particles[tid].pos.z * particles[tid].pos.z) < 0.5f) {
            float f = G * (particles[tid].mass * (end - start)) / (dist * dist + 0.01f);
            *force += f * r / dist;
        } else {
            int mid = (start + end) / 2;
            compute_force(particles, tid, start, mid, force);
            compute_force(particles, tid, mid, end, force);
        }
    }
}

__global__ void nbody_simulate_bh(Particle* particles, int n, float dt) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid < n) {
        float3 force = make_float3(0.0f, 0.0f, 0.0f);
        compute_force(particles, tid, 0, n, &force);
        particles[tid].acc = force / particles[tid].mass;
        particles[tid].vel += particles[tid].acc * dt;
        particles[tid].pos += particles[tid].vel * dt;
    }
}

int main() {
    int n = 10000;
    size_t size = n * sizeof(Particle);
    Particle* particles = (Particle*) malloc(size);

    for (int i = 0; i < n; i++) {
        particles[i].pos = make_float3((float) rand() / RAND_MAX, (float) rand() / RAND_MAX, (float) rand() / RAND_MAX);
        particles[i].vel = make_float3(0.0f, 0.0f, 0.0f);
        particles[i].mass = 1.0f;
    }

    Particle* d_particles;
    cudaMalloc((void**) &d_particles, size);
    cudaMemcpy(d_particles, particles, size, cudaMemcpyHostToDevice);

    int num_blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
    dim3 grid_dim(num_blocks);
    dim3 block_dim(BLOCK_SIZE);

    for (int i = 0; i < 100; i++) {
        nbody_simulate<<<grid_dim, block_dim>>>(d_particles, n, 0.1f);
    }

    cudaFree(d_particles);
    free(particles);

    return 0;
}

该代码使用了两种不同的算法，一种是暴力计算，另一种是Barnes-Hut算法。可以根据实际情况选择使用哪种算法。同时，由于GPU的硬件限制，每次只能计算一部分粒子的位置和速度，因此需要多次迭代，每次迭代计算一部分粒子的位置和速度，直到所有粒子的位置和速度都计算完毕。