Что ученые должны знать о железе для написания быстрого кода

источник изображения

Программирование сегодня используется во многих областях науки, где отдельным ученым часто приходится собственноручно писать код для своих проектов. Для большинства ученых, однако, компьютерные науки не являются их областью знаний; они изучили программирование по необходимости. Я считаю себя одним из них. Хотя мы можем быть достаточно хорошо знакомы с программированием со стороны софта, мы редко имеем даже базовое представление о том, как железо влияет на производительность кода.

Цель этого урока — дать непрофессиональным программистам краткий обзор особенностей современного оборудования, которые нужно понимать, чтобы писать быстрый код. Это будет дистилляция того, что мы узнали за последние несколько лет. Этот учебник будет использовать Julia, потому что она позволяет легко продемонстрировать эти относительно низкоуровневые соображения на высокоуровневом интерактивном языке.

Это не руководство по языку программирования Julia

Чтобы писать быстрый код, вы должны сначала понять свой язык программирования и его особенности. Но это не руководство по языку программирования Julia. Я рекомендую прочитать раздел советы по производительности из документации Julia.

Это не объяснение конкретных структур данных или алгоритмов

Помимо знания вашего языка, вы также должны понимать свой собственный код, чтобы сделать его быстрым. Вы должны понять идею, лежащую в основе нотации big-O, почему некоторые алгоритмы работают быстрее других, и как различные структуры данных работают внутри. Не зная, что такое "массив", как можно оптимизировать код, использующий массивы?

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

• Как двоичное целое число представляется в памяти
• Как число с плавающей запятой представлено в памяти (изучение этого также необходимо для понимания вычислительных ошибок от операций с плавающей запятой, что является обязательным при выполнении научного программирования)
• Расположение в памяти "строк", включая кодировку ASCII и UTF-8
• Основы того, как структурирован "массив", и в чем разница между плотным массивом, например, целых чисел и массивом ссылок на объекты
• Принципы, лежащие в основе того, как работают Dict (т. е. хэш-таблица) и Set

Кроме того, я бы также рекомендовал ознакомиться с такими штуками как:

• Кучи
• Двусторонние очереди
• Кортежи

Это не учебник по бенчмаркингу вашего кода

Чтобы написать быстрый код на практике, необходимо его профилировать, чтобы найти узкие места, где ваша машина проводит большую часть времени. Необходимо сравнить различные функции и подходы, чтобы найти наиболее быстрые на практике. У Джулии (и других языков) есть инструменты именно для этой цели, но я не буду на них заостряться.

Содержание

• Минимизируйте записи на диск
• CPU cache
• Выравнивание
• Ассемблерный код
• Минимизируйте выделения памяти
• Используйте SIMD векторизации
• Структура массивов
• Используйте специализированные CPU инструкции
• Inline для маленьких функций
• Разворачивайте циклы
• Избегайте непредсказуемых ветвей
• Многопоточность
• GPUs

Прежде чем начать, следует установить пакеты:

# using Pkg
# Pkg.add("BenchmarkTools")
# Pkg.add("StaticArrays")

using StaticArrays
using BenchmarkTools

# Print median elapsed time of benchmark
function print_median(trial)
    println("Median time: ", BenchmarkTools.prettytime(median(trial).time))
end;

Базовая структура компьютерного оборудования

Начнем с упрощенной ментальной модели компьютера. Далее, будем добавлять детали к нашей модели по мере необходимости.

На этой простой диаграмме стрелки представляют поток данных в любом направлении. На диаграмме показаны три важные части компьютера:

• Центральный процессор (CPU) представляет собой чип размером со штамп. Именно здесь на самом деле происходят все вычисления — в мозгу компьютера.
• Оперативная память (Random access memory, ОЗУ, или просто "память") — это кратковременная память компьютера. Для поддержания этой памяти требуется электроэнергия, и она теряется при выключении компьютера. Оперативная память служит временным хранилищем данных между диском и процессором. Большая часть времени, затрачиваемого на "загрузку" различных приложений и операционных систем, фактически тратится на перемещение данных с диска в оперативную память и распаковку их там. Типичный потребительский ноутбук имеет около бит оперативной памяти.
• Диск является элементом хранения данных. Данные на диске сохраняются после отключения питания, поэтому диск содержит долговременную память компьютера. Он также намного дешевле в пересчете на гигабайт, чем оперативная память. Потребительский ПК имеет около бит дискового пространства.

Minimize disk writes

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

Поэтому для быстрой работы кода необходимо держать рабочие данные в оперативной памяти и максимально ограничить запись на диск.

Следующий пример служит для иллюстрации разницы в скорости: Первая функция открывает файл, получает доступ к одному байту из файла и снова закрывает его. Вторая функция случайным образом получает доступ к 1 000 000 целых чисел из оперативной памяти.

# Open a file
function test_file(path)
    open(path) do file
        # Go to 1000'th byte of file and read it
        seek(file, 1000)
        read(file, UInt8)
    end
end
@time test_file("README.md")

# Randomly access data N times
function random_access(data::Vector{UInt}, N::Integer)
    n = rand(UInt)
    mask = length(data) - 1
    @inbounds for i in 1:N
        n = (n >>> 7) ⊻ data[n & mask + 1]
    end
    return n
end

data = rand(UInt, 2^24)
@time random_access(data, 1000000);

# 0.000781 seconds (15 allocations: 992 bytes)
# 0.095022 seconds

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

На моем компьютере поиск одного байта в файле (включая открытие и закрытие файла) занимает около 781 мкс, а доступ к 1 000 000 целых чисел из памяти занимает 95 миллисекунд. Таким образом, оперативная память работает примерно в 8000 раз быстрее, чем диск.

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

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

CPU cache

Оперативная память быстрее диска, а процессор, в свою очередь, быстрее оперативной памяти. Процессор тикает, как часы, со скоростью около 3 ГГц, то есть 3 миллиарда тиков в секунду. Один "тик" этих часов называется тактовым циклом. Хотя это не совсем так, вы можете себе представить, что каждый цикл процессор выполняет одну простую команду, называемую инструкцией процессора, которая выполняет одну операцию над небольшим фрагментом данных. Тогда тактовая частота может служить ориентиром для других таймингов в компьютере. Стоит понимать, что за один такт фотон пройдет только около 10 см, и это ставит барьер тому, насколько быстро может работать память (которая расположена на некотором расстоянии от процессора). На самом деле, современные компьютеры настолько быстры, что существенным узким местом в их скорости является задержка, вызванная временем, необходимым для движения электричества по проводам внутри компьютера.

В этом масштабе чтение из оперативной памяти занимает около 100 тактов. Аналогично тому, как медлительность дисков может быть уменьшена путем копирования данных в более быструю оперативную память, данные из оперативной памяти копируются в меньший чип памяти непосредственно на процессоре, называемый кэш. Кэш быстрее, потому что он физически находится на чипе процессора (уменьшая путевые задержки), а также потому, что он использует более быстрый тип оперативной памяти, статическую оперативную память, вместо более медленной (но более дешевой в производстве) динамической оперативной памяти. Поскольку он должен быть размещен на процессоре, ограничивается его размер, и поскольку он более дорог в производстве, типичный кэш процессора содержит только около бит, примерно в 1000 раз меньше, чем оперативная память. На самом деле существует несколько уровней кэша процессора, но здесь мы упрощаем его и просто называем "кэш" как одну единственную вещь:

Когда процессор запрашивает часть данных из оперативной памяти, скажем, один байт, он сначала проверяет, находится ли память уже в кэше. Если да, то он будет читать оттуда. Это происходит намного быстрее чем доступ к оперативной памяти — обычно, всего за несколько тактов. Если нет, то получается промах кэша, и ваша программа будет тормозиться в течение десятков наносекунд, пока ваш компьютер копирует данные из оперативной памяти в кэш.

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

Эффективное использование кэша сводится к локальности, временной и пространственной:

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

Из этой информации можно вывести ряд простых трюков для повышения производительности:

