Résoudre l'équation de Poisson en utilisant FFT avec CUDA

Question

je suis un tutoriel sur l'utilisation de la bibliothèque cuFFT ici: http://gpgpu.org/static/sc2007/SC07_CUDA_3_Libraries.pdf

Après avoir suivi la ligne par ligne de son code, je suis vraiment obtenir des résultats étranges.

J'ai des données d'entrée qui est un tableau NxN de float s. Le programme effectue une transformée directe FFT, résout l'équation de Poisson, puis fait un inverse FFT. Les données d'entrée (et les données de sortie) sont appelées une image carrée avec une longueur latérale N. Lorsque je commente solve_poisson <<<dimGrid, dimBlock>>> (r_complex_d, kx_d, ky_d, N);, il transfère correctement les données vers l'avant puis effectue une transformation inverse, ce qui entraîne les données de sortie à être les mêmes que les données d'entrée. C'est censé arriver.

Voici la sortie sans appeler la méthode solve_poisson.

0  r_initial: 0.00125126 r: 0.00125132 
1  r_initial: 0.563585  r: 0.563585 
2  r_initial: 0.193304  r: 0.193304 
3  r_initial: 0.80874  r: 0.80874 
4  r_initial: 0.585009  r: 0.585009 
5  r_initial: 0.479873  r: 0.479873 
6  r_initial: 0.350291  r: 0.350291 
7  r_initial: 0.895962  r: 0.895962 
8  r_initial: 0.82284  r: 0.82284 
9  r_initial: 0.746605  r: 0.746605 
10  r_initial: 0.174108  r: 0.174108 
11  r_initial: 0.858943  r: 0.858943 
12  r_initial: 0.710501  r: 0.710502 
13  r_initial: 0.513535  r: 0.513535 
14  r_initial: 0.303995  r: 0.303995 
15  r_initial: 0.0149846 r: 0.0149846 
Press any key to continue . . . 

Cependant, quand je décommentez la méthode solve_poisson, les données de sortie est inf ou nan, ce qui me porte à croire que la variable d'échelle était en quelque sorte proche de zéro dans la méthode solve_poisson. J'ai donc changé float scale = -(kx[idx] * kx[idx] + ky[idy] * ky[idy]); en float scale = -(kx[idx] * kx[idx] + ky[idy] * ky[idy]) + 0.00001f. Cette modification ne figure pas dans le didacticiel d'origine. Les résultats calculés ici ne sont pas supposés avoir des valeurs extrêmes positives ou négatives.

0  r_initial: 0.00125126 r: -11448.1 
1  r_initial: 0.563585  r: 11449.3 
2  r_initial: 0.193304  r: -11448.3 
3  r_initial: 0.80874  r: 11449.2 
4  r_initial: 0.585009  r: 11449.4 
5  r_initial: 0.479873  r: -11448.4 
6  r_initial: 0.350291  r: 11449.5 
7  r_initial: 0.895962  r: -11448.6 
8  r_initial: 0.82284  r: -11448.5 
9  r_initial: 0.746605  r: 11449.4 
10  r_initial: 0.174108  r: -11448.3 
11  r_initial: 0.858943  r: 11449.3 
12  r_initial: 0.710501  r: 11449.2 
13  r_initial: 0.513535  r: -11448.4 
14  r_initial: 0.303995  r: 11449.3 
15  r_initial: 0.0149846 r: -11448.1 
Press any key to continue . . . 

Dans le didacticiel, un exemple de calcul sur la diapositive 43 à la page 22 est computed=0.975879 reference=0.975882, mais mes résultats sont complètement différents et vraiment grand.

Le code suivant est ce que j'ai utilisé.

#include <cuda_runtime.h> 
#include <device_launch_parameters.h> 
#include <cufft.h> 

#include <stdlib.h> 
#include <iostream> 

#define N 4 //4 X 4 // N is the sidelength of the image -> 16 pixels in entire image 
#define block_size_x 2 
#define block_size_y 2 

__global__ void real2complex(cufftComplex *c, float *a, int n); 
__global__ void complex2real_scaled(float *a, cufftComplex *c, float scale, int n); 
__global__ void solve_poisson(cufftComplex *c, float *kx, float *ky, int n); 


