Кеш и локальность: транспонирование блоками и кеш-независимые алгоритмы

Если ты знаешь, почему обход по столбцам медленный, как выбрать размер блока для транспонирования и что такое кеш-независимый алгоритм, листай до «Теперь сам».

Как думаешь, какой блок лучше

Матрица 16 × 16 транспонируется в кеше из восьми строк по четыре числа. Обычный двойной цикл даёт 320 промахов при нижней границе 128. Если транспонировать квадратами, какой размер квадрата выбрать: 2, 4 или 8?

Ответь себе до чтения. «Чем больше квадрат, тем лучше» - частый ответ. Квадрат 8 × 8 занимает по 16 строк кеша в каждой из двух матриц (строка квадрата - две строки кеша по четыре числа), а кеш - восемь строк на всё; он не помещается, и промахов снова 320. Лучше всех - 4 × 4.

На что это похоже

Переписать телефонную книгу из столбцов в строки. Если бегать от исходной страницы к новой ради каждого номера, большую часть времени займёт листание. Удобнее положить рядом две открытые страницы - исходную и новую - и переписать всё, что видно на развороте. Но если «разворот» слишком большой и не помещается на столе, страницы снова придётся перекладывать.

Блок транспонирования и есть разворот: обе его части должны поместиться в кеш одновременно.

Механизм: строки кеша, обход и блоки

Двумерный массив в C и numpy лежит по строкам: элемент [i][j] - по адресу i · n + j (в Java каждая строка - отдельный массив, но внутри строки элементы тоже идут подряд). Строка кеша (обычно 64 байта, 8 чисел double) содержит B соседних элементов одной строки матрицы.

• по строкам: соседние обращения - в одну строку кеша; n² / B промахов, минимум;
• по столбцам: каждое обращение - другая строка матрицы, другая строка кеша; если n строк кеша не помещаются, промах на каждое обращение;
• транспонирование T[j][i] = A[i][j] читает A по строкам и пишет T по столбцам - одна из двух матриц всегда «по столбцам».

Блоками: обрабатываем квадрат t × t: из A он занимает t строк кеша (если t ≤ B), из T - тоже t. Если все 2t · ⌈t / B⌉ строк помещаются в кеш (при t ≤ B это 2t ≤ M / B), каждая загруженная строка используется до конца, и промахов около 2n² / B - нижняя граница. Лучший t - чтобы строка блока была целой строкой кеша (t кратно B) и блок помещался.

Кеш-независимое транспонирование (Фриго, Лейзерсон, Прокоп, Рамачандран, 1999): делим длинную сторону матрицы пополам и рекурсивно транспонируем половины. Где-то в глубине рекурсии куски становятся такими, что помещаются в кеш, - каким бы он ни был. Параметры кеша в программе не участвуют, а промахов - столько же, сколько у лучшего блока (при обычном предположении «высокого» кеша, M ≥ B²). И сразу для всех уровней иерархии: L1, L2, L3, памяти и диска.

Правило

Массив лежит по строкам: обход по строкам даёт n² / B промахов, по столбцам - до n². Транспонирование делают блоками t × t, где все 2t · ⌈t/B⌉ строк кеша помещаются в кеш, а строка блока - целая строка кеша; тогда промахов около 2n² / B. Кеш-независимый алгоритм рекурсивно делит задачу пополам и получает тот же результат, не зная параметров кеша.

Разбор: сумма по строкам и по столбцам

Матрица 32 × 32, строка кеша - 8 чисел, кеш на 8 строк. Одна и та же сумма всех элементов, два порядка обхода.

def promahi(adresa, B, M):
    # кеш из M // B строк по B элементов, вытесняется давно не нужная
    kesh = []
    promahov = 0
    for a in adresa:
        blok = a // B
        if blok in kesh:
            kesh.remove(blok)
        else:
            promahov = promahov + 1
            if len(kesh) == M // B:
                kesh.pop(0)
        kesh.append(blok)
    return promahov

n = 32
B = 8
M = 64
po_strokam = []
po_stolbcam = []
for i in range(n):
    for j in range(n):
        po_strokam.append(i * n + j)
        po_stolbcam.append(j * n + i)
