Бывает, возникают задачи, в которых нужно использовать алгоритм обхода многомерного пространства в определённом порядке. Далее представлены некоторые рассуждения об этом вопросе. Рабочий язык Visual C++ 2008.
Рабочее многомерное пространство у нас будет представлять собой одномерный массив, для которого задана конфигурация измерений - массив размерности измерений. Для примера возьмём трёхмерное пространство.
const int DimSizes[] = { 675, 225, 2025 }; // конфигурация пространства - размерности измерений
const int TotalSize = DimSizes[0] * DimSizes[1] * DimSizes[2];
float *space = new float[TotalSize]; // само пространство
// ...
delete[] space;
Вариантом обхода будем считать систему вложенных циклов по измерениям. Например, прямой обход:
for(int i = 0; i < DimSizes[0]; ++i)
for(int j = 0; j < DimSizes[1]; ++j)
for(int k = 0; k < DimSizes[2]; ++k)
{
f(space[(i * DimSizes[0] + j) * DimSizes[1] + k], i, j, k);
}
или обратный обход:
for(int k = 0; k < DimSizes[2]; ++k)
for(int j = 0; j < DimSizes[1]; ++j)
for(int i = 0; i < DimSizes[0]; ++i)
{
f(space[(i * DimSizes[0] + j) * DimSizes[1] + k], i, j, k);
}
здесь функция f - это произвольный код, обрабатывающий очередной элемент пространства. Сигнатура функции:
void f(float &value, const int i, const int j, const int k);
Например, реализуем простое сложение всех элементов пространства с накоплением суммы в глобальной переменной. Замер производительности при прямом обходе даёт примерно 1,5 сек., а при обратном обходе примерно 9 сек. работы циклов. Этот результат неудивителен. Для иллюстрации рассмотрим 2D массив 3x3 { 1, 2, 3, 4, 5, 6, 7, 8, 9 }. При прямом обходе будет прочитана последовательность чисел, в порядке их расположения в массиве. При обратном обходе последовательность будет { 1, 4, 7, 2, 5, 8, 3, 6, 9 }. Во втором случае чтение элементов из памяти происходит "скачками" туда-сюда. На больших массивах это вынуждает процессор часто перезагружать свою кэш-память, что отрицательно сказывается на скорости работы программы. Однако мы оставим эту проблему в стороне. Будем считать описанные выше циклы эталоном производительности для данного порядка обхода и производительность всех полученных решений будем сравнивать с этим эталоном.
Внимательно посмотрим на циклы выше. Разница между ними заключается только в порядке вложения циклов и в порядке индексных аргументов функции f. Мы, разумеется, не может переставлять циклы прямо по ходу программы, поэтому такое решение не является универсальным. Классическим способом абстрагирования в таких случаях является косвенное обращение к переменным.
Определим вспомогательную структуру, которая будет описывать индекс по измерению. В ней, собственно, нужны лишь два поля: ссылка на переменную индекса и размер.
struct Index
{
int &index;
const int size;
Index(int &index, const int size) : index(index), size(size) {}
};
Тогда организация прямого обхода будет описываться следующим кодом:
int i, j, k;
Index x[] = { Index(i, DimSizes[0]), Index(j, DimSizes[1]), Index(k, DimSizes[2]) };
for(int &p = x[0].index = 0, P = x[0].size; p < P; ++p)
for(int &q = x[1].index = 0, Q = x[1].size; q < Q; ++q)
for(int &r = x[2].index = 0, R = x[2].size; r < R; ++r)
{
f(space[(i * DimSizes[0] + j) * DimSizes[1] + k], i, j, k);
}
Обратите внимание, что индексы заменились на косвенные только в заголовках, но не в теле циклов. Для организации обратного обхода достаточно поменять порядок индексных структур в массиве:
Index x[] = { Index(k, DimSizes[2]), Index(j, DimSizes[1]), Index(i, DimSizes[0]) };
Всё остальное остаётся без изменений. В этом решении не возникло никаких дополнительных вычислительных операций или операций с памятью. Ссылки являются полным эквивалентом исходной переменной. В результате производительность полученного решения оказывается идентичной исходному.
Теперь избавимся от ограничения на количество измерений. Помочь нам в этом может замена вложений циклов на вложение вызовов функции - рекурсия. При этом нам понадобится реорганизовать передачу индексов в функцию f.
Для начала упакуем все индексы в массив. Этот же массив будем использовать вместо списка аргументов функции f.
int i[3];
Index x[] = { Index(i[0], DimSizes[0]), Index(i[1], DimSizes[1]), Index(i[2], DimSizes[2]) };
for(int &p = x[0].index = 0, P = x[0].size; p < P; ++p)
for(int &q = x[1].index = 0, Q = x[1].size; q < Q; ++q)
for(int &r = x[2].index = 0, R = x[2].size; r < R; ++r)
{
f(space[(i[0] * DimSizes[0] + i[1]) * DimSizes[1] + i[2]], i);
}
void f(float &value, const int i[]);
Производительность, ожидаемо, не изменилась.
Далее полезно создать класс, который внутри себя будет содержать все вспомогательные структуры, а также производить настройку порядка обхода. Пусть конфигурация пространства задаётся в конструкторе объекта, а собственно порядок обхода в специальном методе, этот обход осуществляющем.
class SpaceRounding
{
const int dimQty;
const int *dimCfg;
int *idx;
float *space;
public:
SpaceRounding(const int dimQty, const int *dimCfg, float *space) :
dimQty(dimQty),
dimCfg(dimCfg),
idx(new int[dimQty]),
space(space)
{}
void Round(const int dimOrder[])
{
for(int i = dimOrder[0], &p = idx[i] = 0, P = dimCfg[i]; p < P; ++p)
for(int j = dimOrder[1], &q = idx[j] = 0, Q = dimCfg[j]; q < Q; ++q)
for(int k = dimOrder[2], &r = idx[k] = 0, R = dimCfg[k]; r < R; ++r)
{
f(space[(idx[0] * dimCfg[0] + idx[1]) * dimCfg[1] + idx[2]], idx);
}
}
~SpaceRounding()
{
delete[] idx;
}
};
Использование этого класса предельно просто:
SpaceRounding sr(3, DimSizes, space);
int order[] = { 0, 1, 2 }; // перечень индексов в порядке обхода
sr.Round(order);
Здесь мы пока не устраняли сами вложенные циклы. Стоит обратить внимание, что, во-первых, индексы как массив динамического размера, размещаются в куче, во-вторых, вместо вспомогательной структуры со ссылкой мы используем прямое обращение к массивам конфигурации пространства и порядка обхода. Оба этих решения обусловлены в общем случае неопределённым на этапе компиляции количеством измерений, с которыми будет работать объект. Каждое дополнительное обращение к массивам приводит к дополнительному обращению к памяти - операция разыменования указателя, - что резко снижает производительность решения примерно в 2 раза: до 2,5 сек. при прямом обходе и до 17 сек. при обратном обходе. Во вспомогательной структуре Index, которой мы пользовались выше, присутствует поле ссылочного типа. Это не позволяет создать массив подобных структур без немедленной инициализации каждого экземпляра, а инициализация в цикле происходит всегда после создания массива - не немедленно. Вариант с динамическими объектами типа Index возможен, но при этом в циклах всё равно будет обращение к объекту через допольнительное разыменование указателя, что и вызывает снижение производительности.
Использование шаблонного класса - это решение, обладающее как достоинствами, так и недостатками. Плюсом такого решения является возможность определить количество измерений на этапе компиляции и использовать это для восстановления производительности. Минусом - трудности с написанием пользовательских универсальных алгоритмов: они, в свою очередь, тоже должны быть шаблонными, что не всегда удобно. Как правило, на практике шаблонный класс оказывается достаточно хорошим решением.
Следующей значимой проблемой является формула расчёта смещения в массиве space. Эту формулу необходимо сделать итерационно вычисляемой для любого количества измерений. Нетрудно заметить, что при передаче индекса верхнего уровня вложенности вниз происходит домножение этого индекса на размер соответствующего измерения, после чего на нижнем уровне прибавляется следующий индекс. Динамическое вычисления смещения в циклах:
for(int i = dimOrder[0], &p = idx[i] = 0, P = dimCfg[i]; p < P; ++p)
for(int j = dimOrder[1], &q = idx[j] = 0, Q = dimCfg[j]; q < Q; ++q)
for(int k = dimOrder[2], &r = idx[k] = 0, R = dimCfg[k]; r < R; ++r)
{
int offset = 0;
for(int g = dimQty - 1, size = 1; g >= 0; size *= dimCfg[g], --g)
{
offset += idx[g] * size;
}
f(space[offset], idx);
}
И за него мы опять расплачиваемся снижением производительности: прямой обход работает уже примерно 4 сек., обратный - примерно 23,5 сек.. Это примерно трёхкратное ухудшение по сравнению с эталонным циклом.
На последнем шаге преобразуем вложенные циклы в рекурсивные вызовы. В основной процедуре будем запускать рекурсивную:
void Round(const int *dimOrder)
{
roundDim(dimQty, dimOrder);
}
Рекурсивная процедура в private-разделе класса:
void roundDim(const int dimCount, const int *dimOrder)
{
if(dimCount == 0)
{
int offset = 0;
for(int g = dimQty - 1, size = 1; g >= 0; size *= dimCfg[g], --g)
{
offset += idx[g] * size;
}
f(space[offset], idx);
}
else
{
for(int i = dimOrder[0], &p = idx[i] = 0, P = dimCfg[i]; p < P; ++p)
{
roundDim(dimCount - 1, dimOrder + 1);
}
}
}
Производительность немного понизилась: прямой обход - 5 сек., обратный - 25 сек. Таким образом, производительность ухудшилась на 1 сек. вне зависимости от порядка обхода - это результат включения между циклами дополнительных вызовов функции. Количество дополнительных вызвов зависит от количества измерений, но не от порядка обхода индексов.
Посмотрим внимательно на алгоритм расчёта положения элемента в пространстве. Он представляет собой попросту сумму индексов, домноженных на полный размер того подпространства, которое выбирается соответствующим индексом. Например, для пространства с конфигурацией 30x20x10 индекс (1, 2, 3) даст смещение 1 * (10 * 20) + 2 * 10 + 3 * 1 = 1 * 200 + 2 * 10 + 3 * 1 = 223. В этом выражении множители всегда постоянны, а индексы меняются каскадно - в порядке обхода. Мы можем воспользоваться свойством коммутативности сложения: не важно, в каком порядке идёт сложение индексов, а важно лишь то, чтобы соответствующий индекс был домножен на правильный размер. Тогда расчёт смещения можно будет распределить по уровням рекурсии, разгрузив внутренние циклы от дополнительной вычислительной операции. Для этого нам понадобится ввести дополнительный массив размеров подпространств соответствующих индексов.
class SpaceRounding
{
const int dimQty;
const int *dimCfg;
int *idx;
int *subspaceSizes;
float *space;
void roundDim(const int dimCount, const int *dimOrder, const int offset)
{
if(dimCount == 0)
{
f(space[offset], idx);
}
else
{
for(int i = dimOrder[0], &p = idx[i] = 0, P = dimCfg[i], O = subspaceSizes[i], o = p * O; p < P; ++p, o = p * O)
{
roundDim(dimCount - 1, dimOrder + 1, o + offset);
}
}
}
public:
SpaceRounding(const int dimQty, const int *dimCfg, float *space) :
dimQty(dimQty),
dimCfg(dimCfg),
idx(new int[dimQty]),
space(space),
subspaceSizes(new int[dimQty])
{
for(int i = dimQty - 1, size = 1; i >= 0; size *= dimCfg[i], --i)
{
subspaceSizes[i] = size;
}
}
void Round(const int *dimOrder)
{
roundDim(dimQty, dimOrder, 0);
}
~SpaceRounding()
{
delete[] subspaceSizes;
delete[] idx;
}
};
Такая модификация позволяет несколько повысить производительность: теперь прямой обход работает около 3,5 сек., а обратный - около 23 сек. Это в 2-3 раза медленнее, чем статические циклы - и это плата за гибкое решение.
Остаётся добавить классу вспомогательную клиентскую инфраструктуру. Сейчас класс обращается к конкретной функции f для обработки каждого элемента пространства. Нужно, чтобы эта функция была настраиваемой. Это означает введение ещё одного указателя внутрь циклов, что опять ухудшит производительность.
class SpaceRounding
{
public:
class Client
{
public:
virtual void Each(float &value, const int idx[]) = 0;
virtual void SetConfiguration(const int dimQty, const int dimCfg[])
{}
};
private:
const int dimQty;
const int *dimCfg;
int *idx;
int *subspaceSizes;
float *space;
Client &client;
void roundDim(const int dimCount, const int dimOrder[], const int offset)
{
if(dimCount == 0)
{
client.Each(space[offset], idx);
}
else
{
for(int i = dimOrder[0], &p = idx[i] = 0, P = dimCfg[i], O = subspaceSizes[i], o = p * O; p < P; ++p, o = p * O)
{
roundDim(dimCount - 1, dimOrder + 1, o + offset);
}
}
}
public:
SpaceRounding(const int dimQty, const int dimCfg[], float space[], Client &client) :
dimQty(dimQty),
dimCfg(dimCfg),
idx(new int[dimQty]),
space(space),
subspaceSizes(new int[dimQty]),
client(client)
{
client.SetConfiguration(dimQty, dimCfg);
for(int i = dimQty - 1, size = 1; i >= 0; size *= dimCfg[i], --i)
{
subspaceSizes[i] = size;
}
}
void Round(const int dimOrder[])
{
roundDim(dimQty, dimOrder, 0);
}
~SpaceRounding()
{
delete[] subspaceSizes;
delete[] idx;
}
};
class Client :
public SpaceRounding::Client
{
float sum;
public:
Client() :
sum(0.0f)
{}
virtual void Each(float &value, const int indexes[])
{
sum += value;
}
float Sum() const
{
return sum;
}
};
int main()
{
const int DimSizes[] = { 675, 225, 2025 };
const int TotalSize = DimSizes[0] * DimSizes[1] * DimSizes[2];
float *space = new float[TotalSize];
Client client;
SpaceRounding sr(3, DimSizes, space, client);
int order[] = { 2, 1, 0 };
sr.Round(order);
printf("%f\n", client.Sum());
delete[] space;
return 0;
}
Интерфейс для клиента имеет два метода. Один метод для обработки элемента пространства - его обязательно должен определить пользователь. Другой метод вспомогательный, для настройки под текущее пространство - его переопределение не обязательно. Производительность при этом падает вдвое от предыдущей и в 5-7 раз от исходной: примерно 11,5 сек. для прямого обхода и почти 50 сек. для обратного. Попытка замены объекта с виртуальной функцией простым указателем на функции даёт тот же самый результат по производительности. Важен сам факт обращения по указателю, а не форма его представления - с этой точки зрения работа с интерфейсом и возможность создать пользовательский объект удобнее для пользователя, чем простой указатель на функцию.
Из всех рассмотренных вариантов на практике во многих случаях достаточен второй - с фиксированным количеством измерений и косвенным обращением к индексам.
В заключении рассмотрим принципиальной иной подход, нацеленный на векторизацию вычислений (т.е. такое распараллеливание, при котором каждый элемент пространства обрабатывается независимо от других). Для этого нам понадобится общий векторизуемый цикл по всем элементам и две операции: разложение позиции в индексы и сборку из них новой позиции. Приём сводится к тому, что разложение позиции в индексы происходит с учётом заданного пользователем порядка обхода, а сборка новой позиции происходит так же, как это описывалось выше.
template<int DimQty>
class SpaceRounding
{
public:
class Client
{
public:
virtual void Each(float &value, const int idx[]) = 0;
virtual void SetConfiguration(const int dimQty, const int dimCfg[])
{}
};
private:
const int *dimCfg;
float *space;
Client &client;
int totalSize;
public:
SpaceRounding(const int dimCfg[DimQty], float space[], Client &client) :
dimCfg(dimCfg),
space(space),
client(client)
{
client.SetConfiguration(DimQty, dimCfg);
totalSize = 1;
for(int i = 0; i < DimQty; ++i)
{
totalSize *= dimCfg[i];
}
}
void Round(const int dimOrder[DimQty])
{
#pragma omp parallel for
for(int i = 0; i < totalSize; ++i)
{
int idx[DimQty];
int rest = i;
for(int j = DimQty - 1; j >= 0; --j)
{
int s = dimCfg[dimOrder[j]];
idx[dimOrder[j]] = rest % s;
rest /= s;
}
int offset = 0;
for(int j = DimQty - 1, size = 1; j >= 0; size *= dimCfg[j], --j)
{
offset += idx[j] * size;
}
client.Each(space[offset], idx);
}
}
};
Для векторизации важно, чтобы вспомогательный массив под индексы был собственным на каждой итерации, а не общим для всех итераций, как это было в предыдущих случаях. Создание серии таких массивов в куче - путь и заведомо очень медленный, и вносящий неопределённость в расход памяти при распараллеливании, и не позволяющий средствами, обеспечивающим параллельность выполнения кода, проанализировать потребность в памяти и эффективно распределить исполнение по рабочим нитям. Поэтому вспомогательный массив внутри цикла должен быть статического размера: заданного либо параметром шаблонного класса (как в примере), либо каким-нибудь разумным строгим ограничением сверху на количество поддерживаемых измерений.
Разумеется, задача векторизации обработки пространства - когда каждый элемент обрабатывается параллельно со всеми остальными - делает бессмысленным сам порядок последовательного обхода индексов из-за отсутствия самой последовательности. Поэтому вышеприведённое решение в таком виде является бессмысленным. Любопытно лишь отметить, что прямой обход здесь выполняется за примерно 21 сек., что существенно медленнее всех предыдущих решений, а обратный за примерно 50 секунд, что полностью соответствует предыдущему варианту.
Одним из вариантов развития этого метода векторизации являются алгоритмы преобразования пространств - например, транспонирование матриц, - когда каждый элемент исходного пространства нужно отобразить в другую позицию нового пространства. Описанное решение по получению исходного индекса, его матричного преобразования и вычисления новой позиции по новому индексу - это векторизуемая задача.
Другим вариантом развития этого метода является задача параллельной обработки подмножеств исходного пространства. Индексы элементов можно разделить на группы, описывающие дерево решений. Каждая группа описывает либо ветвь дерева из параллельно обрабатываемых подпространств, либо лист этого дерева - некое конечное продпространство. Группировка осмысленна, когда необходимо получить некий агрегируемый по группе результат, и агрегация результатов проходит по дереву от листьев к корню.