• Используйте как можно меньше памяти. Когда ваши данные занимают меньше памяти, более вероятно, что ваши данные будут находиться в кэше. Помните, что процессор может выполнить около 100 небольших операций за время, потраченное впустую на один промах кэша.

  
  

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

  
  

Чтобы проиллюстрировать это, давайте сравним время, затраченное на функцию random_access приведенную выше, с новой функцией linear_access. Эта функция выполняет те же вычисления, что и random_access, но использует i вместо n для доступа к массиву, поэтому она обращается к данным линейным способом. Следовательно, единственное различие — это шаблон доступа к памяти.

Обратите внимание на большое расхождение во времени.

function linear_access(data::Vector{UInt}, N::Integer)
    n = rand(UInt)
    mask = length(data) - 1
    for i in 1:N
        n = (n >>> 7) ⊻ data[i & mask + 1]
    end
    return n
end

print_median(@benchmark random_access(data, 4096))
print_median(@benchmark linear_access(data, 4096))

# Median time: 408.149 μs
# Median time: 1.864 μs

Это также имеет последствия для ваших структур данных. Хэш-таблицы, такие как Dict и Set, по своей сути неэффективны в кэше и почти всегда вызывают промахи кэша, в то время как массивы — нет.

Многие из оптимизаций в этом документе косвенно влияют на использование кэша, поэтому это важно иметь в виду.

Alignment

Как уже упоминалось, ваш процессор будет перемещать 512 последовательных бит (64 байта) в основную оперативную память и из нее в кэш одновременно. Эти 64 байта называются строкой кэша. Вся ваша основная память сегментирована на строки кэша. Например, адреса памяти от 0 до 63 — это одна строка кэша, адреса от 64 до 127-следующая, от 128 до 191-следующая и так далее. Ваш процессор может запрашивать только одну из этих строк из памяти, а не, например, 64 байта с адреса 30 до 93.

Это означает, что некоторые структуры данных могут пересекать границы между строками кэша. Если я запрошу 64-битное (8 байтовое) целое число по адресу 60, процессор должен сгенерировать два запроса памяти (а именно, чтобы получить строки кэша 0-63 и 64-127), а затем получить целое число из обеих строк кэша, теряя время.

Время, потраченное впустую, может быть значительным. В ситуации, когда доступ к кэш-памяти оказывается узким местом, замедление может приблизиться к двукратному. В следующем примере я использую указатель для повторного доступа к массиву с заданным смещением от границы строки кэша. Если смещение находится в диапазоне 0:56, то все целые числа помещаются в одну строку кэша, и функция работает быстро. Если смещение находится в 57:63, то все целые числа будут располагаться между строками кэша.

function alignment_test(data::Vector{UInt}, offset::Integer)
    # Jump randomly around the memory.
    n = rand(UInt)
    mask = (length(data) - 9) ⊻ 7
    GC.@preserve data begin # protect the array from moving in memory
        ptr = pointer(data)
        iszero(UInt(ptr) & 63) || error("Array not aligned")
        ptr += (offset & 63)
        for i in 1:4096
            n = (n >>> 7) ⊻ unsafe_load(ptr, (n & mask + 1) % Int)
        end
    end
    return n
end
data = rand(UInt, 256 + 8);

print_median(@benchmark alignment_test(data, 0))
print_median(@benchmark alignment_test(data, 60))

# Median time: 6.707 μs
# Median time: 12.680 μs

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

data = rand(UInt, 1 << 24 + 8)
print_median(@benchmark alignment_test(data, 10))
print_median(@benchmark alignment_test(data, 60))

# Median time: 404.351 μs
# Median time: 484.438 μs

К счастью, компилятор делает несколько трюков, чтобы уменьшить вероятность того, что вы получите доступ к несогласованным данным. Во-первых, Джулия (и другие компилируемые языки) всегда помещает новые объекты в память на границах строк кэша. Когда объект помещается прямо на границе, мы говорим, что он выровнен. Джулия также выравнивает начало больших массивов:

memory_address = reinterpret(UInt, pointer(data))
@assert iszero(memory_address % 64)

Обратите внимание, что если начало массива выровнено, то для 1-, 2-, 4- или 8-байтовых объектов невозможно пересечь границы строк кэша, и все будет выровнено.

По-прежнему было бы возможно, например, чтобы 7-байтовый объект был смещен в массиве. В массиве 7-байтовых объектов 10-й объект будет помещен со смещением байта , и объект перешагнет строку кэша. Однако компилятор обычно не допускает структуру с нестандартным размером. Определим 7-байтовую структуру:

struct AlignmentTest
    a::UInt32 # 4 bytes +
    b::UInt16 # 2 bytes +
    c::UInt8  # 1 byte = 7 bytes?
end

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

function get_mem_layout(T)
    for fieldno in 1:fieldcount(T)
        println("Name: ", fieldname(T, fieldno), "\t",
                "Size: ", sizeof(fieldtype(T, fieldno)), " bytes\t",
                "Offset: ", fieldoffset(T, fieldno), " bytes.")
    end
end

println("Size of AlignmentTest: ", sizeof(AlignmentTest), " bytes.")
get_mem_layout(AlignmentTest)

# Size of AlignmentTest: 8 bytes.
# Name: a   Size: 4 bytes   Offset: 0 bytes.
# Name: b   Size: 2 bytes   Offset: 4 bytes.
# Name: c   Size: 1 bytes   Offset: 6 bytes.

Можно заметить, что, несмотря на то, что AlignmentTest имеет только 4 + 2 + 1 = 7 байт фактических данных, он занимает 8 байт памяти, и доступ к объекту AlignmentTest из массива всегда будет выровнен.

Есть только несколько ситуаций, в которых вы можете столкнуться с проблемами выравнивания, как программист. Я могу придумать два варианта:

1. Если вы вручную создаете объект со странным размером, например, обращаясь к плотному целочисленному массиву с указателями. Это может сэкономить память, но будет времязатратно. Моя реализация фильтра кукушки делает так для экономии места.
2. Во время матричных операций. Если у вас есть матрица, то столбцы иногда не выровнены, потому что она плотно хранится в памяти. например, в матрице 15x15 Float32 выровнен только первый столбец, а все остальные — нет. Это может иметь серьезные последствия при выполнении матричных операций: Я видел бенчмарки где умножение матрицы/вектора 80х80 происходит в 2 раза быстрее, чем умножение матрицы/вектора 79х79 из-за проблем с выравниванием.

Inspect generated assembly

Для запуска любая программа должна быть переведена или скомпилирована в инструкции процессора. Инструкции процессора — это то, что на самом деле работает на вашем компьютере, в отличие от кода, написанного на вашем языке программирования, который является просто описанием программы. Инструкции процессора обычно представляются людям в ассемблер-коде. Ассемблер — это язык программирования, который имеет однозначное соответствие с инструкциями процессора.

Просмотр низкоуровнего кода будет полезен для понимания некоторых из следующих разделов, которые относятся к инструкциям процессора.

В Julia мы можем легко проверить скомпилированный ассемблерный код с помощью функции code_native или эквивалентного макроса @code_native. Проделаем это для простой функции:

# View assembly code generated from this function call
function foo(x)
    s = zero(eltype(x))
    @inbounds for i in eachindex(x)
        s = x[i ⊻ s]
    end
    return s
end

# Actually running the function will immediately crash Julia, so don't.
@code_native foo(data)

section __TEXT,__text,regular,pure_instructions
; ┌ @ In[13]:4 within 'foo'
; │┌ @ abstractarray.jl:212 within 'eachindex'
; ││┌ @ abstractarray.jl:95 within 'axes1'
; │││┌ @ abstractarray.jl:75 within 'axes'
; ││││┌ @ In[13]:3 within 'size'
    movq    24(%rdi), %rax
; ││││└
; ││││┌ @ tuple.jl:157 within 'map'
; │││││┌ @ range.jl:320 within 'OneTo' @ range.jl:311
; ││││││┌ @ promotion.jl:409 within 'max'
    testq   %rax, %rax
; │└└└└└└
    jle L75
    movq    %rax, %rcx
    sarq    $inlinе$63, %rcx
    andnq   %rax, %rcx, %rcx
    movq    (%rdi), %rdx
