Как GPU помогает улучшить итерационные вопросы?

Ваша программа может быть легко распараллелена в цикле while (k<10000000). На самом деле, поскольку условие завершения программы - это неизвестное количество итераций (чтобы достичь 10M хороших наборов), вы можете по существу удалить весь этот код, который вы показали в ядре, и запустить его как есть, с несколькими незначительными изменениями.

#include <curand.h>
#include <curand_kernel.h>

__constant__ double c = 0.24;
__device__ volatile unsigned int k = 0;
#define SCALE 32767.0
#define NUM_GOOD 10000000

__device__ int is_good(double f, double g, double h){
  if (....){
    ...
    return 1;
  }
  return 0;
}

__global__ void initCurand(curandState *state, unsigned long seed){
  int idx = threadIdx.x + blockIdx.x*blockDim.x;
  curand_init(seed, idx, 0, &state[idx]);
}

__global__ void mykernel(curandState *devStates, double *good_f, double *good_g, double *good_h){
  int idx = threadIdx.x + blockDim.x*blockIdx.x;
  double f0, g0, h0, f1, g1, h1;
  curandState localState = devStates[idx];
  while (k<NUM_GOOD){
    // assuming you wanted independent starting values for f, g, h
    f0 = (double)(curand_uniform(&localState)*SCALE);
    g0 = (double)(curand_uniform(&localState)*SCALE);
    h0 = (double)(curand_uniform(&localState)*SCALE);
    for (int i = 0; i< 450000000; i++){
      f1 = g0 + 2*h0 + c;
      g1 = 0.5*f1 - h0;
      h1 = (f1+g1 )/2;
      f0 = g1 + 2*h1 + c;
      g0 = 0.5*f0 - h1;
      h0 = (f0+g0 )/2;}
    if (is_good(f1, g1, h1))
    {
      unsigned int next =  atomicAdd(&k, 1);
      if (next<NUM_GOOD){
        good_f[next] = f1;
        good_g[next] = g1;
        good_h[next] = h1;}
    }
  }
}

Вышеприведенный код является лишь контуром, могут быть некоторые ошибки и, очевидно, не все здесь определено.

Вы можете играть с количеством фактических потоков, которые вы запускаете, чтобы увидеть, что работает быстрее всего. Все запущенные потоки будут работать над заполнением "хорошего" стека, пока он не будет заполнен. Затем каждый поток обнаруживает, что стек заполнен и завершен.

EDIT: ответьте на некоторые вопросы ниже:

кажется, что "int idx = threadIdx.x + blockDim.x * blockIdx.x;" это материал GPU, я думаю, что это связано с потоковой обработкой в GPU, так что это важно для программирования GPU?

Да, эти переменные, такие как threadIdx.x являются "встроенными" переменными в CUDA, которые позволяют каждому потоку делать что-то другое (начинаются с разных случайных значений, в этом случае).

Во-вторых, весь код, который вы указали, выглядит как обычный C++ код. Но вы поставили "критический раздел GPU", так что любой специальный синтаксис, который мне нужно использовать в этом разделе, или он похож на обычный C++ код?

Да, большая часть кода ядра CUDA может быть обычным C++ кодом, часто подобным тому, что вы могли бы написать, чтобы делать то же самое на процессоре. В этом случае я упомянул критический раздел и связал пример, но, подумав об этом, критический раздел (в этом случае использовался для ограничения доступа к области данных, чтобы потоки GPU не переключались друг на друга по мере их обновления "хорошие" значения) здесь избыток. Нужно только использовать атомную операцию для резервирования "пятна" в стеке для каждого потока, который хочет заполнить хорошее значение. Я изменил код соответственно.

В соответствии с

 while (k<10000000)

вы пытаетесь найти 10M good {f, h, g}.

В вашем однопоточном коде процессора вы находите их один за другим, а в GPU легко запускать тысячи потоков, чтобы найти удовлетворительные результаты параллельно, пока общее число не достигнет 10M.

Для итеративной части соединения вам все равно нужно их вычислить традиционным способом. Но вы по-прежнему можете улучшить производительность в этой части, упростив уравнения как

f(n+1) = 1   *g(n) + 2*h(n) +      c;
g(n+1) = 0.5 *g(n)          +  0.5*c;
h(n+1) = 0.75*g(n) + 1*h(n) + 0.75*c;

Матрица преобразования A для вектора [f,g,h,c]' (в коде MATLAB)

A = [ 0 1 2 1 ; 0 .5 0 .5; 0 .75 1 .75 ; 0 0 0 0];

Тогда имеем [f,g,h,c]'{n}=A^n * [f,g,h,c]'{0}. Вы найдете, что A^n сходится к [0 3 2 3; 0 0 0 0; 0 1.5 1 1.5; 0 0 0 0] [0 3 2 3; 0 0 0 0; 0 1.5 1 1.5; 0 0 0 0] [0 3 2 3; 0 0 0 0; 0 1.5 1 1.5; 0 0 0 0] в несколько итераций.