Текст программы для графического процессора

Для того чтобы разделить части приложения, исполняемые на центральном и графическом процессорах, текст вычислительного ядра, предназначенного для графического процессора, размещён в отдельном файле, "Monte_Carlo_kernel.cu". Этот файл подключается к предыдущему директивой #include "Monte_Carlo_kernel.cu" на этапе компиляции.

/* Вычислительное ядро. Исполняется одновременно всеми 4096-ю вычислительными потоками. На то, что процедура является вычислительным ядром, указывает модификатор __global__ перед названием */

__global__ void Monte_Carlo_GPU(float *d_Random)

{

/* Идентификация конкретного вычислительного потока. Используются стандартные (автоматически существующие) системные переменные:

· blockDim.x – количество вычислительных потоков, приходящихся на одну «связку» (в нашем примере blockDim.x = 128);

· blockIdx.x – номер «связки», в нашем примере создаются 32 «связки»;

· threadIdx.x – номер вычислительного потока внутри «связки».

Таким образом, tid – номер конкретного вычислительного потока, изменяющийся через все «связки» от 0 до blockDim.x * blockIdx.x – 1. */

const int tid = blockDim.x * blockIdx.x + threadIdx.x;

/* Описание переменных для генератора Mersenne Twister */

int iState, iState1, iStateM, iOut;

unsigned int mti, mti1, mtiM, x;

unsigned int mt[MT_NN];

/* Модельные константы, описывающие взаимодействие нейтрона с пластиной: */

const float MC_h = 1.0f; // Толщина пластины

const float MC_sigma = 11.0f; // Полное макросечение взаимодействия

const float MC_sigma_S = 10.0f; // Макросечение рассеяния,
// обратное сред ней длине свободного пробега

float MC_s = MC_sigma_S / MC_sigma; // Вероятность рассеяния
// нейтрона

/* Загрузка инициализационных параметров генератора Mersenne Twister для данного конкретного вычислительного потока из device_MT, расположенного в видеопамяти; этот массив был создан процедурой seedMTGPU */

mt_parameters config = device_MT[tid];

/* Инициализация генератора */

mt[0] = config.seed;

for(iState = 1; iState < MT_NN; iState++)

mt[iState] = (1812433253U * (mt[iState - 1] ^ (mt[iState - 1] >> 30)) + iState) & MT_WMASK;

iState = 0;

mti1 = mt[0];

/* Начало генерации последовательности случайных чисел и моделирования прохождения нейтронов сквозь пластину */

float MC_x = 0.0f; // Координата нейтрона вдоль оси x

float MC_w = 1.0f / MC_s; // Начальный вес нейтрона

float MC_m = 1.0f; // Косинус угла отклонения нейтрона

float MC_L; // Длина свободного пробега

/* Начало цикла, каждый шаг которого соответствует одному столкновению нейтрона внутри пластины */

do

{

/* Пересчёт веса нейтрона в результате столкновения. Размещение этой операции в начале цикла исключает использование дополнительного условного оператора, так как возвращение к началу цикла подразумевает столкновение. Для того, чтобы в ходе первой итерации вес нейтрона равнялся 1, перед входом в цикл MC_w = 1 / MC_s */

MC_w *= MC_s;

/* Первый вызов генератора Mersenne Twister, для получения случайной длины свободного пробега MC_L */

iState1 = iState + 1;

iStateM = iState + MT_MM;

if(iState1 >= MT_NN) iState1 -= MT_NN;

if(iStateM >= MT_NN) iStateM -= MT_NN;

mti = mti1;

mti1 = mt[iState1];

mtiM = mt[iStateM];

x = (mti & MT_UMASK) | (mti1 & MT_LMASK);

x = mtiM ^ (x >> 1) ^ ((x & 1) ? config.matrix_a : 0);

mt[iState] = x;

iState = iState1;

x ^= (x >> MT_SHIFT0);

x ^= (x << MT_SHIFTB) & config.mask_b;

x ^= (x << MT_SHIFTC) & config.mask_c;

x ^= (x >> MT_SHIFT1);

/* Преобразование случайного числа x к диапазону (0; 1) и получение случайной длины свободного пробега по формуле (2.5) */

MC_L = -__logf(((float)x + 1.0f) / 4294967296.0f) / MC_sigma;

// Определение нового положения нейтрона

MC_x += MC_L*MC_m;

/* Второй вызов генератора Mersenne Twister, для получения случайного отклонения новой траектории нейтрона от оси x */

iState1 = iState + 1;

iStateM = iState + MT_MM;

if(iState1 >= MT_NN) iState1 -= MT_NN;

if(iStateM >= MT_NN) iStateM -= MT_NN;

mti = mti1;

mti1 = mt[iState1];

mtiM = mt[iStateM];

x = (mti & MT_UMASK) | (mti1 & MT_LMASK);

x = mtiM ^ (x >> 1) ^ ((x & 1) ? config.matrix_a : 0);

mt[iState] = x;

iState = iState1;

//Tempering transformation

x ^= (x >> MT_SHIFT0);

x ^= (x << MT_SHIFTB) & config.mask_b;

x ^= (x << MT_SHIFTC) & config.mask_c;

x ^= (x >> MT_SHIFT1);

/* Преобразование случайного числа x к диапазону (0; 1) и получение случайного косинуса угла отклонения траектории нейтрона от оси x по формуле (2.7) */

MC_m = 2.0f*(((float)x + 1.0f) / 4294967296.0f) - 1.0f;

/* Проверка условия выхода нейтрона из пластины и завершение цикла при выполнении этого условия */

} while ((MC_x < 0.0f) || (MC_x > MC_h))

/* Запись результатов моделирования в видеопамять. В данном примере каждый генератор пишет два числа: 0.0, если нейтрон отразился от пластины либо 1.0, если нейтрон прошёл пластину; второе число – окончательный вес нейтрона; d_Random – указатель на начало массива результатов в видеопамяти*/

int MC_address = 2*tid;

if (MC_x < 0.0f) d_Random[MC_address] = 0.0f;

if (MC_x > MC_h) d_Random[MC_address] = 1.0f;

d_Random[MC_address+1] = MC_w;

} // Завершение вычислительного ядра