; │ @ In[13]:5 within 'foo'
; │┌ @ int.jl:860 within 'xor' @ int.jl:301
    negq    %rcx
    movl    $inlinе$1, %esi
    xorl    %eax, %eax
    nopw    %cs:(%rax,%rax)
    nopl    (%rax)
L48:
    xorq    %rsi, %rax
; │└
; │┌ @ multidimensional.jl:543 within 'getindex' @ array.jl:788
    movq    -8(%rdx,%rax,8), %rax
; │└
; │┌ @ range.jl:597 within 'iterate'
; ││┌ @ promotion.jl:398 within '=='
    leaq    (%rcx,%rsi), %rdi
    addq    $inlinе$1, %rdi
; ││└
    addq    $inlinе$1, %rsi
; ││┌ @ promotion.jl:398 within '=='
    cmpq    $inlinе$1, %rdi
; │└└
    jne L48
; │ @ In[13]:7 within 'foo'
    retq
L75:
    xorl    %eax, %eax
; │ @ In[13]:7 within 'foo'
    retq
    nop
; └

Разберемся с этим:

Строки, начинающиеся с ;, являются комментариями и объясняют, из какого раздела кода берутся следующие инструкции. Они показывают вложенные серии вызовов функций и то, где в исходном коде они находятся. Заметим, что eachindex, вызывает axes1, которая обращается к axes, которая вызывает size. Под строкой комментария, содержащей вызов size, мы видим первую инструкцию процессора. Название инструкции находится в крайнем левом углу, movq. Название состоит из двух частей: mov, типа инструкции (для перемещения содержимого в регистр или из него), и суффикса q, сокращенного от "quad", что означает 64-битное целое число. Существуют следующие суффиксы: b (байт, 8 бит), w (слово, 16 бит), l (длинный, 32 бит) и q (quad, 64 бит).

Следующие два столбца инструкции, 24(%rdi) и %rax, являются аргументами для movq. Это имена регистров (мы вернемся к ним позже), в которых хранятся данные для работы.

Вы также можете заметить (на увеличенном дисплее ассемблер-кода), что код сегментирован на разделы, начинающиеся с имени, с буквой "L", например, есть раздел L48. Переход между этими разделами происходит при использовании операторов if или ветвей. Здесь раздел L48 отмечает фактический цикл. Рассмотрим следующие две инструкции в разделе L48:

; ││┌ @ promotion.jl:401 within '=='
     cmpq    $inlinе$1, %rdi
; │└└
     jne     L48

Первая инструкция cmpq (compare quad) сравнивает данные в реестре rdi, которые содержат данные о количестве оставшихся итераций (плюс одна), с числом 1 и устанавливает определенные флаги (провода) в процессоре на основе результата. Следующая команда jne (jump if not equal) совершает прыжок, если флаг "equal" не установлен в процессоре, что происходит, если остается одна или несколько итераций. Он переходит к "L48", что означает повторение этого раздела.

Быстрая инструкция, медленная инструкция

Не все инструкции процессора одинаково быстры. Ниже приведена таблица избранных инструкций процессора с очень приблизительными оценками того, сколько тактов им требуется для выполнения. Вы можете найти гораздо более подробные таблицы в этом документе. Здесь я кратко опишу скорость выполнения инструкций на современных процессорах Intel. Это очень похоже на все современные процессоры.

Инструкции процессоров обычно занимают несколько циклов процессора для завершения. Однако, если инструкция использует другую часть ЦП во время ее выполнения, ЦП обычно может запустить новую инструкцию до завершения старой: если какая-то операция X занимает, скажем, 4 такта, они могут поставить в очередь одну или даже две операции за такт, используя функцию, называемую конвейеризация инструкций. Следовательно, Инструкция X имеет задержку в 4 цикла, что означает, что для завершения инструкции требуется 4 цикла. Но если процессор может поставить в очередь новую инструкцию в каждый отдельный цикл, он может иметь взаимную пропускную способность 1 такт, что означает в среднем, что требуется только 1 цикл на операцию.

В следующей таблице время измеряется в тактах:

Команда lea принимает три входа, A, B и C, где A должно быть 2, 4 или 8, и вычисляет AB + C. Мы вернемся к тому, что делают инструкции "вектор" позже.

Для сравнения мы можем также добавить некоторые очень грубые оценки других источников задержек:

Если у вас есть внутренний цикл, выполняющийся миллионы раз, вполне резонно проверить сгенерированный низкоуровневый код для цикла и проверить, можете ли вы выразить вычисления в терминах быстрых инструкций процессора. Например, если у вас есть целое число, которое, как вы знаете, равно 0 или выше, и вы хотите разделить его на 8 (отбрасывая любой остаток), вы можете вместо этого сделать битовый сдвиг, так как битовые сдвиги намного быстрее, чем целочисленное деление:

divide_slow(x) = div(x, 8)
divide_fast(x) = x >>> 3;

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

# Calling it with debuginfo=:none removes the comments in the assembly code
code_native(divide_slow, (UInt,), debuginfo=:none)

        .section    __TEXT,__text,regular,pure_instructions
    movq    %rdi, %rax
    shrq    $inlinе$3, %rax
    retq
    nopl    (%rax,%rax)

Minimize allocations

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

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

В таких языках программирования, как Julia, Python, R и Java, освобождение производится автоматически с помощью программы, называемой сборщиком мусора (GC). Эта программа отслеживает, какие объекты оказываются недоступными программисту, и освобождает их. Например:

thing = [1,2,3]
thing = nothing

Нет никакого способа вернуть исходный массив "[1,2,3] " обратно, он недостижим. Поэтому он просто тратит впустую оперативную память и ничего не делает. Это "мусор". Выделение и освобождение объектов иногда приводит к тому, что GC начинает сканирование всех объектов в памяти и освобождает недостижимые, что вызывает значительное отставание. Вы также можете запустить сборщик мусора вручную:

GC.gc()

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

function increment(x::Vector{<:Integer})
    y = similar(x)
    @inbounds for i in eachindex(x)
        y[i] = x[i] + 1
    end
    return y
end

function increment!(x::Vector{<:Integer})
    @inbounds for i in eachindex(x)
        x[i] = x[i] + 1
    end
    return x
end

function print_mean(trial)
    println("Mean time: ", BenchmarkTools.prettytime(BenchmarkTools.mean(trial).time))
end

data = rand(UInt, 2^10);

# Запустите один раз, чтобы скомпилировать функцию - мы не хотим измерять компиляцию
increment(data); increment!(data)

print_mean(@benchmark increment(data));
print_mean(@benchmark increment!(data));

# Mean time: 1.799 μs
# Mean time: 81.788 ns

На моем компьютере функция выделяющая память в среднем работает более чем в 20 раз медленнее. Это связано с несколькими свойствами кода:

• Во-первых, само выделение требует времени
• Во-вторых, выделенные объекты в конечном итоге должны быть освобождены, что также требует времени
• В-третьих, повторные выделения запускают GC, вызывая накладные расходы
• В-четвертых, большее количество выделений иногда означает менее эффективное использование кэша, поскольку вы используете больше памяти

Обратите внимание, что я использовал среднее время вместо медианы, так как для этой функции GC запускает только приблизительно каждый 30-й вызов, но при этом потребляет 30-40 мкс. Все это означает, что производительный код должен сводить выделение ресурсов к минимуму. Обратите внимание, что макрос @btime печатает количество и размер выделений памяти. Эта информация дается потому, что предполагается, что любой программист, который заботится о тестировании своего кода, будет заинтересован в сокращении аллокаций.

Не на все объекты нужно выделять память

Внутри оперативной памяти данные хранятся либо в стеке, либо в куче. Стек — это простая структура данных с началом и концом, похожая на Vector в Julia. Стек может быть изменен только путем добавления или вычитания элементов из конца, аналогично Vector'у только с двумя мутирующими операциями push! и pop!. Эти операции в стеке выполняются очень быстро. Однако, когда мы говорим об аллокациях, мы говорим о данных в куче. В отличие от стека, куча имеет неограниченный размер (ну, она имеет размер оперативной памяти вашего компьютера) и может быть изменена произвольно, при удалении любых объектов.