int main() 
{ 
    float *kx, *ky, *r; 
    kx = (float *)malloc(sizeof(float) * N); 
    ky = (float *)malloc(sizeof(float) * N); 
    r = (float *)malloc(sizeof(float) * N * N); 

    float *kx_d, *ky_d, *r_d; 
    cufftComplex *r_complex_d; 
    cudaMalloc((void **)&kx_d, sizeof(float) * N); 
    cudaMalloc((void **)&ky_d, sizeof(float) * N); 
    cudaMalloc((void **)&r_d, sizeof(float) * N * N); 
    cudaMalloc((void **)&r_complex_d, sizeof(cufftComplex) * N * N); 

    for (int y = 0; y < N; y++) 
     for (int x = 0; x < N; x++) 
      r[x + y * N] = rand()/(float)RAND_MAX; 
      //r[x + y * N] = sin(exp(-((x - N/2.0f) * (x - N/2.0f) + (N/2.0f - y) * (N/2.0f - y))/(20 * 20))) * 255/sin(1); //Here is sample data that will high values at the center of the image and low values as you go farther and farther away from the center. 

    float* r_inital = (float *)malloc(sizeof(float) * N * N); 
    for (int i = 0; i < N * N; i++) 
     r_inital[i] = r[i]; 

    for (int i = 0; i < N; i++) 
    { 
     kx[i] = i - N/2.0f; //centers kx values to be at center of image 
     ky[i] = N/2.0f - i; //centers ky values to be at center of image 
    } 

    cudaMemcpy(kx_d, kx, sizeof(float) * N, cudaMemcpyHostToDevice); 
    cudaMemcpy(ky_d, ky, sizeof(float) * N, cudaMemcpyHostToDevice); 
    cudaMemcpy(r_d, r, sizeof(float) * N * N, cudaMemcpyHostToDevice); 

    cufftHandle plan; 
    cufftPlan2d(&plan, N, N, CUFFT_C2C); 

    /* Compute the execution configuration, block_size_x*block_size_y = number of threads */ 
    dim3 dimBlock(block_size_x, block_size_y); 
    dim3 dimGrid(N/dimBlock.x, N/dimBlock.y); 
    /* Handle N not multiple of block_size_x or block_size_y */ 
    if (N % block_size_x != 0) dimGrid.x += 1; 
    if (N % block_size_y != 0) dimGrid.y += 1; 

    real2complex << < dimGrid, dimBlock >> > (r_complex_d, r_d, N); 

    cufftExecC2C(plan, r_complex_d, r_complex_d, CUFFT_FORWARD); 
    solve_poisson << <dimGrid, dimBlock >> > (r_complex_d, kx_d, ky_d, N); 
    cufftExecC2C(plan, r_complex_d, r_complex_d, CUFFT_INVERSE); 

    float scale = 1.0f/(N * N); 
    complex2real_scaled << <dimGrid, dimBlock >> > (r_d, r_complex_d, scale, N); 

    cudaMemcpy(r, r_d, sizeof(float) * N * N, cudaMemcpyDeviceToHost); 

    for (int i = 0; i < N * N; i++) 
     std::cout << i << "\tr_initial: " << r_inital[i] << "\tr: " << r[i] << std::endl; 
    system("pause"); 

    /* Destroy plan and clean up memory on device*/ 
    free(kx); 
    free(ky); 
    free(r); 
    free(r_inital); 
    cufftDestroy(plan); 
    cudaFree(r_complex_d); 
    cudaFree(kx_d); 
} 

__global__ void real2complex(cufftComplex *c, float *a, int n) 
{ 
    /* compute idx and idy, the location of the element in the original NxN array */ 
    int idx = blockIdx.x * blockDim.x + threadIdx.x; 
    int idy = blockIdx.y * blockDim.y + threadIdx.y; 
    if (idx < n && idy < n) 
    { 
     int index = idx + idy * n; 
     c[index].x = a[index]; 
     c[index].y = 0.0f; 
    } 
} 

__global__ void complex2real_scaled(float *a, cufftComplex *c, float scale, int n) 
{ 
    /* compute idx and idy, the location of the element in the original NxN array */ 
    int idx = blockIdx.x * blockDim.x + threadIdx.x; 
    int idy = blockIdx.y * blockDim.y + threadIdx.y; 
    if (idx < n && idy < n) 
    { 
     int index = idx + idy * n; 
     a[index] = scale * c[index].x; 
    } 
} 


