Самостоятельная работа 2: Оптимизация прикладных программ для Intel Xeon Phi с использованием Intel C/C++ Compiler. Векторизация

Более сложные примеры векторизации

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

Рассмотрим две типичные операции, часто встречающиеся в вычислительных процедурах финансовой математики.

При вычислении справедливой цены опциона Европейского типа в популярной двухфакторной модели HJM возникает код следующего вида:

double* arr = (double*) malloc(size*sizeof(double));
arr[0] = 1.1;
for (i = 0; i < size - 1; i++)
{
    double c0 = sigma1*pow(arr[i], alpha);
    double c1 = sigma2*pow(arr[i], beta);
    arr[i+1] = (c0*z1[i] + c1*z2[i]) + 1;
}

Очевидно, все элементы массива arr зависимы между собой и векторизация цикла for невозможна. Однако в данном случае возможна векторизация вычисления powv внутри цикла. Так как компилятор применяет автоматическую векторизацию лишь для циклов, необходимо трансформировать две скалярные операции в цикл из двух итераций. Кроме того, необходимо использовать #pragma vector always, в противном случае цикл, вероятно, не будет векторизован, так как компилятор сочтет его слишком коротким. Преобразованный код имеет вид:

double* arr = (double*) malloc(size*sizeof(double));
arr[0] = 1.1;
for (i = 0; i < size - 1; i++)
{
    double c[2], sigma[2] = {sigma1, sigma2};
    double power[2] = {alpha, beta};
    #pragma vector always
    for (int k=0; k < 2; k++) 
        c[k] = sigma[k] * pow(arr[i], power[k]);
    arr[i+1] = (c0*z1[i] + c1*z2[i]) + 1;
}

Цикл с вычислением функции pow будет векторизован и превратится в один вызов векторной реализации pow из SVML.

В той же задаче далее необходимо произвести суммирование по всем путям Монте-Карло, для которых выгода положительна, и возникает код следующего вида:

for (int i = 0; i < size; i++) 
{
    double s = coeff * exp(arr[i]);
    double payoff = s - K;
    if (payoff > 0.0)
        sum = sum + payoff;
}

Хотя условные операции в целом нежелательны для векторизации, при простой структуре условия они не являются помехой. Данный цикл векторизуется с использованием векторного сравнения нескольких пар чисел с плавающей запятой. Другим возможным вариантом написания цикла была бы замена условия на явное взятие максимума:

for (int i = 0; i < size; i++) 
{
    double s = coeff * exp(arr[i]);
    double payoff = s - K;
    sum = sum + std::max(payoff, 0.0);
}

В данном случае цикл также векторизуется.

При большой длине цикла size (в данной задаче это количество путей Монте-Карло и, поэтому, является весьма большим) также имеет смысл вычислять значение экспоненты для всех аргументов с использованием VML. При этом код приобретает вид (в предположении, что содержимое массива arr после завершения цикла не используется):

vdExp(n, arr, arr);
for (int i = 0; i < size; i++) 
    sum = sum + std::max(coeff * arr[i] - K, 0.0);

Дополнительные задания

1. Для рассматриваемого в разделе 2 примера сравните время работы исходной (скалярной) и векторизованной версий на центральном процессоре и сопроцессоре Intel Xeon Phi. Сравните время работы версий с использованием различных средств векторизации. Объясните полученный результат.
2. Реализуйте умножение матрицы на вектор. Воспользуйтесь отчетом о векторизации и проверьте, происходит ли векторизация вашей реализации? При необходимости внесите изменения для обеспечения векторизации. Зависит ли векторизация от того, как хранится матрица: по строкам или по столбцам?
3. Проведите исследование, аналогичное предыдущему заданию, для операции вычисления матричного произведения по определению.
4. Рассмотрите пример из раздела 5. Экспериментально определите минимальное количество итераций цикла, при котором предварительное вычисление всех экспонент с использованием VML становится более эффективным по сравнению с использованием SVML.