Интуитивно может показаться очевидным, что все объекты будучи помещены в оперативную память, должны иметь возможность быть удаленными в любое время программой, и поэтому должны быть выделены в куче. И для некоторых языков, таких как Python, это верно. Однако это не относится к Джулии и другим эффективным компилируемым языкам. Целые числа, например, часто могут быть помещены в стек.

Почему некоторые объекты должны выделяться как куча, а другие могут быть в виде стека? Чтобы быть заточенным под стек, компилятор должен точно знать, что:

• Объект имеет достаточно маленький размер и помещается в стек. Это необходимо по техническим причинам для работы стека.
• Компилятор может точно предсказать, когда ему нужно добавить и уничтожить объект, просто открыв стек (аналогично вызову pop! на Vector). Обычно это относится к локальным переменным в компилируемых языках.

Джулия имеет еще больше ограничений на объекты, выделенные стеком.

• Объект должен иметь фиксированный размер, известный во время компиляции.
• Компилятор должен знать, что объект никогда не изменяется. Процессор может свободно копировать объекты, выделенные стеком, а для неизменяемых объектов нет способа отличить копию от оригинала. Это стоит повторить: с неизменяемыми объектами невозможно отличить копию от оригинала. Это дает компилятору и процессору определенные свободы при работе с ним.

В Julia у нас есть понятие bitstype, которое представляет собой объект, который рекурсивно содержит объекты, выделяемые не кучей. Выделенные в куче объекты — это объекты типов String, Array, Symbol, изменяемые объекты или объекты, содержащие любой из предыдущих. Битовые типы более эффективны именно потому, что они неизменяемы, фиксированы по размеру и почти всегда могут быть оформлены как стек. Последний пункт также объясняет, почему объекты неизменяемы по умолчанию в Julia, и приводит к еще одному совету производительности: используйте неизменяемые объекты везде, где это возможно.

Что это означает на практике? В Julia это означает, что если вы хотите быстрые выделяемые как стек объекты, то:

• Ваш объект должен быть создан, использован и уничтожен в полностью скомпилированной функции, чтобы компилятор точно знал, когда ему нужно создать, использовать и уничтожить объект. Если объект возвращается для последующего использования (а не сразу возвращается в другую, полностью скомпилированную функцию), мы говорим, что объект экранируется и должен быть выделен.
• Тип вашего объекта должен быть битовым типом.
• Ваш тип должен быть ограничен по размеру. Я точно не знаю, насколько большим он должен быть, но 100 байт — это нормально.
• Точные требования к памяти вашего типа должно быть известно компилятору.

abstract type AllocatedInteger end

struct StackAllocated <: AllocatedInteger
    x::Int
end

mutable struct HeapAllocated <: AllocatedInteger
    x::Int
end

Сравним подноготную приведенных выше объектов:

@code_native HeapAllocated(1)

    .section    __TEXT,__text,regular,pure_instructions
; ┌ @ In[20]:8 within 'HeapAllocated'
    pushq   %rbx
    movq    %rsi, %rbx
    movabsq $inlinе$jl_get_ptls_states_fast, %rax
    callq   *%rax
    movabsq $inlinе$jl_gc_pool_alloc, %rcx
    movq    %rax, %rdi
    movl    $inlinе$1400, %esi             ## imm = 0x578
    movl    $inlinе$16, %edx
    callq   *%rcx
    movabsq $inlinе$4470743600, %rcx       ## imm = 0x10A7A2230
    movq    %rcx, -8(%rax)
    movq    %rbx, (%rax)
    popq    %rbx
    retq
    nopl    (%rax)
; └

Обратите внимание на инструкции callq в HeapAllocated. Эта инструкция вызывает другие функции, а это означает, что на самом деле для создания объекта HeapAllocated, действительно требуется гораздо больше кода. В отличии от этого, в StackAllocated нужно всего несколько инструкций:

@code_native StackAllocated(1)

    .section    __TEXT,__text,regular,pure_instructions
; ┌ @ In[20]:4 within 'StackAllocated'
    movq    %rsi, %rax
    retq
    nopw    %cs:(%rax,%rax)
    nop
; └

Поскольку битовые типы не должны храниться в куче и могут быть скопированы свободно, то они хранятся inline в массивах. Это означает, что такого рода объекты могут храниться непосредственно в памяти массива. Не-битовые типы имеют уникальную идентичность и расположение в куче. Они отличаются от копий, поэтому не могут быть свободно скопированы, и поэтому массивы содержат ссылку на место памяти в куче, где они хранятся. Доступ к такому объекту из массива тогда означает сначала доступ к массиву, чтобы получить ячейку памяти, а затем доступ к самому объекту, используя эту ячейку памяти. Помимо двойного доступа к памяти, объекты хранятся в куче менее эффективно, а это означает, что больше памяти необходимо скопировать в кэш процессора, что означает больше пропусков кэша. Следовательно, даже при хранении в куче в массиве битовые типы могут храниться более эффективно.

Base.:+(x::Int, y::AllocatedInteger) = x + y.x
Base.:+(x::AllocatedInteger, y::AllocatedInteger) = x.x + y.x

data_stack = [StackAllocated(i) for i in rand(UInt16, 1000000)]
data_heap = [HeapAllocated(i.x) for i in data_stack]

@btime sum(data_stack)
@btime sum(data_heap);

#   130.958 μs (1 allocation: 16 bytes)
#   951.989 μs (1 allocation: 16 bytes)

Мы можем удостовериться, что массив data_stack действительно хранит данные StackAllocated, тогда как data_heap содержит указатели (т. е. адреса памяти):

println("First object of data_stack:         ", data_stack[1])
println("First data in data_stack array:     ", unsafe_load(pointer(data_stack)), '\n')

println("First object of data_heap:          ", data_heap[1])
first_data = unsafe_load(Ptr{UInt}(pointer(data_heap)))
println("First data in data_heap array:      ", repr(first_data))
println("Data at address ", repr(first_data), ": ",
        unsafe_load(Ptr{HeapAllocated}(first_data)))

#=
First object of data_stack:         StackAllocated(10099)
First data in data_stack array:     StackAllocated(10099)

First object of data_heap:          HeapAllocated(10099)
First data in data_heap array:      0x0000000139a7e570
Data at address 0x0000000139a7e570: HeapAllocated(10099)
=#

Exploit SIMD vectorization

Настало время еще раз обновить нашу упрощенную компьютерную схему. Процессор работает только с данными, содержащимися в регистрах. Это небольшие слоты фиксированного размера (например, размером 8 байт) внутри самого процессора. Регистр предназначен для хранения одного фрагмента данных, например целого числа или числа с плавающей запятой. Как указано в разделе, посвященном ассемблерному коду, каждая инструкция обычно относится к одному или двум регистрам, содержащим данные, с которыми работает операция:

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

@code_native UInt128(5) + UInt128(11)

    .section    __TEXT,__text,regular,pure_instructions
; ┌ @ int.jl:53 within '+'
    movq    %rdi, %rax
    addq    %rcx, %rsi
    adcq    %r8, %rdx
    movq    %rsi, (%rdi)
    movq    %rdx, 8(%rdi)
    retq
    nopw    %cs:(%rax,%rax)
    nopl    (%rax,%rax)
; └

Там нет регистра, который может сделать 128-битное сложение. Сначала младшие 64 бита должны быть добавлены с помощью инструкции addq, помещенной в регистр. Затем верхние биты добавляются с помощью инструкции adcxq, которая добавляет цифры, но также использует бит переноса из предыдущей инструкции. Наконец, результаты перемещаются по 64 бита за раз с помощью инструкций movq.