__global__ void solve_poisson(cufftComplex *c, float *kx, float *ky, int n) 
{ 
    /* compute idx and idy, the location of the element in the original NxN array */ 
    int idx = blockIdx.x * blockDim.x + threadIdx.x; 
    int idy = blockIdx.y * blockDim.y + threadIdx.y; 
    if (idx < n && idy < n) 
    { 
     int index = idx + idy * n; 
     float scale = -(kx[idx] * kx[idx] + ky[idy] * ky[idy]) + 0.00001f; 
     if (idx == 0 && idy == 0) scale = 1.0f; 
     scale = 1.0f/scale; 
     c[index].x *= scale; 
     c[index].y *= scale; 
    } 
} 

Y a-t-il quelque chose que j'ai raté? J'apprécierais vraiment que quelqu'un puisse m'aider.

Answer 1

Comme le montre le tutorial, l'implémentation de Matlab sur la diapositive 33 à la page 17 montre que les calculs de Poisson sont basés sur le coin supérieur gauche de l'écran comme origine. Les valeurs de données x et y sont alors x = (0:(N-1))*h; et y = (0:(N-1))*h;, ce qui explique pourquoi le maillage créé à partir de ces valeurs x et y commence et augmente de 0, comme indiqué sur les axes x et y du graphique sur la diapositive 31 de la page 16. Dans ce cas , où la longueur de l'image était de 1 (je me réfère aux données d'entrée de la matrice flottante NxN ou de la grille de maillage comme image), le centre de l'image est en réalité (0.5, 0.5). Je voulais traduire ces points, donc le point central serait à la place (0,0) et suivrait une représentation typique du Cartesian Plane.

Donc dans mon code, au lieu du code Matlab de

x = (0:(N-1))*h; 
y = (0:(N-1))*h; 

qui pourrait être mis en œuvre comme

for (int i = 0; i < N; i++) 
    { 
     kx[i] = i; 
     ky[i] = i; 
    } 

Je l'ai remplacé avec

for (int i = 0; i < N; i++) 
    { 
     kx[i] = i - N/2.0f; //centers kx values to be at center of image 
     ky[i] = N/2.0f - i; //centers ky values to be at center of image 
    } 

Cependant, je l'avais oublié pour changer le calcul de Poisson afin qu'il reconnaisse le centre de l'image comme l'origine au lieu du r coin coin comme origine. Donc, comme M.Robert Crovella dit, je dois

changez cette ligne: if (idx == 0 && idy == 0) scale = 1.0f; à ceci: if (idx == 2 && idy == 2) scale = 1.0f;

pour le cas où la longueur de l'image ou N, est 4. Généraliser ce pour toute image longueur, cette ligne de code pourrait alors être changée en if (idx == n/2 && idy == n/2) scale = 1.0f;

Answer 2

Bien que l'auteur ait trouvé l'erreur lui-même, je veux partager ma propre implémentation du solveur 2D d'équation de Poisson.

La mise en œuvre diffère légèrement de celle liée à l'affiche. La théorie est rapportée à Solve Poisson Equation Using FFT.

Matlab VERSION

Je premier rapport la version Matlab pour référence:

clear all 
close all 
clc 

M  = 64;    % --- Number of Fourier harmonics along x (should be a multiple of 2) 
N  = 32;    % --- Number of Fourier harmonics along y (should be a multiple of 2) 
Lx  = 3;    % --- Domain size along x 
Ly  = 1.5;    % --- Domain size along y 
sigma = 0.1;    % --- Characteristic width of f (make << 1) 

% --- Wavenumbers 
kx = (2 * pi/Lx) * [0 : (M/2 - 1) (- M/2) : (-1)]; % --- Wavenumbers along x 
ky = (2 * pi/Ly) * [0 : (N/2 - 1) (- N/2) : (-1)]; % --- Wavenumbers along y 
[Kx, Ky] = meshgrid(kx, ky); 

% --- Right-hand side of differential equation 
hx    = Lx/M;     % --- Grid spacing along x 
hy    = Ly/N;     % --- Grid spacing along y 
x    = (0 : (M - 1)) * hx; 
y    = (0 : (N - 1)) * hy; 
[X, Y]   = meshgrid(x, y); 
rSquared  = (X - 0.5 * Lx).^2 + (Y - 0.5 * Ly).^2; 
sigmaSquared = sigma^2; 
f    = exp(-rSquared/(2 * sigmaSquared)) .* (rSquared - 2 * sigmaSquared)/(sigmaSquared^2); 
fHat   = fft2(f); 

% --- Denominator of the unknown spectrum 
den    = -(Kx.^2 + Ky.^2); 
den(1, 1)  = 1;   % --- Avoid division by zero at wavenumber (0, 0) 

