Эффективная программа транспонирования матрицы в кеше?
Таким образом, очевидный способ транспонировать матрицу - это использовать:
for( int i = 0; i < n; i++ )
for( int j = 0; j < n; j++ )
destination[j+i*n] = source[i+j*n];
но я хочу что-то, что будет использовать преимущества локальности и блокировки кеша. Я искал его и не могу найти код, который бы это делал, но мне сказали, что это должно быть очень простым изменением оригинала. Есть идеи?
Изменить: у меня есть матрица 2000x2000, и я хочу знать, как я могу изменить код, используя два for циклы, в основном разбивая матрицу на блоки, которые я перемещаю по отдельности, например, блоки 2x2 или 40x40, и вижу, какой размер блока наиболее эффективен.
Edit2: матрицы хранятся в главном порядке столбцов, то есть для матрицы
a1 a2
a3 a4
хранится как a1 a3 a2 a4,
7 ответов
Вы, вероятно, захотите четыре цикла - два для итерации по блокам, а затем еще два для выполнения транспонирования-копии одного блока. Предполагая для простоты размер блока, который делит размер матрицы, я думаю, что-то вроде этого, хотя я бы хотел нарисовать несколько картинок на обороте конвертов, чтобы быть уверенным:
for (int i = 0; i < n; i += blocksize) {
for (int j = 0; j < n; j += blocksize) {
// transpose the block beginning at [i,j]
for (int k = i; k < i + blocksize; ++k) {
for (int l = j; l < j + blocksize; ++l) {
dst[k + l*n] = src[l + k*n];
}
}
}
}
Важным дальнейшим пониманием является то, что на самом деле для этого есть алгоритм, не обращающий внимания на кеш (см. http://en.wikipedia.org/wiki/Cache-oblivious_algorithm, в котором в качестве примера используется именно эта проблема). Неформальное определение "не обращающего внимания на кэш" состоит в том, что вам не нужно экспериментировать с настройкой каких-либо параметров (в данном случае размера блоков), чтобы достичь хорошей / оптимальной производительности кэша. Решение в этом случае состоит в том, чтобы транспонировать путем рекурсивного деления матрицы пополам и перемещения половин в их правильное положение в месте назначения.
Каким бы ни был размер кэша, эта рекурсия использует его в своих интересах. Я ожидаю, что по сравнению с вашей стратегией возникнут некоторые дополнительные издержки на управление, заключающиеся в том, чтобы использовать эксперименты с производительностью, чтобы, по сути, перейти прямо к той точке рекурсии, в которой кэш действительно срабатывает, и не идти дальше. С другой стороны, ваши эксперименты с производительностью могут дать вам ответ, который работает на вашей машине, но не на машинах ваших клиентов.
У меня была точно такая же проблема вчера. Я закончил с этим решением:
void transpose(double *dst, const double *src, size_t n, size_t p) noexcept {
THROWS();
size_t block = 32;
for (size_t i = 0; i < n; i += block) {
for(size_t j = 0; j < p; ++j) {
for(size_t b = 0; b < block && i + b < n; ++b) {
dst[j*n + i + b] = src[(i + b)*p + j];
}
}
}
}
Это в 4 раза быстрее, чем очевидное решение на моей машине.
Это решение заботится о прямоугольной матрице с размерами, не кратными размеру блока.
если dst и src - одна и та же квадратная матрица, вместо нее действительно должна использоваться функция на месте:
void transpose(double*m,size_t n)noexcept{
size_t block=0,size=8;
for(block=0;block+size-1<n;block+=size){
for(size_t i=block;i<block+size;++i){
for(size_t j=i+1;j<block+size;++j){
std::swap(m[i*n+j],m[j*n+i]);}}
for(size_t i=block+size;i<n;++i){
for(size_t j=block;j<block+size;++j){
std::swap(m[i*n+j],m[j*n+i]);}}}
for(size_t i=block;i<n;++i){
for(size_t j=i+1;j<n;++j){
std::swap(m[i*n+j],m[j*n+i]);}}}
Я использовал C++11, но это можно легко перевести на другие языки.
Вместо того, чтобы переставлять матрицу в памяти, почему бы не свернуть операцию транспонирования в следующую операцию, которую вы собираетесь выполнить над матрицей?
Стив Джессоп упомянул забытый в кэше алгоритм транспонирования матрицы. Для записи, я хочу поделиться возможной реализацией транспонирования матрицы кеша.
public class Matrix {
protected double data[];
protected int rows, columns;
public Matrix(int rows, int columns) {
this.rows = rows;
this.columns = columns;
this.data = new double[rows * columns];
}
public Matrix transpose() {
Matrix C = new Matrix(columns, rows);
cachetranspose(0, rows, 0, columns, C);
return C;
}
public void cachetranspose(int rb, int re, int cb, int ce, Matrix T) {
int r = re - rb, c = ce - cb;
if (r <= 16 && c <= 16) {
for (int i = rb; i < re; i++) {
for (int j = cb; j < ce; j++) {
T.data[j * rows + i] = data[i * columns + j];
}
}
} else if (r >= c) {
cachetranspose(rb, rb + (r / 2), cb, ce, T);
cachetranspose(rb + (r / 2), re, cb, ce, T);
} else {
cachetranspose(rb, re, cb, cb + (c / 2), T);
cachetranspose(rb, re, cb + (c / 2), ce, T);
}
}
}
Более подробную информацию о алгоритмах кеширования можно найти здесь.
Умножение матриц приходит на ум, но проблема с кешем там гораздо более выражена, потому что каждый элемент читается N раз.
С помощью транспонирования матрицы вы читаете за один линейный проход, и нет способа оптимизировать это. Но вы можете одновременно обрабатывать несколько строк, чтобы написать несколько столбцов и заполнить все строки кэша. Вам понадобится всего три петли.
Или сделайте это наоборот и читайте в колонках при линейной записи.
Вот пример кода C++ для транспонирования распараллеленной матрицы, удобной для кэширования:
constexpr auto range(unsigned max)
{ return std::views::iota(0u, max); }
constexpr auto range(unsigned start, unsigned max)
{ return std::views::iota(start, max); }
constexpr auto SZ = 1000u;
std::array<std::array<int, SZ>, SZ> mx;
void testtranspose(void)
{
for (int i = 0; i < SZ; i++)
for (int j = 0; j < SZ; j++)
mx[i][j] = j;
// ***************
// cache friendly in place parallel transposition
constexpr unsigned block_size = 10; assert(SZ % block_size == 0);
std::for_each(std::execution::par_unseq,
range(SZ / block_size).begin(), range(SZ / block_size).end(),
[&](const unsigned iib)
{ //cache friendly transposition
const unsigned ii = iib * block_size;
//transpose the diagonal block
for (unsigned i = ii; i < ii + block_size; i++)
for (unsigned j = i + 1; j < ii + block_size; j++)
std::swap(mx[i][j], mx[j][i]);
//transpose the rest of blocks
for (unsigned jj = ii + block_size; jj < SZ; jj += block_size)
for (unsigned i = ii; i < ii + block_size; i++)
for (unsigned j = jj; j < jj + block_size; j++)
std::swap(mx[i][j], mx[j][i]);
});
// ***************
for (int i = 0; i < SZ; i++)
for (int j = 0; j < SZ; j++)
assert(mx[i][j] == i);
}
С большой матрицей, возможно, большой разреженной матрицей, может быть идея разложить ее на более мелкие кеш-блоки (скажем, 4x4 подматрицы). Вы также можете пометить подматрицы как идентификаторы, которые помогут вам в создании оптимизированных путей кода.