Малый размер регистров служит узким местом для пропускной способности процессора: он может работать только с одним int/float за раз. Чтобы обойти это, современные процессоры содержат специализированные 256-битные регистры (или 128-битные в старых процессорах, или 512-битные в совершенно новых), которые могут содержать 4 64-битных int/float одновременно, или 8 32-битных целых чисел и т. д. Смущает то, что данные в таких широких регистрах называются "векторами". ЦП имеет доступ к инструкциям, которые могут выполнять различные операции над векторами, оперируя 4 64-битными целыми числами в одной инструкции. Это называется "одна инструкция, несколько данных", SIMD или векторизация. Примечательно, что 4x64-битная операция не то же самое, что 256-битная операция, например, нет переноса между 4 64-битными целыми числами при добавлении двух векторов. Вместо этого 256-битная векторная операция эквивалентна 4 отдельным 64-битным операциям.

Мы можем проиллюстрировать это на следующем примере:

# Create a single statically-sized vector of 8 32-bit integers
# I could also have created 4 64-bit ones, etc.
a = @SVector Int32[1,2,3,4,5,6,7,8]

# Don't add comments to output
code_native(+, (typeof(a), typeof(a)), debuginfo=:none)

    .section    __TEXT,__text,regular,pure_instructions
    movq    %rdi, %rax
    vmovdqu (%rdx), %ymm0
    vpaddd  (%rsi), %ymm0, %ymm0
    vmovdqu %ymm0, (%rdi)
    vzeroupper
    retq
    nopw    %cs:(%rax,%rax)
    nopl    (%rax)

Здесь два 8*32-битных вектора суммируются вместе одной инструкцией. Здесь процессор использует одну инструкцию vpaddd (vector packed add double) для суммирования 8 32-битных целых чисел, а также соответствующую инструкцию перемещения vmovdqu. Обратите внимание, что векторные инструкции для процессора начинаются с буквы v.

Стоит упомянуть о взаимодействии между SIMD и выравниванием: если серия 256-битных (32-байтовых) нагрузок SIMD не выровнена, то до половины нагрузок могут пересекать границы линий кэша, в отличие от всего лишь 1/8 из 8-байтовых нагрузок. Таким образом, выравнивание является гораздо более серьезной проблемой при использовании SIMD. Поскольку начало массива всегда выровнено, это обычно не является проблемой, но в тех случаях, когда вы не гарантируете, что начнете с выровненной начальной точки, например при матричных операциях, это может иметь существенное значение. В совершенно новых процессорах с 512-битными регистрами проблемы еще хуже, так как размер SIMD совпадает с размером строки кэша, поэтому все нагрузки будут смещены, если будет смещена начальная нагрузка.

Векторизация SIMD, например, 64-битных целых чисел может увеличить пропускную способность почти в 4 раза, поэтому она имеет огромное значение в высокопроизводительном программировании. Компиляторы автоматически векторизуют операции, если это возможно. Что может помешать этой автоматической векторизации?

SIMD нуждается в непрерывной итерации фиксированной длины

Поскольку векторизованные операции работают с несколькими данными одновременно, невозможно прервать цикл в произвольной точке. Например, если 4 64-битных целых числа обрабатываются за один такт, то невозможно остановить цикл SIMD после обработки 3 целых чисел. Предположим, у вас есть такой цикл:

for i in 1:8
    if foo()
        break
    end
    # do stuff with my_vector[i]
end

Здесь цикл может закончиться на любой итерации из-за оператора break. Поэтому любая SIMD-инструкция, загруженная для нескольких целых чисел, может работать с данными после того, как цикл разрывается, то есть с данными, которые никогда не должны быть прочитаны. Это было бы неправильным поведением, и поэтому компилятор не может использовать инструкции SIMD.

Хорошее эмпирическое правило заключается в том, что для simd предпочтительны:

• Цикл с заданной длиной, поэтому она знает, когда остановиться, и
• Цикл без ветвей (т. е. операторов if) в цикле

Фактически, даже boundschecking, то есть проверка того, что вы не индексируете вне границ вектора, вызывает ветвление. В конце концов, если предполагается, что код вызывает ошибку границ после 3 итераций, даже одна операция SIMD будет неправильной! Чтобы достичь векторизации SIMD, все boundschecks должны быть отключены. Мы можем использовать это для демонстрации влияния SIMD:

function sum_nosimd(x::Vector)
    n = zero(eltype(x))
    for i in eachindex(x)
        n += x[i]
    end
    return n
end

function sum_simd(x::Vector)
    n = zero(eltype(x))
    # By removing the boundscheck, we allow automatic SIMD
    @inbounds for i in eachindex(x)
        n += x[i]
    end
    return n
end

# Make sure the vector is small enough to fit in cache so we don't time cache misses
data = rand(UInt64, 4096);

@btime sum_nosimd(data)
@btime sum_simd(data);

#  1.777 μs (1 allocation: 16 bytes)
#  182.045 ns (1 allocation: 16 bytes)

На моем компьютере код SIMD работает в 10 раз быстрее, чем код без SIMD. Только на SIMD приходится всего около 4-кратных улучшений (с тех пор как мы перешли от 64 бит на итерацию к 256 битам на итерацию). Остальная часть выигрыша происходит от того, что вы не тратите время на проверку границ, и от автоматического развертывания цикла (объясненного позже), что также стало возможным благодаря аннотации @inbounds.

SIMD нуждается в цикле, где порядок циклов не имеет значения

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

Представьте, что мы хотим суммировать некоторые 64-битные целые числа в массиве с помощью SIMD. Для простоты предположим, что массив состоит из 8 элементов: A, B, C… H. В обычном цикле, отличном от SIMD, сложения будут выполняться следующим образом:

Тогда как при загрузке целых чисел с помощью SIMD четыре 64-битных целых числа будут загружены в один вектор <A, B, C, D>, а остальные четыре — в другой <E, F, G, H>. Эти два вектора будут суммированы: <A+E, B+F, C+G, D+H>. После цикла будут суммированы четыре целых числа в результирующем векторе. Таким образом, общий порядок будет следующим:

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

x = eps(1.0) * 0.4
1.0 + (x + x) == (1.0 + x) + x
# false

по этой причине суммирование для float не будет автоматически векторизовано:

data = rand(Float64, 4096)
@btime sum_nosimd(data)
@btime sum_simd(data);

#   3.510 μs (1 allocation: 16 bytes)
#   3.836 μs (1 allocation: 16 bytes)

Однако высокопроизводительные языки программирования обычно предоставляют команду, сообщающую компилятору, что можно переупорядочить цикл, даже для неассоциативных циклов. В Julia этим занимается макрос @simd :

function sum_simd(x::Vector)
    n = zero(eltype(x))
    # Here we add the '@simd' macro to allow SIMD of floats
    @inbounds @simd for i in eachindex(x)
        n += x[i]
    end
    return n
end

data = rand(Float64, 4096)
@btime sum_nosimd(data)
@btime sum_simd(data);

#   3.510 μs (1 allocation: 16 bytes)
#   247.368 ns (1 allocation: 16 bytes)

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

Struct of arrays

Если мы создадим массив, содержащий четыре объекта AlignmentTest A, B, C и D, то объекты будут лежать в массиве впритык, вот так:

Objects: |      A        |       B       |       C       |        D      |
Fields:  |   a   | b |c| |   a   | b |c| |   a   | b |c| |   a   | b |c| |
Byte:     1               9              17              25              33

Еще раз обратите внимание, что байты № 8, 16, 24 и 32 пусты для сохранения выравнивания, что приводит к потере памяти.
Теперь предположим, что вы хотите выполнить операцию со всеми полями структур .a (broadcasting). Поскольку поля .a разбросаны на 8 байт друг от друга, операции SIMD гораздо менее эффективны (загрузка до 4 полей одновременно), чем если бы все поля .a хранились вместе (где 8 полей могли бы поместиться в 256-битном регистре). При работе только с полями .a будут считываться все 64-байтовые строки кэша, из которых только половина или 32 байта будут полезны. Мало того, что это приводит к большему количеству пропусков кэша, нам также потребуются инструкции, чтобы выделить половину нужных данных из регистров SIMD.

Структура памяти, которую мы имеем выше, называется "массив структур", потому что, ну, это массив, заполненный структурами. Вместо этого мы можем структурировать наши 4 объекта от A до D как "структуру массивов". Концептуально это может выглядеть так:

struct AlignmentTestVector
    a::Vector{UInt32}
    b::Vector{UInt16}
    c::Vector{UInt8}
end

Со следующим распределением памяти для каждого поля:

Object: AlignmentTestVector
.a |   A   |   B   |   C   |   D   |
.b | A | B | C | D |
.c |A|B|C|D|

Выравнивание больше не является проблемой, никакое пространство не тратится впустую на компоновку. При прохождении через поля а все строки кэша содержат полные 64 байта релевантных данных, поэтому операции SIMD не нуждаются в дополнительных операциях для выбора релевантных данных:

Base.rand(::Type{AlignmentTest}) = AlignmentTest(rand(UInt32), rand(UInt16), rand(UInt8))

N  = 1_000_000
array_of_structs = [rand(AlignmentTest) for i in 1:N]
struct_of_arrays = AlignmentTestVector(rand(UInt32, N), rand(UInt16, N), rand(UInt8, N));

@btime sum(x -> x.a, array_of_structs)
@btime sum(struct_of_arrays.a);

#   296.427 μs (1 allocation: 16 bytes)
#   86.857 μs (1 allocation: 16 bytes)

Use specialized CPU instructions

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

Джулия предлагает несколько специальных инструкций, таких как:

• Количество ненулевых битов в целочисленном числе эффективно подсчитывается с помощью инструкции popcnt, предоставляемой через функцию count_ones.
• Инструкции tzcnt подсчитывают количество конечных нулей в битах целых чисел, выставляемое с помощью функции trailing_zeros.
• Порядок отдельных байтов в многобайтовом целом может быть изменен с помощью инструкции bswap, предоставляемой через функцию bswap. Это может быть полезно, когда приходится иметь дело с endianness.

Следующий пример иллюстрирует разницу в производительности между ручной реализацией функции count_ones и встроенной версией, которая использует инструкцию popcnt:

function manual_count_ones(x)
    n = 0
    while x != 0
        n += x & 1
        x >>>= 1
    end
    return n
end

data = rand(UInt, 10000)
@btime sum(manual_count_ones, data)
@btime sum(count_ones, data);

#  218.341 μs (1 allocation: 16 bytes)
#  2.009 μs (1 allocation: 16 bytes)

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

Вызов инструкций процессора

Джулия позволяет напрямую вызывать инструкции процессора. Это обычно не рекомендуется, так как не все ваши пользователи будут иметь доступ к одному и тому же процессору с одинаковыми инструкциями.

Последние процессоры содержат специальные инструкции для AES-шифрования и SHA256-хэширования. Если вы хотите вызвать эти инструкции, вы можете использовать бэкэнд-компилятор Julia, LLVM, напрямую. В приведенном ниже примере я создаю функцию, которая непосредственно вызывает инструкцию vaesenc (один проход шифрования AES) :

# This is a 128-bit CPU "vector" in Julia
const __m128i = NTuple{2, VecElement{Int64}}

# Define the function in terms of LLVM instructions
aesenc(a, roundkey) = ccall("llvm.x86.aesni.aesenc", llvmcall, __m128i, (__m128i, __m128i), a, roundkey);

Мы можем проверить его работу, просмотрев сборку функции, которая должна содержать только одну инструкцию vaesenc, а также инструкцию retq (return) и nopw (do nothing, используемый в качестве наполнителя для выравнивания инструкций процессора в памяти) :

@code_native aesenc(__m128i((213132, 13131)), __m128i((31231, 43213)))

    .section    __TEXT,__text,regular,pure_instructions
; ┌ @ In[36]:5 within 'aesenc'
    vaesenc %xmm1, %xmm0, %xmm0
    retq
    nopw    %cs:(%rax,%rax)
; └

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

Inline small functions

Рассмотрим ассемблирование такой функции:

# Simply throw an error
f() = error()
@code_native f()

    .section    __TEXT,__text,regular,pure_instructions
; ┌ @ In[38]:2 within 'f'
    pushq   %rax
    movabsq $inlinе$error, %rax
    callq   *%rax
    nopl    (%rax)
; └

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

Однако глянем на низкоуровневый код этой функции:

call_plus(x) = x + 1
code_native(call_plus, (Int,), debuginfo=:none)

    .section    __TEXT,__text,regular,pure_instructions
    leaq    1(%rdi), %rax
    retq
    nopw    %cs:(%rax,%rax)
    nop

Функция call_plus вызывает + и компилируется в одну инструкцию leaq (а также некоторые заполнители retq и nopw). Но + — это нормальная функция Julia, поэтому call_plus — это пример того, как одна обычная функция Julia вызывает другую. Почему нет инструкции callq для вызова +?

Компилятор решил встроить функцию + в call_plus. Это означает, что вместо вызова + он скопировал содержимое + непосредственно в call_plus. Преимущества этого заключаются в следующем:

• Нет никаких накладных расходов от вызова функции
• Нет необходимости создавать "Кортеж" для хранения аргументов функции +.
• Все вычисления, происходящие в +, компилируются вместе с call_plus, что позволяет компилятору использовать информацию из одного в другом и, возможно, упростить некоторые вычисления.

Так почему же тогда не все функции встроены? Инлайнинг копирует код, увеличивает его размер и потребляет оперативную память. Кроме того, сами инструкции ЦП также должны помещаться в кэш ЦП (хотя инструкции ЦП имеют свой собственный кэш), чтобы быть эффективно извлеченными. Если бы все было встроено, программы стали бы огромными. Инлайнинг хорош, когда встроенная функция мала.

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

Если код содержит чувствительный ко времени раздел, например внутренний цикл, важно посмотреть на код сборки, чтобы убедиться, что небольшие функции в цикле встроены. Например, в этой строке в моем хешировании kmer code, общая производительность minhashing падает в два раза, если эта аннотация @inline удаляется.

Крайняя разница между инлайнингом и отсутствием инлайнинга может быть продемонстрирована таким образом:

@noinline noninline_poly(x) = x^3 - 4x^2 + 9x - 11
inline_poly(x) = x^3 - 4x^2 + 9x - 11

function time_function(F, x::AbstractVector)
    n = 0
    for i in x
        n += F(i)
    end
    return n
end;

@btime time_function(noninline_poly, data)
@btime time_function(inline_poly, data);

#   13.140 μs (1 allocation: 16 bytes)
#   5.951 μs (1 allocation: 16 bytes)

Unroll tight loops

Рассмотрим функцию, которая суммирует вектор 64-разрядных целых чисел. Если смещение памяти вектора данных хранится в регистре %r9, длина вектора сохраняется в регистре %r8, текущий индекс вектора в %rcx, а накапливание в %rax, то код вложенного цикла может выглядеть следующим образом:

L1:
    ; add the integer at location %r9 + %rcx * 8 to %rax
    addq   (%r9,%rcx,8), %rax

    ; increment index by 1
    addq   $inlinе$1, %rcx

    ; compare index to length of vector
    cmpq   %r8, %rcx

    ; repeat loop if index is smaller
    jb     L1

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

Однако, если функция написана следующим образом:

function sum_vector(v::Vector{Int})
    n = 0
    i = 1
    for chunk in 1:div(length(v), 4)
        n += v[i + 0]
        n += v[i + 1]
        n += v[i + 2]
        n += v[i + 3]
        i += 4
    end
    return n
end

Результат, очевидно, тот же самый, если мы предположим, что длина вектора делится на четыре. Если длина не делится на четыре, мы можем просто использовать приведенную выше функцию для суммирования первых N — rem(N, 4) элементов и добавления последних нескольких элементов в другой цикл. Несмотря на функционально идентичный результат, код цикла отличается и может выглядеть примерно так:

L1:
    addq   (%r9,%rcx,8), %rax
    addq   8(%r9,%rcx,8), %rax
    addq   16(%r9,%rcx,8), %rax
    addq   24(%r9,%rcx,8), %rax
    addq   $inlinе$4, %rcx
    cmpq   %r8, %rcx
    jb     L1