% --- Unknown determination 
uHat   = ifft2(fHat ./ den); 
% uHat(1, 1)  = 0;   % --- Force the unknown spectrum at (0, 0) to be zero 
u    = real(uHat); 
u    = u - u(1,1); % --- Force arbitrary constant to be zero by forcing u(1, 1) = 0 

% --- Plots 
uRef = exp(-rSquared/(2 * sigmaSquared)); 
err  = 100 * sqrt(sum(sum(abs(u - uRef).^2))/sum(sum(abs(uRef).^2))); 
errMax = norm(u(:)-uRef(:),inf) 
fprintf('Percentage root mean square error = %f\n', err); 
fprintf('Maximum error = %f\n', errMax); 
surf(X, Y, u) 
xlabel('x') 
ylabel('y') 
zlabel('u') 
title('Solution of 2D Poisson equation by spectral method') 

CUDA VERSION

Voici la version CUDA correspondant:

#include <stdio.h> 
#include <fstream> 
#include <iomanip> 

// --- Greek pi 
#define _USE_MATH_DEFINES 
#include <math.h> 

#include <cufft.h> 

#define BLOCKSIZEX  16 
#define BLOCKSIZEY  16 

#define prec_save 10 

/*******************/ 
/* iDivUp FUNCTION */ 
/*******************/ 
int iDivUp(int a, int b){ return ((a % b) != 0) ? (a/b + 1) : (a/b); } 

/********************/ 
/* CUDA ERROR CHECK */ 
/********************/ 
// --- Credit to http://stackoverflow.com/questions/14038589/what-is-the-canonical-way-to-check-for-errors-using-the-cuda-runtime-api 
void gpuAssert(cudaError_t code, const char *file, int line, bool abort = true) 
{ 
    if (code != cudaSuccess) 
    { 
     fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line); 
     if (abort) { exit(code); } 
    } 
} 

void gpuErrchk(cudaError_t ans) { gpuAssert((ans), __FILE__, __LINE__); } 

/**************************************************/ 
/* COMPUTE RIGHT HAND SIDE OF 2D POISSON EQUATION */ 
/**************************************************/ 
__global__ void computeRHS(const float * __restrict__ d_x, const float * __restrict__ d_y, 
          float2 * __restrict__ d_r, const float Lx, const float Ly, const float sigma, 
          const int M, const int N) { 

    const int tidx = threadIdx.x + blockIdx.x * blockDim.x; 
    const int tidy = threadIdx.y + blockIdx.y * blockDim.y; 

    if ((tidx >= M) || (tidy >= N)) return; 

    const float sigmaSquared = sigma * sigma; 

    const float rSquared = (d_x[tidx] - 0.5f * Lx) * (d_x[tidx] - 0.5f * Lx) + 
          (d_y[tidy] - 0.5f * Ly) * (d_y[tidy] - 0.5f * Ly); 

    d_r[tidy * M + tidx].x = expf(-rSquared/(2.f * sigmaSquared)) * (rSquared - 2.f * sigmaSquared)/(sigmaSquared * sigmaSquared); 
    d_r[tidy * M + tidx].y = 0.f; 

} 

/****************************************************/ 
/* SOLVE 2D POISSON EQUATION IN THE SPECTRAL DOMAIN */ 
/****************************************************/ 
__global__ void solvePoisson(const float * __restrict__ d_kx, const float * __restrict__ d_ky, 
           float2 * __restrict__ d_r, const int M, const int N) 
{ 
    const int tidx = threadIdx.x + blockIdx.x * blockDim.x; 
    const int tidy = threadIdx.y + blockIdx.y * blockDim.y; 

    if ((tidx >= M) || (tidy >= N)) return; 

    float scale = -(d_kx[tidx] * d_kx[tidx] + d_ky[tidy] * d_ky[tidy]); 

    if (tidx == 0 && tidy == 0) scale = 1.f; 

    scale = 1.f/scale; 
    d_r[M * tidy + tidx].x *= scale; 
    d_r[M * tidy + tidx].y *= scale; 

} 

