Привет, у меня есть математическая функция (многомерная, что означает, что есть индекс, который я передаю в C ++ - функцию, в которую я хочу вернуть одну математическую функцию. Например, скажем, у меня есть такая математическая функция:
f = Vector(x^2*y^2 / y^2 / x^2*z^2)
Я бы реализовал это так:
double myFunc(int function_index)
{
switch(function_index)
{
case 1:
return PNT[0]*PNT[0]*PNT[1]*PNT[1];
case 2:
return PNT[1]*PNT[1];
case 3:
return PNT[2]*PNT[2]*PNT[1]*PNT[1];
}
}
, тогда как PNT
определяется глобально так: double PNT[ NUM_COORDINATES ]
. Теперь я хочу реализовать производные каждой функции для каждой координатытаким образом генерируя производную матрицу (столбцы = координаты; строки = отдельные функции). Я уже написал свое ядро, которое работает до сих пор и которое вызывает myFunc ().
Проблема : Для расчетапроизводная математической подфункции i, касающейся координаты j, я бы использовал в последовательном режиме (например, на процессорах) следующий код (тогда как это упрощено, потому что обычно вы уменьшаете h, пока не достигнете определенной точности своей производной):
f0 = myFunc(i);
PNT[ j ] += h;
derivative = (myFunc(j)-f0)/h;
PNT[ j ] -= h;
Теперь, когда я хочу сделать это на GPU параллельно, проблемаподходит: что делать с PNT?Поскольку мне нужно увеличить определенные координаты на h, вычислить значение и затем уменьшить его снова, возникает проблема: как это сделать, не «мешая» другим потокам?Я не могу изменить PNT
, потому что другие потоки нуждаются в «исходной» точке для изменения их собственной координаты.
Вторая идея, которая у меня возникла, состояла в том, чтобы сохранить одну измененную точку для каждого потока, но я отбросил эту идею довольно быстро, потому что при использовании нескольких тысяч потоков параллельно это довольно плохо и, вероятно, медленно (возможно, вообще не реализуемо)из-за ограничений памяти) идея.
'ЗАКЛЮЧИТЕЛЬНОЕ' РЕШЕНИЕ Итак, как я это делаю в настоящее время, так это добавление значения 'add' во время выполнения (без его хранения где-либо) через макрос препроцессора к координате, обозначенной coord_index
.
#define X(n) ((coordinate_index == n) ? (PNT[n]+add) : PNT[n])
__device__ double myFunc(int function_index, int coordinate_index, double add)
{
//*// Example: f[i] = x[i]^3
return (X(function_index)*X(function_index)*X(function_index));
// */
}
Это работает довольно хорошо и быстро.При использовании производной матрицы с 10000 функциями и 10000 координатами это займет всего 0,5 сек.PNT
определяется либо глобально, либо как постоянная память, такая как __constant__ double PNT[ NUM_COORDINATES ];
, в зависимости от переменной препроцессора USE_CONST
.Строка return (X(function_index)*X(function_index)*X(function_index));
- это просто пример, где каждая подфункция выглядит одинаково, математически говоря:
f = Vector(x0^3 / x1^3 / ... / xN^3)
СЕЙЧАС БОЛЬШОЙ ПРОБЛЕМА :
myFunc
- это математическая функция, которую пользователь должен выполнять по своему усмотрению.Например, он мог бы также реализовать следующую математическую функцию:
f = Vector(x0^2*x1^2*...*xN^2 / x0^2*x1^2*...*xN^2 / ... / x0^2*x1^2*...*xN^2)
, таким образом, каждая функция выглядела одинаково.Вы, как программист, должны кодировать только один раз, и не в зависимости от реализованной математической функции.Поэтому, когда вышеуказанная функция реализуется в C ++, она выглядит следующим образом:
__device__ double myFunc(int function_index, int coordinate_index, double add)
{
double ret = 1.0;
for(int i = 0; i < NUM_COORDINATES; i++)
ret *= X(i)*X(i);
return ret;
}
И теперь доступ к памяти очень «странный» и плохо сказывается на производительности, поскольку каждому потоку необходим доступ к каждому элементуPNT
дважды.Конечно, в таком случае, когда каждая функция выглядит одинаково, я мог бы переписать полный алгоритм, который заключает в себе вызовы myFunc
, но, как я уже сказал: я не хочу кодировать в зависимости от пользовательской функции myFunc
...
Может кто-нибудь придумать, как решить эту проблему ??Спасибо!