开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 3D图形编程 > 查看源码
kernel.cu
资源名称：CUDAseg.zip　[点击查看]
上传用户：liuping58
上传日期：2022-06-05
资源大小：105k
文件大小：7k
源码类别：

3D图形编程

开发平台：

Visual C++

kernel.cu：源码内容
  1. #include <stdio.h>
  2. #include <cutil.h>
  4. #define DT  0.1
  6. #define max(x,y)    ((x>y) ? x : y )
  7. #define min(x,y)    ((x<y) ? x : y )
  9. #define INDEX(i,j,j_off)  (i +__mul24(j,j_off))
  11. #define BLOCKDIM_X  32
  12. #define BLOCKDIM_Y  4
  13. #define BLOCKDIM_Z  1
  15. __global__ void updatephi( float *d_phi, float *d_phi1, float *d_D, int NX, int NY, int NZ, float alpha, int pitch)
  16. {
  18. float dx,dy,dz;
  19. float dxplus, dyplus, dzplus, dxminus, dyminus, dzminus;
  20. float dxplusy, dxminusy, dxplusz, dxminusz, dyplusx, dyminusx, dyplusz, dyminusz, dzplusx, dzminusx, dzplusy, dzminusy;
  22. float gradphimax, gradphimin, nplusx, nplusy, nplusz, nminusx, nminusy, nminusz, curvature;
  23. float F, gradphi;
  25. int   indg, indg_h, indg0;
  26. int   i, j, k, ind, ind_h, halo, active;
  28. #define IOFF  1
  29. #define JOFF  (BLOCKDIM_X+2)
  30. #define KOFF  (BLOCKDIM_X+2)*(BLOCKDIM_Y+2)
  32. int NXM1 = NX-1; int NYM1 = NY-1; int NZM1 = NZ-1;
  34. __shared__ float s_data[3*(BLOCKDIM_X+2)*(BLOCKDIM_Y+2)];
  36. k    =  threadIdx.y*BLOCKDIM_X + threadIdx.x;
  37. halo = k < 2*(BLOCKDIM_X+BLOCKDIM_Y+2);
  39. if (halo) {
  40. if (threadIdx.y<2) {               // y-halos (coalesced)
  41. i = threadIdx.x;
  42. j = threadIdx.y*(BLOCKDIM_Y+1) - 1;
  43. }
  44. else {                             // x-halos (not coalesced)
  45. i = (k%2)*(BLOCKDIM_X+1) - 1;
  46. j =  k/2 - BLOCKDIM_X - 1;
  47. }
  49. ind_h  = INDEX(i+1,j+1,BLOCKDIM_X+2);
  51. i      = INDEX(i,blockIdx.x,BLOCKDIM_X);   // global indices
  52. j      = INDEX(j,blockIdx.y,BLOCKDIM_Y);
  53. indg_h = INDEX(i,j,pitch);
  55. halo   =  (i>=0) && (i<NX) && (j>=0) && (j<NY);
  56. }
  58. //
  59. // then set up indices for main block
  60. //
  62. i    = threadIdx.x;
  63. j    = threadIdx.y;
  64. ind  = INDEX(i+1,j+1,BLOCKDIM_X+2) ;
  66. i    = INDEX(i,blockIdx.x,BLOCKDIM_X);     // global indices
  67. j    = INDEX(j,blockIdx.y,BLOCKDIM_Y);
  68. indg = INDEX(i,j,pitch);
  70. active = (i<NX) && (j<NY);
  72. //
  73. // read initial plane of u1 array
  74. //
  76. if (active) s_data[ind+KOFF+KOFF] = d_phi1[indg]; if (halo) s_data[ind_h+KOFF+KOFF] = d_phi1[indg_h];
  80. for(int k=0;k<NZ;k++){
  82. if (active) { indg0 = indg; indg  = INDEX(indg,NY,pitch); s_data[ind-KOFF+KOFF] = s_data[ind+KOFF]; s_data[ind+KOFF]      = s_data[ind+KOFF+KOFF]; if (k<NZ-1) s_data[ind+KOFF+KOFF] = d_phi1[indg]; } if (halo) { indg_h = INDEX(indg_h,NY,pitch); s_data[ind_h-KOFF+KOFF] = s_data[ind_h+KOFF]; s_data[ind_h+KOFF]      = s_data[ind_h+KOFF+KOFF]; if (k<NZ-1) s_data[ind_h+KOFF+KOFF] = d_phi1[indg_h]; }
  84. if (active) {
  86. int ind2=ind+KOFF;
  88. if(i==0||i==NXM1){dx=0;} else {dx=(s_data[ind2+IOFF]-s_data[ind2-IOFF])/2;}
  89. if(j==0||j==NYM1){dy=0;} else {dy=(s_data[ind2-JOFF]-s_data[ind2+JOFF])/2;}
  90. if(k==0||k==NZM1){dz=0;} else {dz=(s_data[ind2+KOFF]-s_data[ind2-KOFF])/2;}
  92. if(i==NXM1){dxplus=0;}   else {dxplus =(s_data[ind2+IOFF]-s_data[ind2     ]);}
  93. if(j==0){dyplus=0;}  else {dyplus =(s_data[ind2-JOFF]-s_data[ind2     ]);}
  94. if(k==NZM1){dzplus=0;}   else {dzplus =(s_data[ind2+KOFF]-s_data[ind2     ]);}
  95. if(i==0){dxminus=0;}     else {dxminus=(s_data[ind2     ]-s_data[ind2-IOFF]);}
  96. if(j==NYM1){dyminus=0;}  else {dyminus=(s_data[ind2     ]-s_data[ind2+JOFF]);}
  97. if(k==0){dzminus=0;}     else {dzminus=(s_data[ind2     ]-s_data[ind2-KOFF]);}
  99. if(i==0||i==NXM1||j==0){dxplusy=0;}  else {dxplusy =(s_data[ind2-JOFF+IOFF]-s_data[ind2-JOFF-IOFF])/2;}
  100. if(i==0||i==NXM1||j==NYM1){dxminusy=0;}  else {dxminusy=(s_data[ind2+JOFF+IOFF]-s_data[ind2+JOFF-IOFF])/2;}
  101. if(i==0||i==NXM1||k==NZM1) {dxplusz=0;}  else {dxplusz =(s_data[ind2+KOFF+IOFF]-s_data[ind2+KOFF-IOFF])/2;}
  102. if(i==0||i==NXM1||k==0) {dxminusz=0;}  else {dxminusz=(s_data[ind2-KOFF+IOFF]-s_data[ind2-KOFF-IOFF])/2;}
  103. if(j==0||j==NYM1||i==NXM1){dyplusx=0;}   else {dyplusx =(s_data[ind2-JOFF+IOFF]-s_data[ind2+JOFF+IOFF])/2;}
  104. if(j==0||j==NYM1||i==0){dyminusx=0;}  else {dyminusx=(s_data[ind2-JOFF-IOFF]-s_data[ind2+JOFF-IOFF])/2;}
  105. if(j==0||j==NYM1||k==NZM1) {dyplusz=0;}  else {dyplusz =(s_data[ind2+KOFF-JOFF]-s_data[ind2+KOFF+JOFF])/2;}
  106. if(j==0||j==NYM1||k==0) {dyminusz=0;}  else {dyminusz=(s_data[ind2-KOFF-JOFF]-s_data[ind2-KOFF+JOFF])/2;}
  107. if(k==0||k==NZM1||i==NXM1) {dzplusx=0;}  else {dzplusx =(s_data[ind2+IOFF+KOFF]-s_data[ind2+IOFF-KOFF])/2;}
  108. if(k==0||k==NZM1||i==0) {dzminusx=0;}  else {dzminusx=(s_data[ind2-IOFF+KOFF]-s_data[ind2-IOFF-KOFF])/2;}
  109. if(k==0||k==NZM1||j==0) {dzplusy=0;}  else {dzplusy =(s_data[ind2-JOFF+KOFF]-s_data[ind2-JOFF-KOFF])/2;}
  110. if(k==0||k==NZM1||j==NYM1) {dzminusy=0;} else {dzminusy=(s_data[ind2+JOFF+KOFF]-s_data[ind2+JOFF-KOFF])/2;}
  113. gradphimax=sqrt((sqrt(max(dxplus,0)*max(dxplus,0)+max(-dxminus,0)*max(-dxminus,0)))*(sqrt(max(dxplus,0)*max(dxplus,0)+max(-dxminus,0)*max(-dxminus,0)))
  114. +(sqrt(max(dyplus,0)*max(dyplus,0)+max(-dyminus,0)*max(-dyminus,0)))*(sqrt(max(dyplus,0)*max(dyplus,0)+max(-dyminus,0)*max(-dyminus,0)))
  115. +(sqrt(max(dzplus,0)*max(dzplus,0)+max(-dzminus,0)*max(-dzminus,0)))*(sqrt(max(dzplus,0)*max(dzplus,0)+max(-dzminus,0)*max(-dzminus,0))));
  117. gradphimin=sqrt((sqrt(min(dxplus,0)*min(dxplus,0)+min(-dxminus,0)*min(-dxminus,0)))*(sqrt(min(dxplus,0)*min(dxplus,0)+min(-dxminus,0)*min(-dxminus,0)))
  118. +(sqrt(min(dyplus,0)*min(dyplus,0)+min(-dyminus,0)*min(-dyminus,0)))*(sqrt(min(dyplus,0)*min(dyplus,0)+min(-dyminus,0)*min(-dyminus,0)))
  119. +(sqrt(min(dzplus,0)*min(dzplus,0)+min(-dzminus,0)*min(-dzminus,0)))*(sqrt(min(dzplus,0)*min(dzplus,0)+min(-dzminus,0)*min(-dzminus,0))));
  121. nplusx = dxplus / sqrt(1.192092896e-07F + (dxplus*dxplus) + ((dyplusx + dy)*(dyplusx + dy)*0.25) + ((dzplusx + dz)*(dzplusx + dz)*0.25));
  122. nplusy = dyplus / sqrt(1.192092896e-07F + (dyplus*dyplus) + ((dxplusy + dx)*(dxplusy + dx)*0.25) + ((dzplusy + dz)*(dzplusy + dz)*0.25));
  123. nplusz = dzplus / sqrt(1.192092896e-07F + (dzplus*dzplus) + ((dxplusz + dz)*(dxplusz + dz)*0.25) + ((dyplusz + dy)*(dyplusz + dy)*0.25));
  125. nminusx=dxminus / sqrt(1.192092896e-07F + (dxminus*dxminus) + ((dyminusx + dy)*(dyminusx + dy)*0.25) + ((dzminusx + dz)*(dzminusx + dz)*0.25));
  126. nminusy=dyminus / sqrt(1.192092896e-07F + (dyminus*dyminus) + ((dxminusy + dx)*(dxminusy + dx)*0.25) + ((dzminusy + dz)*(dzminusy + dz)*0.25));
  127. nminusz=dzminus / sqrt(1.192092896e-07F + (dzminus*dzminus) + ((dxminusz + dz)*(dxminusz + dz)*0.25) + ((dyminusz + dy)*(dyminusz + dy)*0.25));
  129. curvature= ((nplusx-nminusx)+(nplusy-nminusy)+(nplusz-nminusz))/2;
  131. F = (-alpha * d_D[indg0]) + ((1-alpha) * curvature);
  132. if(F>0) {gradphi=gradphimax;} else {gradphi=gradphimin;}
  133. d_phi[indg0]=s_data[ind2] + (DT * F * gradphi);
  134. }
  136. __syncthreads();
  138. }
  140. }