Уравнение Пуассона на cuda(проблема с вызовом функции)

Всем дорого времени суток!

Решаю сейчас в рамках курсовой задачу решения уравнения Пуассона, так вот код написал, по идее должен работать но выбивает драйвера на функции Poisson2DSOR. Может кто подскажет, может там с памятью проблемы, или много раз вызываю данную функцию подряд... Да и память я асинхронно выделяю. И кстати функции Init_f и NormaF тоже не считают. Не могу понять что делаю не верно?

Nvidia 310M

  1. #include "cuda_runtime.h"
  2. #include "device_launch_parameters.h"
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <conio.h>
  6. #include <math.h>
  7. #include <locale.h>
  8. #include <iostream>
  9. #include <fstream>
  10. #define M_PI 3.14159265358979323846
  11. using namespace std;
  12.  
  13. #define BlocksPerGrid 1                         //Блоков в сетке
  14. #define ThreadsPerBlock 4
  15.  
  16. __global__ void Init_f(float *dev_f, int N,int dx_f, float h)
  17. {
  18.         int i=threadIdx.x+blockIdx.x*blockDim.x;
  19.         int j=threadIdx.y+blockIdx.y*blockDim.y;
  20.         if( ((i>(N/2-dx_f)) && (i<(N/2+dx_f)))  &&  (j>(N/2-dx_f)) && (j<(N/2+dx_f)) )
  21.                         dev_f[i*N+j] = h*h * 1.0 / (h*h);
  22.        
  23. }
  24.  
  25. __global__ void NormaF(float *dev_f, int N,float *norma_f)
  26. {
  27.         __shared__ float dev_norma_f;
  28.         int i=threadIdx.x+blockIdx.x*blockDim.x;
  29.         int j=threadIdx.y+blockIdx.y*blockDim.y;
  30.         while(i<N && j<N)
  31.         {
  32.                 dev_norma_f += fabs(dev_f[i*N+j]);
  33.         }
  34.         *norma_f=dev_norma_f;
  35. }
  36.  
  37. __global__ void Poisson2DSOR(float *dev_f,float *dev_psi,int N,float nev,float w,float *norma,float h)
  38. {
  39.         int i=threadIdx.x+blockIdx.x*blockDim.x;
  40.         int j=threadIdx.y+blockIdx.y*blockDim.y;
  41.         while( (i>1 && i<N-1) && (j>1 && j<N-1) )
  42.         {
  43.                 nev=dev_psi[(i - 1)*N+j]+dev_psi[(i + 1)*N+j]+dev_psi[i*N+(j - 1)]+dev_psi[i*N+(j + 1)]-4 * dev_psi[i*N+j] - 4*M_PI* h * h * dev_f[i*N+j]; //невязка
  44.                 *norma += fabs(nev);
  45.                 dev_psi[i*N+j]+=w*nev*0.25;
  46.         }
  47. }
  48. int main()
  49. {    
  50.         setlocale(LC_ALL,"Russian");
  51.  
  52.         int N =100;                                             // Число ячеек вдоль оси
  53.         int dx_f = N/5;
  54.         int i, iter=0;                                  //итеративные переменные
  55.         int maxit;                                              //максимальное кол-во итераций
  56.         int a, b;                                               //для определения четности шага(обход массива)
  57.         int memSIZE=sizeof(float)*N*N;
  58.         float eps, norma = 0.0, w = 1.0, nev=0.0;               //точность, норма, коэффициент, невязка
  59.         float h = 1.0 / (N - 1.0);                                          //шаг
  60.         float norma_f=0.0, sum = 0.0;
  61.         float *psi,*dev_psi; cudaMallocHost(&psi,memSIZE);//выделение памяти под psi на HOST
  62.         float *f,*dev_f; cudaMallocHost(&f,memSIZE);            //аналогично для f
  63.  
  64.         if(cudaSuccess != cudaMalloc((void**)&dev_f,memSIZE))
  65.         {
  66.                 printf("Error malloc dev_f!\n");getch();return 0;
  67.         }
  68.         if(cudaSuccess != cudaMalloc((void**)&dev_psi,memSIZE))
  69.         {
  70.                 printf("Error malloc dev_psi!\n");getch();return 0;
  71.         }
  72.  
  73.         dim3 gridsize=(BlocksPerGrid,BlocksPerGrid,1);
  74.         dim3 blocksize=(ThreadsPerBlock,ThreadsPerBlock,1);
  75.  
  76.         cudaStream_t stream[2];
  77.         for (int i=0;i<2;i++)
  78.         {
  79.                 if(cudaSuccess !=cudaStreamCreate(&stream[i]))
  80.                 {
  81.                         printf("Error streamCreate!\n");getch();return 0;
  82.                 }
  83.         }
  84.         if(cudaSuccess ==cudaMemcpyAsync( dev_f,
  85.                                          f,
  86.                                          memSIZE,
  87.                                          cudaMemcpyHostToDevice,
  88.                                          stream[0]))
  89.                 {
  90.                         printf("Success async copy dev_f!\n");getch();
  91.                 }
  92.         if(cudaSuccess ==cudaMemcpyAsync(dev_psi,
  93.                                         psi,
  94.                                         memSIZE,
  95.                                         cudaMemcpyHostToDevice,
  96.                                         stream[1]))
  97.                 {
  98.                         printf("Success acync copy dev_psi!\n");getch();
  99.                 }
  100.         if(cudaSuccess ==cudaDeviceSynchronize())
  101.         {
  102.                 printf("Success of cudaDeviceSynchronize!\n");getch();
  103.         }
  104.         Init_f <<<BlocksPerGrid,ThreadsPerBlock>>> (dev_f,N,dx_f,h);/*)*/
  105.          NormaF <<<BlocksPerGrid,ThreadsPerBlock>>> (dev_f,N,&norma_f);/*)*/
  106.         norma_f*=(4*M_PI*h*h);
  107.         printf("\nНорма f=%lf\n",norma_f);
  108.         maxit=100000;
  109.         eps=0.0001;
  110.  
  111.         double usl=eps*norma_f;
  112.         cout<<"USL="<<usl;
  113.         ofstream ofs("Test.txt");
  114.         cout<<"Вычисления....";
  115.  
  116.         for(iter = 0; iter < maxit; iter++)
  117.         {
  118.                 //cout<<"Итерация "<<iter<<endl;
  119.                 /*if(cudaSuccess !=*/ Poisson2DSOR <<<BlocksPerGrid,ThreadsPerBlock>>> (dev_f,dev_psi,N,nev,w,&norma,h);/*)*/
  120.                 /*{
  121.                         printf("Error Poisson2DSOR!\n");getch(); return 0;
  122.                 }*/
  123.                 cudaDeviceSynchronize();
  124.                 ofs<<iter<<"\t"<<norma<<"\t"<<endl;
  125.                 if(norma<usl)
  126.                         {
  127.                                 printf("\n\nВычисления закончены \n\n");
  128.                                 /*ofstream ofs1("FI.txt");
  129.                                 for(int i=0;i<N;i++)
  130.                                         for(int j=0;j<N;j++)
  131.                                                 {
  132.                                                         ofs1<<i<<" "<<j<<" "<<psi[i][j]*-1<<"\n";
  133.                                                 }*/
  134.                                 break;
  135.                         }
  136.     }
  137.         getch();
  138.  
  139.         return 0;
  140. }

Forums: