ajout correction equation de la chaleur

exercices/03-equation-chaleur/correction/Makefile

+NVCC := nvcc
+CFLAGS := -O3
+SRC := heat_equation.cu
+EXEC := heat_equation
+
+all: $(EXEC)
+
+$(EXEC): $(SRC)
+	$(NVCC) $(CFLAGS) -o $(EXEC) $(SRC)
+
+clean:
+	rm -f $(EXEC)
+
+.PHONY: all clean
\ No newline at end of file

exercices/03-equation-chaleur/correction/heat_equation.cu

+#include<iostream>
+#include<assert.h>
+#include <chrono>
+
+using namespace std;
+
+typedef struct _matrix{
+    int m,n;
+    double** data;
+} matrix_t;
+
+matrix_t matrix_create(int m, int n){
+    matrix_t mat;
+    mat.m = m;
+    mat.n = n;
+    mat.data = (double**)calloc(m, sizeof(double*));
+    mat.data[0] = (double*)calloc(m*n, sizeof(double));
+    for(int i=0; i<mat.m; i++){
+        mat.data[i] = &mat.data[0][i*n];
+    }
+    return mat;
+}
+
+void matrix_destroy(matrix_t mat){
+    free(mat.data[0]);
+    free(mat.data);
+}
+
+void matrix_swap(matrix_t *m1, matrix_t *m2){
+    int temp_m = m1->m; 
+    int temp_n = m1->n;
+    double** temp_data = m1->data;
+    
+    m1->m = m2->m;
+    m1->n = m2->n;
+    m1->data = m2->data;
+    
+    m2->m = temp_m;
+    m2->n = temp_n;
+    m2->data = temp_data;
+}
+
+void set_boundary_conditions(matrix_t U){
+    for(int i=0; i<U.m; i++){
+        U.data[i][0] = 1.0;
+        U.data[i][U.n-1] = 1.0;
+    }
+    for(int j=0; j<U.n; j++){
+        U.data[0][j] = 1.0;
+        U.data[U.m-1][j] = 1.0;
+    }
+}
+
+void cpu_step(matrix_t U, matrix_t Unext){
+    for(int i=0; i<U.m; i++){
+        for(int j=0; j<U.n; j++){
+            if(i==0 || j==0 || i==U.m-1 || j==U.n-1) Unext.data[i][j] = U.data[i][j];
+            else Unext.data[i][j] = (U.data[i+1][j]+U.data[i-1][j]+U.data[i][j+1]+U.data[i][j-1])/4.0;
+        }
+    }
+}
+
+void matrix_print(matrix_t U){
+    for(int i=0; i<U.m; i++){
+        for(int j=0; j<U.n; j++){
+            std::cout << U.data[i][j] << "  ";
+        }
+        std::cout << std::endl;
+    }
+}
+
+
+int double_compare(double* a1, double* a2, int s, double e){
+    int dif = 0; 
+    for(int i=0; i<s; i++){
+        if(fabs(a1[i]-a2[i])>e) dif++;
+    }
+    return dif;
+}
+
+__global__ void gpu_step(double *U, double *Unext, int m, int n){
+    int i = (blockIdx.x*blockDim.x) + threadIdx.x;
+    int j = (blockIdx.y*blockDim.y) + threadIdx.y;
+    if(i>0 && i<m-1 && j>0 && j<n-1){
+        Unext[i*n+j] = (U[(i+1)*n+j]+U[(i-1)*n+j]+U[i*n+j+1]+U[i*n+j-1])/4.0;
+    }
+    if(i==0 || i==m-1 || j==0 || j==n-1){
+        Unext[i*n+j] = U[i*n+j];
+    }
+}
+
+int main(){
+    int m = 1000;
+    int n = 1000;
+    int Ttot = 10000;
+
+    cout << "m : " << m << ", n : " << n << ", Ttot : " << Ttot << endl;
+
+    // alloc and init on CPU
+    matrix_t U = matrix_create(m, n);
+    matrix_t Unext = matrix_create(m, n);
+    set_boundary_conditions(U);
+
+    // alloc and init on GPU
+    double *d_U, *d_Unext;
+    cudaMalloc(&d_U, m*n*sizeof(double));
+    cudaMalloc(&d_Unext, m*n*sizeof(double));
+    cudaMemcpy(d_U, U.data[0], m*n*sizeof(double), cudaMemcpyHostToDevice);
+
+    auto start_cpu = chrono::high_resolution_clock::now();
+    // perform several iterations on CPU
+    cout << "start CPU work" << endl;
+    for(int t=0; t<Ttot; t++){
+         cpu_step(U, Unext);
+         matrix_swap(&U, &Unext);
+    }
+    auto end_cpu = chrono::high_resolution_clock::now();
+
+    
+    auto start_gpu = chrono::high_resolution_clock::now();
+
+    // GPU exec config
+    dim3 dimBlock(10, 10, 1);
+    dim3 dimGrid((m + dimBlock.x - 1) / dimBlock.x, (n + dimBlock.y - 1) / dimBlock.y, 1);  
+
+    // perform several iterations on GPU
+    cout << "start GPU work" << endl;
+    for(int t=0; t<Ttot; t++){
+        gpu_step<<<dimGrid, dimBlock>>>(d_U, d_Unext, m, n); 
+        swap(d_U, d_Unext);
+    }
+
+    // use Unext memory space to gather GPU res and compare results
+    double* gpu_res = Unext.data[0];
+    cudaMemcpy(gpu_res, d_U, m*n*sizeof(double), cudaMemcpyDeviceToHost);
+    auto end_gpu = chrono::high_resolution_clock::now();
+
+    // matrix_print(Unext);
+
+    // ensure CPU and GPU res are close
+    assert(double_compare(gpu_res, U.data[0], m*n, 0.0000001) == 0);
+
+    cout << "OK" << endl;
+    auto duration_cpu = chrono::duration_cast<chrono::milliseconds>(end_cpu - start_cpu);
+    auto duration_gpu = chrono::duration_cast<chrono::milliseconds>(end_gpu - start_gpu);
+    cout << "time on CPU : " << duration_cpu.count() << "ms" << endl;
+    cout << "time on GPU : " << duration_gpu.count() << "ms" << endl;
+
+    matrix_destroy(U);
+    matrix_destroy(Unext);
+    cudaFree(d_U);
+    cudaFree(d_Unext);
+
+}
\ No newline at end of file