title: Cuda version not working while serial working
body: <p>my following minimalist Cuda code returns an incorrect result (all polygons have 0 vertices at the end) while the same code running in serial in C++ is working well. The problem is embarrassingly parallel : no communication, no syncthreads etc., and the Cuda memory allocations are sucessful. Even my dummy variable that stores the content of the input array for debug purpose is 0 for the Cuda version. There is no access out of bounds since my arrays are largely large enough. Replacing the memcpy by a loop in Cuda doesn't change anything.<br>
  I really don't understand what happens... any idea ?  Thanks!</p>
  
  <p>Cuda code:</p>
  
  <pre><code>    #include &lt;stdio.h&gt;
      #include &lt;iostream&gt;
      #include &lt;stdlib.h&gt;
      #include &lt;cuda.h&gt;
  
      class Point2D {
       public:
       __device__ Point2D(double xx=0, double yy=0):x(xx),y(yy){};
       double x, y;
      };
  
      __device__ double dot(const Point2D &amp;A, const Point2D &amp;B) {
       return A.x*B.x + A.y*B.y;
      }
      __device__ Point2D operator*(double a, const Point2D &amp;P) {
       return Point2D(a*P.x, a*P.y);
      }
      __device__ Point2D operator+(Point2D A, const Point2D &amp;B) {
       return Point2D(A.x + B.x, A.y + B.y);
      }
      __device__ Point2D operator-(Point2D A, const Point2D &amp;B) {
       return Point2D(A.x - B.x, A.y - B.y);
      }
      __device__ Point2D inter(const Point2D &amp;A, const Point2D &amp;B, const Point2D &amp;C, const Point2D &amp;D) { //intersects AB by *the mediator* of CD 
        Point2D M = 0.5*(C+D);
        return A - (dot(A-M, D-C)/dot(B-A, D-C)) * (B-A);
      }
  
      class Polygon {
      public:
        __device__ Polygon():nbpts(0){};
        __device__ void addPts(Point2D pt) {
          pts[nbpts] = pt;
          nbpts++;
        }; 
        __device__ Polygon&amp; operator=(const Polygon&amp; rhs) {
          nbpts = rhs.nbpts;
          dummy = rhs.dummy;
          memcpy(pts, rhs.pts, nbpts*sizeof(Point2D));
          return *this;
        }
        __device__ void cut(const Point2D &amp;inside_pt, const Point2D &amp;outside_pt) {
  
          int new_nbpts = 0;
          Point2D newpts[128];
          Point2D AB(outside_pt-inside_pt);
          Point2D M(0.5*(outside_pt+inside_pt));
          double ABM = dot(AB, M);
  
          Point2D S = pts[nbpts-1];
  
          for (int i=0; i&lt;nbpts; i++) {
  
            Point2D E = pts[i];
  
            double ddot = -ABM + dot(AB, E);
            if (ddot&lt;0) { // E inside clip edge
              double ddot2 = -ABM + dot(AB, S);
              if (ddot2&gt;0) {
                 newpts[new_nbpts] = inter(S,E, inside_pt, outside_pt);
                 new_nbpts++;
              }
              newpts[new_nbpts] = E;
              new_nbpts++;
            } else {
              double ddot2 = -ABM + dot(AB, S);
              if (ddot2&lt;0) {
                 newpts[new_nbpts] = inter(S,E, inside_pt, outside_pt);
                 new_nbpts++;
              }       
            }
            S = E;
          }
  
          memcpy(pts, newpts, min(128, new_nbpts)*sizeof(Point2D));
          nbpts = new_nbpts;
        }
  
      //private:
       Point2D pts[128];
       int nbpts;
       float dummy;
      };
  
  
      __global__ void cut_poly(float *a, Polygon* polygons, int N)
      {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx&gt;=N/2) return;
  
        Polygon pol;
        pol.addPts(Point2D(0.,0.));
        pol.addPts(Point2D(1.,0.));
        pol.addPts(Point2D(1.,1.));
        pol.addPts(Point2D(0.,1.));
  
        Point2D curPt(a[2*idx], a[2*idx+1]);
  
        for (int i=0; i&lt;N/2; i++) {
          Point2D other_pt(a[2*i], a[2*i+1]);
          pol.cut(curPt, other_pt);
        }
        pol.dummy = a[idx];
  
        polygons[idx] = pol;
      }
  
  
  
      int main(int argc, unsigned char* argv[])
      {
  
        const int N = 100; 
        float a_h[N], *a_d; 
        Polygon p_h[N/2], *p_d;
  
        size_t size = N * sizeof(float);
        size_t size_pol = N/2 * sizeof(Polygon);
  
        cudaError_t err  = cudaMalloc((void **) &amp;a_d, size);   
        cudaError_t err2 = cudaMalloc((void **) &amp;p_d, size_pol);  
  
        for (int i=0; i&lt;N; i++) a_h[i] = (float)(rand()%1000)*0.001;
        cudaMemcpy(a_d, a_h, size, cudaMemcpyHostToDevice);
  
        int block_size = 4;
        int n_blocks = N/block_size + (N%block_size == 0 ? 0:1);
        cut_poly &lt;&lt;&lt; n_blocks, block_size &gt;&gt;&gt; (a_d, p_d, N);
  
        cudaMemcpy(a_h, a_d, sizeof(float)*N, cudaMemcpyDeviceToHost);
        cudaMemcpy(p_h, p_d, sizeof(Polygon)*N/2, cudaMemcpyDeviceToHost);
  
        for (int i=0; i&lt;N/2; i++)
         printf("%f \t %f \t %u\n", a_h[i], p_h[i].dummy, p_h[i].nbpts);
  
        cudaFree(a_d);
        cudaFree(p_d);
  
  
          return 0;
      }
  </code></pre>
  
  <p>Same code in C++ that works properly:</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;iostream&gt;
  #include &lt;stdlib.h&gt;
  
  class Point2D {
   public:
   Point2D(double xx=0, double yy=0):x(xx),y(yy){};
   double x, y;
  };
  
  double dot(const Point2D &amp;A, const Point2D &amp;B) {
   return A.x*B.x + A.y*B.y;
  }
  Point2D operator*(double a, const Point2D &amp;P) {
   return Point2D(a*P.x, a*P.y);
  }
  Point2D operator+(Point2D A, const Point2D &amp;B) {
   return Point2D(A.x + B.x, A.y + B.y);
  }
  Point2D operator-(Point2D A, const Point2D &amp;B) {
   return Point2D(A.x - B.x, A.y - B.y);
  }
  Point2D inter(const Point2D &amp;A, const Point2D &amp;B, const Point2D &amp;C, const Point2D &amp;D) { //intersects AB by *the mediator* of CD 
    Point2D M = 0.5*(C+D);
    return A - (dot(A-M, D-C)/dot(B-A, D-C)) * (B-A);
  }
  
  class Polygon {
  public:
    Polygon():nbpts(0){};
    void addPts(Point2D pt) {
      pts[nbpts] = pt;
      nbpts++;
    }; 
    Polygon&amp; operator=(const Polygon&amp; rhs) {
      nbpts = rhs.nbpts;
      dummy = rhs.dummy;
      memcpy(pts, rhs.pts, nbpts*sizeof(Point2D));
      return *this;
    }
    void cut(const Point2D &amp;inside_pt, const Point2D &amp;outside_pt) {
  
      int new_nbpts = 0;
      Point2D newpts[128];
      Point2D AB(outside_pt-inside_pt);
      Point2D M(0.5*(outside_pt+inside_pt));
      double ABM = dot(AB, M);
  
      Point2D S = pts[nbpts-1];
  
      for (int i=0; i&lt;nbpts; i++) {
  
        Point2D E = pts[i];
  
        double ddot = -ABM + dot(AB, E);
        if (ddot&lt;0) { // E inside clip edge
          double ddot2 = -ABM + dot(AB, S);
          if (ddot2&gt;0) {
             newpts[new_nbpts] = inter(S,E, inside_pt, outside_pt);
             new_nbpts++;
          }
          newpts[new_nbpts] = E;
          new_nbpts++;
        } else {
          double ddot2 = -ABM + dot(AB, S);
          if (ddot2&lt;0) {
             newpts[new_nbpts] = inter(S,E, inside_pt, outside_pt);
             new_nbpts++;
          }
        }
          S = E;
      }
  
      memcpy(pts, newpts, std::min(128, new_nbpts)*sizeof(Point2D));
      /*for (int i=0; i&lt;128; i++) {
        pts[i] = newpts[i];
      }*/
      nbpts = new_nbpts;
    }
  
  //private:
   Point2D pts[128];
   int nbpts;
   float dummy;
  };
  
  
  void cut_poly(int idx, float *a, Polygon* polygons, int N)
  {
    if (idx&gt;=N/2) return;
  
    Polygon pol;
    pol.addPts(Point2D(0.,0.));
    pol.addPts(Point2D(1.,0.));
    pol.addPts(Point2D(1.,1.));
    pol.addPts(Point2D(0.,1.));
  
    Point2D curPt(a[2*idx], a[2*idx+1]);
  
    for (int i=0; i&lt;N/2; i++) {
      if (idx==i) continue;
      Point2D other_pt(a[2*i], a[2*i+1]);
      pol.cut(curPt, other_pt);
    }
    pol.dummy = a[idx];
  
    polygons[idx] = pol;
  }
  
  
  
  int main(int argc, unsigned char* argv[])
  {
  
    const int N = 100;  // Number of elements in arrays
    float a_h[N], *a_d;  // Pointer to host &amp; device arrays
    Polygon p_h[N/2], *p_d;
  
    for (int i=0; i&lt;N; i++) a_h[i] = (float)(rand()%1000)*0.001;
  
    for (int idx=0; idx&lt;N; idx++)
      cut_poly(idx, a_h, p_h, N);
  
    for (int i=0; i&lt;N/2; i++)
     printf("%f \t %f \t %u\n", a_h[i], p_h[i].dummy, p_h[i].nbpts);
  
     return 0;
  }
  </code></pre>