/****************************************************************************/ 
/* SOLVE 2D POISSON EQUATION IN THE SPECTRAL DOMAIN - SHARED MEMORY VERSION */ 
/****************************************************************************/ 
__global__ void solvePoissonShared(const float * __restrict__ d_kx, const float * __restrict__ d_ky, 
    float2 * __restrict__ d_r, const int M, const int N) 
{ 
    const int tidx = threadIdx.x + blockIdx.x * blockDim.x; 
    const int tidy = threadIdx.y + blockIdx.y * blockDim.y; 

    if ((tidx >= M) || (tidy >= N)) return; 

    // --- Use shared memory to minimize multiple access to same spectral coordinate values 
    __shared__ float kx_s[BLOCKSIZEX], ky_s[BLOCKSIZEY]; 

    kx_s[threadIdx.x] = d_kx[tidx]; 
    ky_s[threadIdx.y] = d_ky[tidy]; 
    __syncthreads(); 

    float scale = -(kx_s[threadIdx.x] * kx_s[threadIdx.x] + ky_s[threadIdx.y] * ky_s[threadIdx.y]); 

    if (tidx == 0 && tidy == 0) scale = 1.f; 

    scale = 1.f/scale; 
    d_r[M * tidy + tidx].x *= scale; 
    d_r[M * tidy + tidx].y *= scale; 

} 

/******************************/ 
/* COMPLEX2REAL SCALED KERNEL */ 
/******************************/ 
__global__ void complex2RealScaled(float2 * __restrict__ d_r, float * __restrict__ d_result, const int M, const int N, float scale) 
{ 
    const int tidx = threadIdx.x + blockIdx.x * blockDim.x; 
    const int tidy = threadIdx.y + blockIdx.y * blockDim.y; 

    if ((tidx >= M) || (tidy >= N)) return; 

    d_result[tidy * M + tidx] = scale * (d_r[tidy * M + tidx].x - d_r[0].x); 
} 

/******************************************/ 
/* COMPLEX2REAL SCALED KERNEL - OPTIMIZED */ 
/******************************************/ 
__global__ void complex2RealScaledOptimized(float2 * __restrict__ d_r, float * __restrict__ d_result, const int M, const int N, float scale) 
{ 
    const int tidx = threadIdx.x + blockIdx.x * blockDim.x; 
    const int tidy = threadIdx.y + blockIdx.y * blockDim.y; 

    if ((tidx >= M) || (tidy >= N)) return; 

    __shared__ float d_r_0[1]; 

    if (threadIdx.x == 0) d_r_0[0] = d_r[0].x; 

    volatile float2 c; 
    c.x = d_r[tidy * M + tidx].x; 
    c.y = d_r[tidy * M + tidx].y; 

    d_result[tidy * M + tidx] = scale * (c.x - d_r_0[0]); 
} 

/**************************************/ 
/* SAVE FLOAT2 ARRAY FROM GPU TO FILE */ 
/**************************************/ 
void saveGPUcomplextxt(const float2 * d_in, const char *filename, const int M) { 

    float2 *h_in = (float2 *)malloc(M * sizeof(float2)); 

    gpuErrchk(cudaMemcpy(h_in, d_in, M * sizeof(float2), cudaMemcpyDeviceToHost)); 

    std::ofstream outfile; 
    outfile.open(filename); 
    for (int i = 0; i < M; i++) { 
     //printf("%f %f\n", h_in[i].c.x, h_in[i].c.y); 
     outfile << std::setprecision(prec_save) << h_in[i].x << "\n"; outfile << std::setprecision(prec_save) << h_in[i].y << "\n"; 
    } 
    outfile.close(); 

} 

/*************************************/ 
/* SAVE FLOAT ARRAY FROM GPU TO FILE */ 
/*************************************/ 
template <class T> 
void saveGPUrealtxt(const T * d_in, const char *filename, const int M) { 

    T *h_in = (T *)malloc(M * sizeof(T)); 

    gpuErrchk(cudaMemcpy(h_in, d_in, M * sizeof(T), cudaMemcpyDeviceToHost)); 

    std::ofstream outfile; 
    outfile.open(filename); 
    for (int i = 0; i < M; i++) outfile << std::setprecision(prec_save) << h_in[i] << "\n"; 
    outfile.close(); 

} 