В общей сложности 7 инструкций на 4 суммирования, или 1,75 инструкции на операцию. Это меньше половины количества инструкций на число integer! Увеличение скорости происходит от того, что проверка выхода заграницы цикла происходит реже. Мы называем этот процесс разворачиванием цикла, в данном случае в четыре раза. Естественно, разворачивание может быть выполнено только в том случае, если мы заранее знаем число итераций, поэтому мы не "перескакиваем" число итераций. Часто компилятор автоматически разворачивает циклы для повышения производительности, но иногда стоит посмотреть на ассемблерный код. Например, вот сборка для самого внутреннего цикла, сгенерированного на моем компьютере для sum([1]):

L144:
    vpaddq  16(%rcx,%rax,8), %ymm0, %ymm0
    vpaddq  48(%rcx,%rax,8), %ymm1, %ymm1
    vpaddq  80(%rcx,%rax,8), %ymm2, %ymm2
    vpaddq  112(%rcx,%rax,8), %ymm3, %ymm3
    addq    $inlinе$16, %rax
    cmpq    %rax, %rdi
    jne L144

Там где можно, он развернут в четыре раза и использует 256-битные инструкции SIMD, в общей сложности 128 байт, 16 целых чисел, суммированных за итерацию, или 0.44 инструкции на целое число.

Обратите также внимание, что компилятор выбирает использование 4 различных SIMD-регистров ymm, от ymm0 до ymm3, тогда как в моем примере ассемблерного кода я использовал только один регистр rax. Это происходит потому, что если вы используете 4 независимых регистра, то вам не нужно ждать завершения одного vpaddq (помните, что у него была задержка ~3 такта), прежде чем процессор сможет начать следующий.

Avoid unpredictable branches

Как упоминалось ранее, инструкции ЦП занимают по несколько циклов, но могут быть поставлены в очередь в ЦП до того, как предыдущая инструкция завершит вычисление. Итак, что же происходит, когда процессор сталкивается с ветвлением (то есть с инструкцией перехода)? Он не может знать, какую инструкцию поставить в очередь следующей, потому что это зависит от инструкции, которую он только что поставил в очередь и которая еще не выполнена.

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

Для программиста это означает, что скорость if-оператора зависит от того, насколько легко его угадать. Если это тривиально, предсказатель ветвления будет правильным почти все время, и оператор if займет не больше, чем простая инструкция, обычно 1 такт. В ситуации, когда ветвление является случайным, оно будет неверным примерно в 50% случаев, и каждое неверное предсказание может стоить около 10 тактов.

Ветви, вызванные циклами, являются одними из самых простых для угадывания. Если у вас есть цикл с 1000 элементами, код повторит цикл 999 раз и вырвется из цикла только один раз. Следовательно, предсказатель ветвей может просто всегда предсказать "опять по циклу" и получить точность 99,9%.

Мы можем продемонстрировать накладки неправильного предсказания ветвей с помощью простой функции:

# Copy all odd numbers from src to dst.
function copy_odds!(dst::Vector{UInt}, src::Vector{UInt})
    write_index = 1
    @inbounds for i in eachindex(src) # <--- this branch is trivially easy to predict
        v = src[i]
        if isodd(v)  # <--- this is the branch we want to predict
            dst[write_index] = v
            write_index += 1
        end
    end
    return dst
end

dst = rand(UInt, 5000)
src_random = rand(UInt, 5000)
src_all_odd = [2i+1 for i in src_random];

@btime copy_odds!(dst, src_random)
@btime copy_odds!(dst, src_all_odd);

#  10.977 μs (0 allocations: 0 bytes)
#  1.700 μs (0 allocations: 0 bytes)

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

Обратите внимание, что если вы используете меньшие векторы и повторяете вычисления много раз, как это делает макрос @btime, предсказатель ветвей способен выучить шаблон малых случайных векторов наизусть и достигнет гораздо лучшего результата, чем случайное предсказание. Это особенно заметно в самых современных процессорах (и в частности, в процессорах, продаваемых AMD, как я слышал), где предсказатели ветвей стали намного лучше. Это "обучение наизусть" является артефактом цикла в процессе бенчмаркинга. Вы не ожидали бы многократного выполнения одних и тех же вычислений на реальных данных:

src_random = rand(UInt, 100)
src_all_odd = [2i+1 for i in src_random];

@btime copy_odds!(dst, src_random)
@btime copy_odds!(dst, src_all_odd);

#   63.647 ns (0 allocations: 0 bytes)
#   43.105 ns (0 allocations: 0 bytes)

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

Если ветви не могут быть легко предсказаны, часто можно перефразировать функцию, чтобы избежать ветвей. Например, в приведенном выше примере copy_odds! мы могли бы вместо этого написать так:

function copy_odds!(dst::Vector{UInt}, src::Vector{UInt})
    write_index = 1
    @inbounds for i in eachindex(src)
        v = src[i]
        dst[write_index] = v
        write_index += isodd(v)
    end
    return dst
end

dst = rand(UInt, 5000)
src_random = rand(UInt, 5000)
src_all_odd = [2i+1 for i in src_random];

@btime copy_odds!(dst, src_random)
@btime copy_odds!(dst, src_all_odd);

#   2.019 μs (0 allocations: 0 bytes)
#   1.985 μs (0 allocations: 0 bytes)

Тут нет никаких других ветвей, кроме той, которая вызвана самим циклом (что легко предсказуемо), и приводит к скорости несколько хуже, чем идеально предсказанная, но гораздо лучше для случайных данных.

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

Переменная тактовая частота

Современный процессор ноутбука, оптимизированный для низкого энергопотребления, потребляет примерно 25 Вт мощности на чипе размером со штамп (и тоньше человеческого волоса). Без надлежащего охлаждения это приведет к тому, что температура процессора взлетит до небес и расплавит пластик чипа, разрушив его. Как правило, процессоры имеют максимальную рабочую температуру около 100 C. Энергопотребление, а следовательно, и тепловыделение зависит, помимо прочего, от тактовой частоты — более высокие тактовые частоты генерируют больше тепла.

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

Как программист, вы мало что можете сделать, чтобы учесть температуру процессора, но это полезно знать. В частности, колебания температуры процессора часто объясняют наблюдаемую разницу в производительности:

• Процессоры обычно работают быстрее всего в начале рабочей нагрузки, а затем снижают производительность по мере достижения максимальной температуры
• Инструкции SIMD обычно требуют больше энергии, чем обычные инструкции, генерируя больше тепла и снижая тактовую частоту. Это может компенсировать некоторый прирост производительности SIMD, но SIMD все равно всегда будет более эффективным, когда это применимо

Multithreading

В старые добрые времена тактовая частота процессора увеличивалась с каждым годом по мере появления на рынке новых процессоров. Частично из-за тепловыделения это ускорение замедлилось, как только процессоры достигли отметки в 3 ГГц. Теперь мы видим лишь незначительные приращения тактовой частоты в каждом поколении процессоров. Вместо необработанной скорости выполнения, фокус сместился на получение большего количества вычислений за такт. Кэши ЦП, конвейеризация ЦП, прогнозирование ветвей и инструкции SIMD — все это важные достижения в этой области, и все они были рассмотрены выше.

Еще одна важная область, в которой процессоры улучшились — это просто их число: почти все процессорные чипы содержат несколько меньших процессоров или ядер внутри. Каждое ядро имеет свой собственный небольшой кэш процессора и выполняет вычисления параллельно. Кроме того, многие процессоры имеют функцию под названием hyper-threading, где два потока (т.е. потоки инструкций) могут работать на каждом ядре. Идея состоит в том, что всякий раз, когда один процесс останавливается (например, из-за того, что он испытывает промах кэша или неправильное предсказание), другой процесс может продолжаться на том же ядре. Процессор "притворяется", что у него в два раза больше процессоров. Например, я пишу это на ноутбуке с процессором Intel Core i9-9880H. Этот процессор имеет 8 ядер, но различные операционные системы, такие как Windows или Linux, будут показывать 16 "процессоров" в программе системного монитора.

