Nested loop with variable ranges of the inner loop in CUDA C++
<p>I need to parallelize with CUDA a nested loop:</p>

<pre><code>for(int ix=0; ix&lt;Nx; ix++) {
   for(int iy=0; iy&lt;Ny[ix]; iy++) {
      SomeFunction(ix, iy);
      ...
</code></pre>

<p>where <code>Ny[]</code> is an array in CUDA global memory. This loop is called several times in a single run, <code>Nx</code> and the elements of <code>Ny[]</code> change in different calls, and they can be large (<code>Nx</code> from 0 to about 1 million, <code>Ny</code> from 0 to about 10000).
In principle I could use a CUDA kernel with threads arranged in a two-dimensional grid of size <code>Nx*max(Ny)</code> indexed by <code>ix</code> and <code>iy</code>, so the computational cost of the kernel would be <code>O[Nx*max(Ny)/Ncores]</code>. The problem is that in my case <code>max(Ny)</code> can be larger that the average value of <code>Ny</code> by orders of magnitudes. In some calls a few elements of <code>Ny</code> can be very large (say around 1000) while most of the other elements are very small. 
I would expect a good parallel implementation of the above loop to have a computational cost of <code>O[Nx*average(Ny)/Ncores]</code>, which in my case would be much smaller than <code>O[Nx*max(Ny)/Ncores]</code>, but I do not know I to do it. I have a vague idea that a good solution could be based on first sorting the elements of <code>Ny</code>.
Just for testing purpose, here is a piece of code that produces values of Ny that have a similar distribution as mine.</p>

<pre><code>#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;

int UpdateArray(int *array, int *Nx);

int main()
{
  int Nx_max=1000000;
  int *Ny;
  int Nx, i;

  Ny=(int*)malloc(Nx_max*sizeof(int));
  UpdateArray(Ny, &amp;Nx);

  for(i=0; i&lt;Nx; i++) {
    printf("%d\n", Ny[i]);
  }

  return 0;
}

int UpdateArray(int *array, int *Nx)
{
  int Nx0_min=500000, Nx0_max=1000000;
  int Nx1_min=50000, Nx1_max=100000;
  int Nx2_min=5000, Nx2_max=10000;

  int Ny0_min=1, Ny0_max=10;
  int Ny1_min=10, Ny1_max=100;
  int Ny2_min=100, Ny2_max=1000;

  int nx0, nx1, nx2, i, ix;

  nx0 = Nx0_min + rand()%(Nx0_max-Nx0_min); 
  for(i=0; i&lt;nx0; i++) {
    array[i] = Ny0_min + rand()%(Ny0_max-Ny0_min);
  }
  nx1 = Nx1_min + rand()%(Nx1_max-Nx1_min); 
  for(i=0; i&lt;nx1; i++) {
    ix = rand()%nx0;
    array[ix] = Ny1_min + rand()%(Ny1_max-Ny1_min);
  }
  nx2 = Nx2_min + rand()%(Nx2_max-Nx2_min); 
  for(i=0; i&lt;nx2; i++) {
    ix = rand()%nx0;
    array[ix] = Ny2_min + rand()%(Ny2_max-Ny2_min);
  }
  *Nx = nx0;

  return 0;
}
</code></pre>