print("матрица", n, "x", n, "по строкам в памяти, строка кеша", B,
      "кеш на", M // B, "строк")
print("сумма по строкам:", promahi(po_strokam, B, M))
print("сумма по столбцам:", promahi(po_stolbcam, B, M))

матрица 32 x 32 по строкам в памяти, строка кеша 8 кеш на 8 строк
сумма по строкам: 128
сумма по столбцам: 1024

По строкам - 128 промахов, ровно 1024 / 8. По столбцам - 1024, промах на каждое обращение: пока программа спускается по столбцу через 32 строки матрицы, восемь строк кеша сменяются четыре раза, и к следующему столбцу ничего полезного не остаётся. Восьмикратная разница - ровно B.

Разбор: какой блок выбрать

Матрица 16 × 16, строка кеша - 4 числа, кеш на 8 строк. Транспонируем обычно, квадратами 2, 4 и 8 и рекурсивно.

def promahi(adresa, B, M):
    # кеш из M // B строк по B элементов, вытесняется давно не нужная
    kesh = []
    promahov = 0
    for a in adresa:
        blok = a // B
        if blok in kesh:
            kesh.remove(blok)
        else:
            promahov = promahov + 1
            if len(kesh) == M // B:
                kesh.pop(0)
        kesh.append(blok)
    return promahov

def obychno(n):
    # A[i][j] -> T[j][i]; A по адресам 0..n*n-1, T - следом за ней
    adresa = []
    for i in range(n):
        for j in range(n):
            adresa.append(i * n + j)
            adresa.append(n * n + j * n + i)
    return adresa

def blokami(n, t):
    adresa = []
    for bi in range(0, n, t):
        for bj in range(0, n, t):
            for i in range(bi, bi + t):
                for j in range(bj, bj + t):
                    adresa.append(i * n + j)
                    adresa.append(n * n + j * n + i)
    return adresa

def rekursivno(n, i0, j0, h, w, adresa):
    # делим длинную сторону пополам, пока кусок не станет маленьким;
    # размер кеша нигде не участвует
    if h * w <= 4:
        for i in range(i0, i0 + h):
            for j in range(j0, j0 + w):
                adresa.append(i * n + j)
                adresa.append(n * n + j * n + i)
    elif h >= w:
        rekursivno(n, i0, j0, h // 2, w, adresa)
        rekursivno(n, i0 + h // 2, j0, h - h // 2, w, adresa)
    else:
        rekursivno(n, i0, j0, h, w // 2, adresa)
        rekursivno(n, i0, j0 + w // 2, h, w - w // 2, adresa)

n = 16
B = 4
M = 32
print("n =", n, "B =", B, "строк кеша", M // B,
      "- нижняя граница 2n² / B =", 2 * n * n // B)
print("обычно:", promahi(obychno(n), B, M))
for t in [2, 4, 8]:
    print("блоками", t, "x", t, ":", promahi(blokami(n, t), B, M))
r = []
rekursivno(n, 0, 0, n, n, r)
print("рекурсивно:", promahi(r, B, M))

n = 16 B = 4 строк кеша 8 - нижняя граница 2n² / B = 128
обычно: 320
блоками 2 x 2 : 192
блоками 4 x 4 : 128
блоками 8 x 8 : 320
рекурсивно: 128

Обычно - 320. Квадраты 2 × 2 - 192: блок помещается, но использует половину каждой строки кеша, и строка загружается дважды. Квадраты 4 × 4 - 128, ровно нижняя граница: 4 строки из A и 4 из T - восемь, весь кеш, и каждая строка блока - целая строка кеша. Квадраты 8 × 8 - снова 320: 32 строки не помещаются в 8, блок вытесняет сам себя. Рекурсия - 128, как лучший блок, хотя про B и M она не знает ничего.

Что ломается без этого

Циклы не в том порядке. for j: for i: s += a[i][j] в C или numpy по массиву, лежащему по строкам, - в B раз больше промахов. Обмен двух строк for - часто самая дешёвая оптимизация в программе.

Блок, подобранный на одной машине. Размер, идеальный для кеша в 32 КБ, на машине с другим кешем или при другом размере элемента теряет выигрыш. Кеш-независимые алгоритмы от этого свободны.

Матрица со стороной, кратной большой степени двойки. Строки матрицы 4096 × 4096 попадают в одни и те же наборы ассоциативного кеша и вытесняют друг друга, хотя кеш почти пуст. Лечат дополнением строки несколькими лишними элементами.

Зачем это в работе

Библиотеки линейной алгебры (OpenBLAS, MKL, BLIS) - это в первую очередь подбор блоков под каждый уровень кеша: умножение матриц блоками в десятки раз быстрее наивного тройного цикла при той же арифметике. numpy.ascontiguousarray и порядок C/F в numpy существуют ради того же. Сортировки и поиск в стандартных библиотеках выбирают структуры, удобные кешу, - об этом раздел «Структуры данных на практике».

Теперь сам

1. Строка кеша - 64 байта, числа по 8 байт. Сколько промахов даст сумма матрицы 1000 × 1000 по строкам?
2. Кеш - 32 строки по 16 чисел. Каким может быть наибольший блок t × t при t, кратном 16, для транспонирования?
3. Почему рекурсивному транспонированию не нужен параметр B, но ему всё же нужна какая-то остановка рекурсии?