Гиперпоточность действительно имеет значение только тогда, когда ваши потоки иногда не могут выполнять работу. Помимо внутренних причин процессора, таких как промахи кэша, поток также может быть приостановлен, потому что он ожидает внешнего ресурса, такого как веб-сервер или данные с диска. Если вы пишете программу, в которой некоторые потоки проводят значительное время на холостом ходу, ядро может быть использовано другим потоком, и hyperthreading может показать свою ценность.

Давайте посмотрим нашу первую параллельную программу в действии. Во-первых, нам нужно убедиться, что Джулия действительно была запущена с правильным количеством потоков. Вы можете установить переменную окружения JULIA_NUM_THREADS перед запуском Julia. У меня есть 8 ядер на этом процессоре, все с гиперпоточностью, поэтому я установил количество потоков на 16:

Threads.nthreads() # 16

# Spend about half the time waiting, half time computing
function half_asleep(start::Bool)
    a, b = 1, 0
    for iteration in 1:5
        start && sleep(0.06)
        for i in 1:100000000
            a, b = a + b, a
        end
        start || sleep(0.06)
    end
    return a
end

function parallel_sleep(n_jobs)
    jobs = []
    for job in 1:n_jobs
        push!(jobs, Threads.@spawn half_asleep(isodd(job)))
    end
    return sum(fetch, jobs)
end

parallel_sleep(1); # run once to compile it

for njobs in (1, 4, 8, 16, 32)
    @time parallel_sleep(njobs);
end

  0.595490 seconds (48 allocations: 1.969 KiB)
  0.560908 seconds (371 allocations: 21.062 KiB)
  0.573709 seconds (422 allocations: 20.875 KiB)
  0.649239 seconds (697 allocations: 32.984 KiB)
  2.157162 seconds (1.25 k allocations: 57.156 KiB)

Видно, что с помощью этой задачи мой компьютер может выполнять 16 заданий параллельно почти так же быстро, как он может выполнять 1. Но 32 задания занимают гораздо больше времени.

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

Распараллеливаемость

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

Необходимым условием эффективного использования многопоточности является то, что ваши вычисления могут быть разбиты на несколько блоков, над которыми можно работать независимо. К счастью, большинство вычислительных задач (по крайней мере, в моей области работы, биоинформатике) содержат подзадачи, которые до неприличия параллельны. Это означает, что существует естественный и простой способ разбить его на подзадачи, которые могут быть обработаны независимо. Например, если для 100 генов требуется определенное независимое вычисление, естественно использовать один поток для каждого гена.

Давайте возьмем пример небольшой, до неприличия параллельной задачи. Мы хотим построить Julia set. Оно названо в честь Гастона Жюлиа и не имеет ничего общего с языком Джулия. Множество Жюлиа — это (часто) фрактальные наборы комплексных чисел. Сопоставляя реальный и комплексный компонент членов множества со значениями пикселей X и Y экрана, можно генерировать ЛСД-Триппи-изображения, связанные с фракталами.

Множество Жюлиа, который я создаю ниже, определяется следующим образом: мы определяем функцию , где — некоторая константа. Затем мы записываем количество раз, когда может быть применено к любому заданному комплексному числу до . Количество итераций соответствует яркости одного пикселя на изображении. Мы просто повторяем это для диапазона реальных и мнимых значений в сетке, чтобы создать изображение.

Во-первых, давайте посмотрим на непараллельное решение:

const SHIFT = Complex{Float32}(-0.221, -0.713)

f(z::Complex) = z^2 + SHIFT

"Set the brightness of a particular pixel represented by a complex number"
function mandel(z)
    n = 0
    while ((abs2(z) < 4) & (n < 255))
        n += 1
        z = f(z)
    end
    return n
end

"Set brightness of pixels in one column of pixels"
function fill_column!(M::Matrix, x, real)
    for (y, im) in enumerate(range(-1.0f0, 1.0f0, length=size(M, 1)))
        M[y, x] = mandel(Complex{Float32}(real, im))
    end
end

"Create a Julia fractal image"
function julia()
    M = Matrix{UInt8}(undef, 20000, 5000)
    for (x, real) in enumerate(range(-1.0f0, 1.0f0, length=size(M, 2)))
        fill_column!(M, x, real)
    end
    return M
end;

@time M = julia();
#   9.565873 seconds (2 allocations: 95.368 MiB, 0.09% gc time)

Заняло около 10 секунд на моем компьютере. Теперь перейдем к параллельному:

function recursive_fill_columns!(M::Matrix, cols::UnitRange)
    F, L = first(cols), last(cols)
    # If only one column, fill it using fill_column!
    if F == L
        r = range(-1.0f0,1.0f0,length=size(M, 1))[F]
        fill_column!(M, F, r)
    # Else divide the range of columns in two, spawning a new task for each half
    else
        mid = div(L+F,2)
        p = Threads.@spawn recursive_fill_columns!(M, F:mid)
        recursive_fill_columns!(M, mid+1:L)
        wait(p)
    end
end

function julia()
    M = Matrix{UInt8}(undef, 20000, 5000)
    recursive_fill_columns!(M, 1:size(M, 2))
    return M
end;

@time M = julia();
#   0.564578 seconds (35.57 k allocations: 98.856 MiB, 3.72% gc time)

Почти в 16 раз быстрее! С 16 потоками это близко к наилучшему сценарию, возможному только для почти идеальных смущающе параллельных задач.

Несмотря на потенциал больших выгод, на мой взгляд, многопоточность должна быть одним из последних средств повышения производительности по трем причинам:

1. реализация многопоточности во многих случаях сложнее, чем другие методы оптимизации. В приведенном примере это было очень просто. В сложном рабочем процессе можно быстро запутаться.
2. многопоточность может привести к труднодиагностируемым и неустойчивым ошибкам. Они почти всегда связаны с несколькими потоками, считывающими и записывающими в одну и ту же память. Например, если два потока одновременно увеличивают целое число со значением N, то оба потока будут считывать N из памяти и записывать N+1 обратно в память, где правильный результат двух приращений должен быть N+2! Бесит, что эти ошибки появляются и исчезают непредсказуемо, так как они вызваны неудачным временем. Эти ошибки, конечно, имеют решения, но это сложная тема, выходящая за рамки данного опуса.
3. наконец, достижение производительности с помощью нескольких потоков — это действительно достижение производительности за счет потребления большего количества ресурсов, а не получение чего-то из ничего. Часто вы платите за использование большего количества потоков, либо буквально покупая время облачных вычислений, либо оплачивая счет за повышенное потребление электроэнергии из нескольких ядер процессора, либо метафорически претендуя на большее количество ресурсов процессора ваших пользователей, которые они могли бы использовать где-то еще. Напротив, более эффективные вычисления ничего не стоят.

GPUs

До сих пор мы рассматривали только самый важный вид вычислительных чипов — центральный процессор. Но есть много других видов чипов. Наиболее распространенным видом альтернативного чипа является графический процессор (GPU).

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

Графические процессоры пожертвовали многими наворотами процессоров, описанных в этом документе, такими как специализированные инструкции, SIMD и прогнозирование ветвей. Они также обычно работают на более низких частотах, чем CPU. Это означает, что их необработанная вычислительная мощность во много раз медленнее, чем у процессора. Чтобы компенсировать это, у них есть большое количество ядер. Например, высокопроизводительный игровой графический процессор NVIDIA RTX 2080Ti имеет 4.352 ядра. Следовательно, некоторые задачи могут испытывать 10-кратное или даже 100-кратное ускорение с помощью графического процессора. В частности, для научных приложений матричные и векторные операции обладают высокой степенью распараллеливаемости.

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

Существуют также более эзотерические чипы, такие как TPU (явно предназначенные для низкоточных тензорных операций, распространенных в глубоком обучении) и ASICs (общий термин для узкоспециализированных чипов, предназначенных для одного приложения). На момент написания статьи эти чипы были необычными, дорогими, плохо поддерживаемыми и имели ограниченное применение, и поэтому не представляли никакого интереса для исследователей, не связанных с компьютерными науками.