// 循环展开
#include "check.h"
#include "timer.h"
#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>

#define SOFTENING 1e-9f
#define BLOCK_SIZE 512  // 或256、或128，但512最佳
#define MULTIPLE 8      // 线程倍增数，一个body由MULTIPLE个线程负责，使得nbodies = BLOCK_SIZE * MULTIPLE

/*
 * Each body contains x, y, and z coordinate positions,
 * as well as velocities in the x, y, and z directions.
 */

typedef struct {
    float x, y, z, vx, vy, vz;
} Body;

/*
 * Do not modify this function. A constraint of this exercise is
 * that it remain a host function.
 */

void randomizeBodies(float* data, int n)
{
    for (int i = 0; i < n; i++) {
        data[i] = 2.0f * (rand() / (float)RAND_MAX) - 1.0f;
    }
}

/*
 * This function calculates the gravitational impact of all bodies in the system
 * on all others, but does not update their positions.
 */

__global__ void bodyForce(Body* p, float dt, int n, int nBlocks)
{
    // index是该线程负责的body序号
    int index = (threadIdx.x + blockIdx.x * blockDim.x) % n;
    Body body = p[index];

    float Fx = 0.0f;
    float Fy = 0.0f;
    float Fz = 0.0f;

    // 使用线程块的共享内存
    __shared__ float3 tile[BLOCK_SIZE];

    // 该线程应该考虑受力的第一个body
    int start = blockIdx.x / nBlocks * BLOCK_SIZE;

    // 将全局内存的数据拷贝到共享内存
    Body temp = p[threadIdx.x + start];
    tile[threadIdx.x] = make_float3(temp.x, temp.y, temp.z);

    __syncthreads();

    // 循环展开
#pragma unroll 64
    for (int j = 0; j < BLOCK_SIZE; j++) {
        float dx = tile[j].x - body.x;
        float dy = tile[j].y - body.y;
        float dz = tile[j].z - body.z;
        float distSqr = dx * dx + dy * dy + dz * dz + SOFTENING;
        float invDist = rsqrtf(distSqr);
        float invDist3 = invDist * invDist * invDist;
        
        Fx += dx * invDist3;
        Fy += dy * invDist3;
        Fz += dz * invDist3;
    }

    __syncthreads();

    // 计算速度
    atomicAdd(&p[index].vx, dt * Fx);
    atomicAdd(&p[index].vy, dt * Fy);
    atomicAdd(&p[index].vz, dt * Fz);
}

__global__ void integratePosition(Body* p, float dt, int n) {
    int index = threadIdx.x + blockIdx.x * blockDim.x;

    // 计算位置
    p[index].x += p[index].vx * dt;
    p[index].y += p[index].vy * dt;
    p[index].z += p[index].vz * dt;
}

int main(const int argc, const char** argv)
{

    /*
     * Do not change the value for `nBodies` here. If you would like to modify it,
     * pass values into the command line.
     */

    int nBodies = 2 << 11;
    int salt = 0;
    if (argc > 1)
        nBodies = 2 << atoi(argv[1]);

    /*
     * This salt is for assessment reasons. Tampering with it will result in automatic failure.
     */

    if (argc > 2)
        salt = atoi(argv[2]);

    const float dt = 0.01f; // time step
    const int nIters = 10; // simulation iterations

    // 原来线程块的数量nBlocks，倍增线程块的数量allBlocks
    int nBlocks = nBodies / BLOCK_SIZE;
    int allBlocks = nBlocks * MULTIPLE;

    // 所需存储空间字节大小
    int bytes = nBodies * sizeof(Body);

    // malloc申请host内存
    float* buf;
    cudaMallocHost(&buf, bytes);

    // cudaMalloc申请device内存
    float* gpu_buf;
    cudaMalloc(&gpu_buf, bytes);

    // 指针转类型
    Body* p = (Body*)buf;
    Body* gpu_p = (Body*)gpu_buf;

    /*
     * As a constraint of this exercise, `randomizeBodies` must remain a host function.
     */

    randomizeBodies(buf, 6 * nBodies); // Init pos / vel data

    // 将host拷贝到device
    cudaMemcpy(gpu_p, p, bytes, cudaMemcpyHostToDevice);

    double totalTime = 0.0;

    /*
     * This simulation will run for 10 cycles of time, calculating gravitational
     * interaction amongst bodies, and adjusting their positions to reflect.
     */

    /*******************************************************************/
    // Do not modify these 2 lines of code.
    for (int iter = 0; iter < nIters; iter++) {
        StartTimer();
        /*******************************************************************/

        /*
         * You will likely wish to refactor the work being done in `bodyForce`,
         * as well as the work to integrate the positions.
         */

        bodyForce<<<allBlocks, BLOCK_SIZE>>>(gpu_p, dt, nBodies, nBlocks);

        /*
         * This position integration cannot occur until this round of `bodyForce` has completed.
         * Also, the next round of `bodyForce` cannot begin until the integration is complete.
         */

        integratePosition<<<nBlocks * 3, BLOCK_SIZE>>>(gpu_p, dt, nBodies);
        

        if (iter == nIters - 1) {
            // 将device拷贝到host
            cudaMemcpy(buf, gpu_buf, bytes, cudaMemcpyDeviceToHost);
        }

        /*******************************************************************/
        // Do not modify the code in this section.
        const double tElapsed = GetTimer() / 1000.0;
        totalTime += tElapsed;
    }

    double avgTime = totalTime / (double)(nIters);
    float billionsOfOpsPerSecond = 1e-9 * nBodies * nBodies / avgTime;

#ifdef ASSESS
    checkPerformance(buf, billionsOfOpsPerSecond, salt);
#else
    checkAccuracy(buf, nBodies);
    printf("%d Bodies: average %0.3f Billion Interactions / second\n", nBodies, billionsOfOpsPerSecond);
    salt += 1;
#endif
    /*******************************************************************/

    /*
     * Feel free to modify code below.
     */

    cudaFreeHost(buf);
    cudaFree(gpu_buf);
}