/********/ 
/* MAIN */ 
/********/ 
int main() 
{ 
    const int M  = 64;    // --- Number of Fourier harmonics along x (should be a multiple of 2) 
    const int N  = 32;    // --- Number of Fourier harmonics along y(should be a multiple of 2) 
    const float Lx  = 3.f;    // --- Domain size along x 
    const float Ly  = 1.5f;   // --- Domain size along y 
    const float sigma = 0.1f;   // --- Characteristic width of f(make << 1) 

    // --- Wavenumbers on the host 
    float *h_kx = (float *)malloc(M * sizeof(float)); 
    float *h_ky = (float *)malloc(N * sizeof(float)); 
    for (int k = 0; k < M/2; k++) h_kx[k]  = (2.f * M_PI/Lx) * k; 
    for (int k = -M/2; k < 0; k++) h_kx[k + M] = (2.f * M_PI/Lx) * k; 
    for (int k = 0; k < N/2; k++) h_ky[k]  = (2.f * M_PI/Ly) * k; 
    for (int k = -N/2; k < 0; k++) h_ky[k + N] = (2.f * M_PI/Ly) * k; 

    // --- Wavenumbers on the device 
    float *d_kx; gpuErrchk(cudaMalloc(&d_kx, M * sizeof(float))); 
    float *d_ky; gpuErrchk(cudaMalloc(&d_ky, N * sizeof(float))); 
    gpuErrchk(cudaMemcpy(d_kx, h_kx, M * sizeof(float), cudaMemcpyHostToDevice)); 
    gpuErrchk(cudaMemcpy(d_ky, h_ky, N * sizeof(float), cudaMemcpyHostToDevice)); 

    // --- Domain discretization on the host 
    float *h_x = (float *)malloc(M * sizeof(float)); 
    float *h_y = (float *)malloc(N * sizeof(float)); 
    for (int k = 0; k < M; k++) h_x[k] = Lx/(float)M * k; 
    for (int k = 0; k < N; k++) h_y[k] = Ly/(float)N * k; 

    // --- Domain discretization on the device 
    float *d_x;  gpuErrchk(cudaMalloc(&d_x, M * sizeof(float))); 
    float *d_y;  gpuErrchk(cudaMalloc(&d_y, N * sizeof(float))); 
    gpuErrchk(cudaMemcpy(d_x, h_x, M * sizeof(float), cudaMemcpyHostToDevice)); 
    gpuErrchk(cudaMemcpy(d_y, h_y, N * sizeof(float), cudaMemcpyHostToDevice)); 

    // --- Compute right-hand side of the equation on the device 
    float2 *d_r; gpuErrchk(cudaMalloc(&d_r, M * N * sizeof(float2))); 
    dim3 dimBlock(BLOCKSIZEX, BLOCKSIZEY); 
    dim3 dimGrid(iDivUp(M, BLOCKSIZEX), iDivUp(N, BLOCKSIZEY)); 
    computeRHS << <dimGrid, dimBlock >> >(d_x, d_y, d_r, Lx, Ly, sigma, M, N); 
    gpuErrchk(cudaPeekAtLastError()); 
    gpuErrchk(cudaDeviceSynchronize()); 

    // --- Create plan for CUDA FFT 
    cufftHandle plan; 
    cufftPlan2d(&plan, N, M, CUFFT_C2C); 

    // --- Compute in place forward FFT of right-hand side 
    cufftExecC2C(plan, d_r, d_r, CUFFT_FORWARD); 

    // --- Solve Poisson equation in Fourier space 
    //solvePoisson << <dimGrid, dimBlock >> > (d_kx, d_ky, d_r, M, N); 
    solvePoissonShared << <dimGrid, dimBlock >> > (d_kx, d_ky, d_r, M, N); 

    // --- Compute in place inverse FFT 
    cufftExecC2C(plan, d_r, d_r, CUFFT_INVERSE); 

    //saveGPUcomplextxt(d_r, "D:\\Project\\poisson2DFFT\\poisson2DFFT\\d_r.txt", M * N); 

    // --- With cuFFT, an FFT followed by an IFFT will return the same array times the size of the transform 
    // --- Accordingly, we need to scale the result. 
    const float scale = 1.f/((float)M * (float)N); 
    float *d_result; gpuErrchk(cudaMalloc(&d_result, M * N * sizeof(float))); 
    //complex2RealScaled << <dimGrid, dimBlock >> > (d_r, d_result, M, N, scale); 
    complex2RealScaledOptimized << <dimGrid, dimBlock >> > (d_r, d_result, M, N, scale); 

    //saveGPUrealtxt(d_result, "D:\\Project\\poisson2DFFT\\poisson2DFFT\\d_result.txt", M * N); 

    // --- Transfer data from device to host 
    float *h_result = (float *)malloc(M * N * sizeof(float)); 
    gpuErrchk(cudaMemcpy(h_result, d_result, M * N * sizeof(float), cudaMemcpyDeviceToHost)); 

    return 0; 

}