Самостоятельная работа 1: Оптимизация вычислительно трудоемкого программного модуля для архитектуры Intel Xeon Phi. Метод Монте-Карло

Оптимизация 2: отказ от использования дублирующих массивов

Основной недостаток текущей версии программы состоит в том, что мы не можем использовать все возможности сопроцессора – число потоков, на которых мы можем запустить программу, ограничено объемом памяти сопроцессора. И если при работе на CPU использование дублирующих массивов – вполне нормальная и часто применяемая практика, то при переходе на Intel Xeon Phi эта техника не всегда себя оправдывает.

В нашем случае дублирование выполнено для двух типов результатов: общей фотонной карты траекторий и фотонных карт для каждого из детекторов. Объем операций доступа к этим массивам существенно различен. Общая фотонная карта траекторий обновляется каждым фотоном, в то время как фотонные карты детекторов обновляются реже, так как далеко не каждый фотон попадает в тот или иной детектор. С другой стороны, суммарный размер фотонных карт для всех детекторов на порядок превосходит размер общей фотонной карты (в примере это 40 и 4 МБ соответственно).

Отказ от дублирования общей фотонной карты с одной стороны приводит к существенным затратам на синхронизацию, а с другой – выигрыш по памяти в этом случае не так велик.

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

Заметим, что в данной версии мы отказываемся от дублирования еще одного массива – массива сигналов на детекторах (output->weightInDetector). Целесообразность такого решения предлагается оценить слушателю в рамках дополнительных заданий.

Для реализации предложенной идеи прежде всего следует избавиться от дублирования данных. Основные изменения коснутся функции параллельного запуска процесса моделирования LaunchOMP() (./Opt2/xmcmlLauncher/launcher_omp.cpp):

void LaunchOMP(InputInfo* input, OutputInfo* output,
    MCG59* randomGenerator, int numThreads)
{
    …
    #pragma omp parallel 
    {
        int threadId = omp_get_thread_num();
        InitThreadOutput(&(threadOutputs[threadId]),
            output);

        PhotonTrjectory trajectory;

        for (uint64 i = threadId;
            i < input->numberOfPhotons; i += numThreads)
        {
            ComputePhoton(specularReflectance, input,
                &(threadOutputs[threadId]),
                &(randomGenerator[threadId]), &trajectory);
        }
    }

    output->specularReflectance = specularReflectance;
    for (int i = 0; i < numThreads; ++i)
    {
        for (int j = 0; j < output->gridSize; ++j)
        {
            output->commonTrajectory[j] +=
                threadOutputs[i].commonTrajectory[j];
        }
    }

    for (int i = 0; i < numThreads; ++i)
    {
        FreeThreadOutput(threadOutputs + i);
    }
    delete[] threadOutputs;
}

Финальное суммирование осталось только для общей фотонной карты, вместо функций InitOutput() и FreeOutput() используются (./Opt2/xmcmlLauncher/launcher_omp.cpp):

void InitThreadOutput(OutputInfo* dst, OutputInfo* src)
{
    dst->gridSize = src->gridSize;
    dst->detectorTrajectory = src->detectorTrajectory;
    dst->numberOfDetectors = src->numberOfDetectors;
    dst->numberOfPhotons = src->numberOfPhotons;
    dst->specularReflectance = src->specularReflectance;
    dst->weightInDetector = src->weightInDetector;
    
    dst->commonTrajectory = new uint64[dst->gridSize];
    memset(dst->commonTrajectory, 0, dst->gridSize*sizeof(uint64));
}

void FreeThreadOutput(OutputInfo* output)
{
    if (output->commonTrajectory != NULL)
    {
        delete[] output->commonTrajectory;
    }
}

Избавившись от дублирования, необходимо позаботиться о синхронизации. Для этого прежде нужно изменить функцию UpdateDetectorTrajectory() (./Opt2/xmcml/mcml_kernel.cpp):

void UpdateDetectorTrajectory(OutputInfo* output, Area* area, PhotonTrjectory* trajectory, int detectorId)
{
    int index;
    int ny = area->partitionNumber.y;
    int nz = area->partitionNumber.z;

    uint64* detectorTrajectory =
        output->detectorTrajectory[detectorId].trajectory;

    for (int i = 0; i < trajectory->size; ++i)
    {
        index = trajectory->x[i]*ny*nz +
            trajectory->y[i]*nz + trajectory->z[i];

        #pragma omp atomic
        ++(detectorTrajectory[index]);
    }

    #pragma omp atomic
    ++(output->detectorTrajectory[detectorId].
        numberOfPhotons);
}

Также изменению подвергнется функция UpdateWeightInDetector() (./Opt2/xmcml/mcml_kernel.cpp):

void UpdateWeightInDetector(OutputInfo* output,
    double photonWeight, int detectorId)
{
    #pragma omp atomic
    output->weightInDetector[detectorId] += photonWeight; 
}

Время работы программы на центральном процессоре после применения оптимизации практически не изменилось. Время работы программы на сопроцессоре в 120 потоков уменьшилось до 36 секунд. Более того, появилась возможность запуска программы в 240 потоков. Однако время моделирования в этом случае составило 40 секунд. Итого, данная оптимизация позволила сократить время вычислений на сопроцессоре еще на 35%.