Parsing Integers with SIMD - HackMD

<style> .reveal { font-size: 24px; } .reveal section img { background:none; border:none; box-shadow:none; } pre.graphviz { box-shadow: none; } .slide-background-content { background: #fff; } code { color: #900 !important; background: #eee !important; border: 1px solid #ddd; border-radius: 3px; /* padding: 0 2px; */ } .code-wrapper { font-size: 13px !important; box-shadow: none !important; } a { font-size: 0.9em; } h1, h2, h3 { /* text-transform: none !important; */ } .hljs{display:block;overflow-x:auto;padding:.5em;color:#333;background:#eee}.hljs-comment,.hljs-quote{color:#998;font-style:italic}.hljs-keyword,.hljs-selector-tag,.hljs-subst{color:#333;font-weight:700}.hljs-literal,.hljs-number,.hljs-tag .hljs-attr,.hljs-template-variable,.hljs-variable{color:teal}.hljs-doctag,.hljs-string{color:#d14}.hljs-section,.hljs-selector-id,.hljs-title{color:#900;font-weight:700}.hljs-subst{font-weight:400}.hljs-class .hljs-title,.hljs-type{color:#458;font-weight:700}.hljs-attribute,.hljs-name,.hljs-tag{color:navy;font-weight:400}.hljs-link,.hljs-regexp{color:#009926}.hljs-bullet,.hljs-symbol{color:#990073}.hljs-built_in,.hljs-builtin-name{color:#0086b3}.hljs-meta{color:#999;font-weight:700}.hljs-deletion{background:#fdd}.hljs-addition{background:#dfd}.hljs-emphasis{font-style:italic}.hljs-strong{font-weight:700} </style> ## Parsing Integers with SIMD Sergey Slotin C++ Zero Cost Conf 2023 Note: Привет, меня зовут Сергей Слотин Некоторые из вас знают меня или возможно читали мои статьи по алгоритмам в интернете. Я занимаюсь тем, что беру базовые алгоритмы и пытаюсь их соптимизировать. В прошлом году я имел часть выступать на этой замечательной конференции с рассказом о том как соптимизировать бинарный поиск. Моя была не рассказать про сам бинарный поиск, а рассказать про то как работает пайплайн процессора, память и всё такое. А сегодня я хотел бы погрузить вас в замечательный мир SIMD-программирования ---  ```cpp // parse n integers from stdin into a void parse_integers(int *a, int n) { // ... } ``` ```python # generate 100M non-negative 32-bit integers (~1GB) from random import randint n = 10**8 for _ in range(n): print(randint(0, 2**31 - 1)) ``` ```bash python generator.py > input.txt g++ -std=c++17 -O3 -march=native -o run time ./run < input.txt ``` All benchmarking will be done with GCC 12.2 on Zen 2 @ 2GHz The input file is cached (RAM bandwidth ~8GB) Note: И рассказывать я буду на примере другой кажущейся абсурдной задаче — мы будем парсить числа. --- ```cpp void parse_integers(int *a, int n) { for (int i = 0; i < n; i++) std::cin >> a[i]; } ``` 50s total → 500ns per number 20MB/s → 100 CPU cycles per byte --- ```cpp void parse_integers(int *a, int n) { std::ios_base::sync_with_stdio(false); std::cin.tie(NULL); for (int i = 0; i < n; i++) std::cin >> a[i]; } ``` 24 cycles per byte --- ```cpp void parse_integers(int *a, int n) { std::ios_base::sync_with_stdio(false); std::cin.tie(NULL); for (int i = 0; i < n; i++) { unsigned int x = 0; std::cin >> x; a[i] = x; } } ``` 24 → 22 cycles per byte --- ```cpp void parse_integers(int *a, int n) { for (int i = 0; i < n; i++) scanf("%d", &a[i]); } ``` 40 cycles per byte --- ```cpp void parse_integers(int *a, int n) { for (int i = 0; i < n; i++) scanf("%u", &a[i]); // ^ unsigned } ``` 40 → 30 cycles per byte --- ```cpp void parse_integers(int *a, int n) { int k = 0; // how many we've parsed int x = 0; // current number while (k < n) { int c = getchar(); if (c >= '0'/* && c <= '9'*/) { x = 10 * x + (c - '0'); } else { a[k++] = x; x = 0; } } } ``` 12.6 cycles per byte --- ```text $ man getchar_unlocked NAME getc_unlocked, getchar_unlocked, putc_unlocked, putchar_unlocked - non‐locking stdio functions DESCRIPTION Each of these functions has the same behavior as its counterpart with‐ out the "_unlocked" suffix, except that they do not use locking (they do not set locks themselves, and do not test for the presence of locks set by others) and hence are thread-unsafe. ``` Note: В зависимости от реализации и вашей удачи, что-то может потеряться, не прочитаться или прочитаться дважды. Но у нас всего один поток, и можно не бояться. --- ```cpp void parse_integers(int *a, int n) { int k = 0; int x = 0; while (k < n) { int c = getchar_unlocked(); if (c >= '0') { x = 10 * x + (c - '0'); } else { a[k++] = x; x = 0; } } } ``` 11.2 cycles per byte   --- ```text $ man fread NAME fread, fwrite - binary stream input/output SYNOPSIS size_t fread(void ptr[restrict .size * .nmemb], size_t size, size_t nmemb, FILE *restrict stream); DESCRIPTION The function fread() reads nmemb items of data, each size bytes long, from the stream pointed to by stream, storing them at the location given by ptr. ``` There is also `fread_unlocked`, but taking lock once per block does not make a difference --- ```cpp void parse_integers(int *a, int n) { const int BUFFER_SIZE = (1 << 14); // 16K char buffer[BUFFER_SIZE]; int k = 0; int x = 0; while (k < n) { int parsed = fread(buf, 1, BUFFER_SIZE, stdin); for (int i = 0; i < parsed; i++) { char c = buf[i]; if (c >= '0') { x = 10 * x + (c - '0'); } else { a[k++] = x; x = 0; } } } } ``` 4.6 cycles per byte --- ```cpp $ man mmap NAME mmap — map pages of memory SYNOPSIS void *mmap(void *addr, size_t len, int prot, int flags, int fildes, off_t off); DESCRIPTION The mmap() function shall establish a mapping between an address space of a processand a memory object. ``` --- ```cpp void parse_integers(int *a, int n) { int fsize = lseek(0, 0, SEEK_END); auto input = (char*) mmap(0, fsize, PROT_READ, MAP_PRIVATE, 0, 0); int k = 0; int x = 0; for (int i = 0; i < fsize; i++) { char c = input[i]; if (c >= '0') x = 10 * x + (c - '0'); else { a[k++] = x; x = 0; } } } ``` 3.8 cycles per byte (multiplication + addition takes 3+1 = 4 cycles) Note, this won't work unless stdin is a file redirect --- ![](https://hackmd.io/_uploads/ryTzweh93.png =800x) Note: Ну ладно, с тем, как читать файлы разобрались, давайте перейдем к содержательной части --- ![](https://upload.wikimedia.org/wikipedia/commons/thumb/2/21/SIMD.svg/1200px-SIMD.svg.png =350x) SIMD: <u>S</u>ingle <u>I</u>nstruction, <u>M</u>ultiple <u>D</u>ata Vector registers holding 128, 256, or 512 bits of data logically split into blocks of 8, 16, 32, or 64 bits  Note: SIMD это сокращение от Single Instruction Multiple Data. Это такой подход в архитектуре компьютеров, где одна операция применяется параллельно сразу к блоку данных. (click) Он особенно популярен на GPU, где задачи очень простые и параллельные, но последние лет 20 есть и на обычных процессорах — как мощных серверных, так и десктопных и мобильных, в немного разных редакциях. --- Layers of SIMD abstractions: *  x86 assembly *  C/C++ intrinsics *  Built-in vector types *  SIMD libraries (VCL, EVE, Highway) *  Sometime in the next decade: std::simd *  SPMD compilers (OpenMP, ISPC) *  Auto-vectorization Later are easier, former give more control Note: Когда программа переписывается с обычного скалярного кода на код, использующий векторные инструкцией, это называют векторизацией, и есть несколько способов это сделать, отличающиеся по трудоёмкости и широте возможностей. Во-первых, векторные инструкции можно использовать так же, как и обычные — напрямую через ассемблер. Никому последние 50 лет не нравится программировать на ассемблере кроме как для совсем performance-critical вещей, и поэтому в достаточно низкоуровневых языках вроде C и C++ есть более элегантные способы их использовать. Есть специальные векторные типы данных, в которые вмещаются 128, 256 или 512 бит и интринзики (сишные функции, принимающие и возвращающие векторные типы), которые дают почти прямой доступ к этим инструкциям. Они используют конкретный набор инструкций и не очень хорошо портируются между архитектурами. В отдельных компиляторах есть чуть более высокоуровневая абстракция для векторных типов, которая позволяет применять простые операции вроде поэлементно сложить два вектора из 8 чисел с помощью перегруженных операторов (я дальше покажу на примере). Дальше, есть библиотеки, которые реализуют свои портируемые векторные типы и по-разному пытаются предоставить доступ к инструкциям. И на ещё более высоком уровне есть разные API основанные на парадигме SPMD (single program, multiple data) -- это когда ты пишешь обычный скалярный код который что-то делает с единичными элементами, а не целыми векторами, но как бы подразумневая что эту операцию нужно будет применить ко всему массиву. Но чаще всего программисты сталкиваются с SIMD, даже это не осознавая. В достаточно простых и часто встречающихся случаях, оптимзирующий компилятор может векторизовать код самостоятельно. --- ```cpp const int n = 100000; int a[n], s = 0; int main() { for (int t = 0; t < 100000; t++) for (int i = 0; i < n; i++) s += a[i]; return 0; } ``` `g++ -O3 sum.cc -o run && time ./run` 1.26 seconds to perform $10^{10}$ operations 8B additions per second, or 4 per cycle Note: Для примера рассмотрим такой код, в котором в цикле сто тысяч раз суммируют сто тысяч интов. Если его скомпилировать с GCC с самым агрессивным уровнем оптимизаций, от завершается за 1.26 секунды. Сто тысяч на сто тысяч это суммарно 10^10 операций, если поделить это на время исполнения, то получается 8 миллиардов операций в секунду, или, так как у нас частота процессора фиксированная 2 гигагерца, то есть 2 миллиарда циклов в секунду, получается 4 операции сложения на цикл. Как вы уже догадались, это получается за счёт того, что компилятор автоматически векторизовал цикл, используя инструкции, которые позволяют за раз обработать блок из 128 бит, в который помещаются 4 инта. --- ![](https://i0.wp.com/www.urtech.ca/wp-content/uploads/2017/11/Intel-mmx-sse-sse2-avx-AVX-512.png =500x) Backwards-compatible up until AVX-512 Note: ... --- ```cpp #pragma GCC target("avx2") // (the rest is the same) const int n = 100000; int a[n], s = 0; int main() { for (int t = 0; t < 100000; t++) for (int i = 0; i < n; i++) s += a[i]; return 0; } ``` You can also add `-mavx` or `-march=native` flags 0.63 seconds — twice as fast (!) Note: ... Но компилятор может так делать только в самых простых случаях. На самом деле, он даже не справляется соптимизировать сумму массива, но об этом чуть попозже. Чтобы делать чтото сложное, нужно научиться пользоваться интринзиками. --- All C/C++ intrinsics can be included with `x86intrin.h` ```cpp #include <x86intrin.h> ``` (make sure you add relevant compilation flags or it will crash) Note: ... --- ```cpp double a[100], b[100], c[100]; // iterate in blocks of 4, // because that's how many doubles can fit into a 256-bit register for (int i = 0; i < 100; i += 4) { // load two 256-bit segments into registers __m256d x = _mm256_loadu_pd(&a[i]); __m256d y = _mm256_loadu_pd(&b[i]); // add 4+4 64-bit numbers together __m256d z = _mm256_add_pd(x, y); // write the 256-bit result into memory, starting with c[i] _mm256_storeu_pd(&c[i], z); } ``` Note: ... --- ```cpp for (int i = 0; i < n; i += B) { // do something with a block of B elements } ``` If $n$ is not divisible by the block size $B$, we have two options: 1.  Pad the array with neutal elements (e. g. zeros) 2.  Break the loop before the last block and proceed normally Humans prefer #1, compilers prefer #2  Note: One large disadvantage of SIMD is that you need to get data in vectors, and sometimes the memory layout or the specifics of the problem make it hard. In the array addition example, and for loop vectorization in general, we split the array into small blocks that we can process with SIMD instructions and iterate over these blocks. We have actually ignored the fact that the array length may not be perfectly divisible by the SIMD block size, and this may be a problem. In this case, we have two options. --- ```cpp typedef __m256i reg; // <- this typedef will be frequently used void sum(int *a, int *b, int *c, int n) { // vectorized part for (int i = 0; i + 7 < n; i += 8) { // ^ careful to not overchop reg x = _mm256_loadu_pd(&a[i]); reg y = _mm256_loadu_pd(&b[i]); reg z = _mm256_add_epi32(x, y); _mm256_storeu_pd(&c[i], z); } // scalar part for (int i = n / 8 * 8; i < n; i++) c[i] = a[i] + b[i]; } ``` (We will omit scalar remainders for brevity) Note: Я эту часть просто буду пропускать. --- - 128-bit `__m128`, `__m128d` and `__m128i` types for `float`, `double` and `int` - 256-bit `__m256`, `__m256d`, `__m256i` - 512-bit `__m512`, `__m512d`, `__m512i` - Typedefs and templates can be handy: `typedef __m256i reg` Note: There are different flavors, that specify the size of a vector and what kind of data is inside: a 32-bit single float, a 64-bit double float or some kind of integer data. These types exist only for type checking -- the actual physical registers are the same and the processor doesn't really differentiate between them. You can freely convert between these types with C-style casting without any cost. It is tedious to type these names, so we will mostly be using a typedef from here on and only work with 256-bit integer vectors unless explicitly specified otherwise. --- Most SIMD intrinsics follow the `_mm<size>_<action>_<type>` naming convention: - `_mm_add_epi16`: add two 128-bit vectors of 16-bit *extended packed integers* - `_mm256_acos_pd`: calculate elementwise $\arccos$ for 4 *packed doubles* - `_mm256_ceil_ps`: round up each of 8 `float`s to the nearest integer - `_mm256_broadcast_sd`: copy a `double` from memory to 4 elements of a vector - `_mm256_cmpeq_epi32`: compare 8+8 packed `int`s and return a 256-bit mask - `_mm256_blendv_ps`: pick elements from one of two vectors according to a mask - `_mm256_permutevar8x32_epi32`: pick elements by indices from a vector Note: Most SIMD intrinsics follow a naming convention similar to `_mm<size>_<action>_<type>` and correspond to a single similarly named assembly instruction. There are a lot of them; here are some examples. (explain intrinsics) The naming can be a bit confusing. I've spent hundreds of hours writing SIMD code with intrinsics and I still can't remember whether it is underscore-m-m (used for intrinsics) or underscore-underscore-m (used for types). Also, even when accounting for versions of an instruction with different data types, there is still *a lot* of them… --- ![](https://hackmd.io/_uploads/rJfWughc2.png =300x) https://db.in.tum.de/~finis/x86%20intrinsics%20cheat%20sheet%20v1.0.pdf --- ![](https://i.imgur.com/ZIzDidV.png =600x) software.intel.com/sites/landingpage/IntrinsicsGuide/ Note: …so intel created a very helpful reference, called intel intrinsics guide. (describe reference) Intel being Intel doesn't include timings for AMD CPUs, but you can look them up in other instruciton tables Если вы в России, то этот сайт не откроется, потому что санкции, интел ушел из России и показывают заглушку на своем сайте, включая девтулзы. Но вы можете нагуглить оффлайн-версию или открывать сайт через vpn. --- ### Reductions ```cpp int sum(int *a, int n) { int s = 0; for (int i = 0; i < n; i++) s += a[i]; // <- dependency between iterations return s; } ``` * Calculating A + B is easy because there are no data dependencies * Calculating array sum is different: we need the accumulator from the previous step Note: Окей, разобрались с семантикой, научились скалыдвать два массива, давайте теперь сделаем что-то более сложное — посчитаем, собственно, сумму массива. Это сложнее, потому что здесь цикл зависит от предыдущей итерации. Переменная s должна обновиться до того, как мы можем начать следующую итерацию. --- ```cpp // the sum of elements in an 8-int vector int hsum(reg x) { // ... } int sum(int *a, int n) { reg s = _mm256_setzero_si256(); // "vertical summation" for (int i = 0; i < n; i += 8) { reg x = _mm256_load_si256((reg*) &a[i]); s = _mm256_add_epi32(s, x); } return hsum(s); // "horizontal summation" } ``` We can calculate $B$ partial sums $\{i+kB\}$ for each $i<B$ and sum them up This trick works with any other commutative operator (e. g. `min`)  Note: We can solve that by logically splitting the array into 8 interleaved partitions, calculate the 8 sums independently with SIMD, and then sum up the partial sums. Here we create an accumulator variable that is an 8-element vector (initially filled with zeros) that we use to sum up each block. This part is called vertical summation, and then we need to sum up the partial sums by just dumping the vector somewhere in the memory and performing the 8 additions manually -- this is called horizontal summation. There are faster ways to do the horizontal summation, but this part is executed only once so it doesn't really matter that much. The loop can be further optimized, but we need to focus on a different aspect. --- ```cpp int hsum(reg x) { int t[8], s = 0; _mm256_storeu_si256((reg*) t, x); for (int i = 0; i < 8; i++) s += t[i]; return s; } ``` --- ![](https://www.codeproject.com/KB/cpp/874396/Fig1.jpg =400x) "Horizontal addition" --- ```cpp int hsum(__m256i x) { __m128i l = _mm256_extracti128_si256(x, 0); __m128i h = _mm256_extracti128_si256(x, 1); l = _mm_add_epi32(l, h); l = _mm_hadd_epi32(l, l); return _mm_extract_epi32(l, 0) + _mm_extract_epi32(l, 1); } ``` Note: ... --- ### Instruction-level parallelism ```cpp int s = 0; for (int i = 0; i < n; i++) s += a[i]; ``` (Pretend vectorization doesn't exist) Note: Modern CPUs are superscalar: they don't just execute an instruction sequence one by one. Instead, they read the sequence of instructions ahead of time and try to schedule them in an optimal way, most of the time actually executing several instructions concurrently if that is possible. Having this in mind, in general, when optimizing loops and other procedures for throughput, you need first to make sure that there is nothing stalling the execution, and then you need to remove pressure from the bottleneck, whatever it may be. Here, if we for simplicity look at how the scalar code works and pretend that vectorization doesn't exist, we can notice the following. ---- ```cpp int s = 0; s += a[0]; s += a[1]; s += a[2]; s += a[3]; // ... ``` The next iteration depends on the previous one, so we can't execute more than one iteration per cycle (despite `mov`/`add` having a throughput of two) Note: Even though all the operations involved have high enough throughput, the loop still can't go faster than one iteration per CPU cycle because the next iteration of the loop depends on the previous one ---- ```cpp int s0 = 0, s1 = 0; s0 += a[0]; s1 += a[1]; s0 += a[2]; s1 += a[3]; // ... int s = s0 + s1; ``` We can now proceed in two concurrent "threads", doubling the throughput Note: What we can do to resolve this contention is to similarly split the array into two parts -- odd and even -- and have two independent accumulators that we sum up only at the very end. This way we can do two additions per cycle and maximize (in this case double) the throughput of the entire loop ---- ```cpp int sum(int *a, int n) { reg s0 = _mm256_setzero_si256(); reg s1 = _mm256_setzero_si256(); for (int i = 0; i < N; i += 16) { reg x1 = _mm256_load_si256((reg*) &a[i]); reg x2 = _mm256_load_si256((reg*) &a[i + 8]); s0 = _mm256_add_epi32(s0, x1); s1 = _mm256_add_epi32(s1, x2); } return hsum(_mm256_add_epi32(s0, s1)); } ``` ~32B operations per second 2x faster than `std::accumulate` and naive loops (on GCC; recent LLVM vectorizes optimally) Note: We can apply the same trick in the vectorized version. Instead of one accumulator, we can use 2 vector accumulators. Our version is correspondingly two times faster than the previous and the auto-vectorized version (well, on GCC at least; the latest LLVM versions seem to be able to vectorize it optimally). This instruction-level parallelism trick is very powerfull and will be a recurring theme in this talk. --- ### Throughput bottleneck ![](https://hackmd.io/_uploads/rJK5O1n92.png =600x) Intel Intrinsics Guide Note: Можно на самом деле доказать --- ![](https://hackmd.io/_uploads/rkqCdy2ch.png =500x) https://www.agner.org/optimize/instruction_tables.pdf --- ![](https://hackmd.io/_uploads/rJp-tJ2ch.png =600x) https://uops.info/table.html --- ```cpp int sum(int *a, int n) { reg s0 = _mm256_setzero_si256(); reg s1 = _mm256_setzero_si256(); for (int i = 0; i < N; i += 16) { reg x1 = _mm256_load_si256((reg*) &a[i]); reg x2 = _mm256_load_si256((reg*) &a[i + 8]); s0 = _mm256_add_epi32(s0, x1); s1 = _mm256_add_epi32(s1, x2); } return hsum(_mm256_add_epi32(s0, s1)); } ``` 16 numbers processed per cycle, 32B operations per second Perfectly bottlenecked by loading (2/2), adding (2/2), and overall number of instructions (4/4)  --- ### Masking ```cpp for (int i = 0; i < N; i++) a[i] = rand() % 100; int s = 0; // branch: for (int i = 0; i < N; i++) if (a[i] < 50) s += a[i]; // no branch: for (int i = 0; i < N; i++) s += (a[i] < 50) * a[i]; // also no branch: for (int i = 0; i < N; i++) s += (a[i] < 50 ? a[i] : 0); ``` The main downside of SIMD is that there is no branching: we need to use predication Note: Before continuing doing interesting things, I need to teach you an important concept. Other than the need to store your data continuously and process it by chunks, SIMD has another major downside in that there is no branching: we need to use a technique called predication to replace it. Consider this example. Here we create a random array and then calculate the sum of all its elements that are under 50. If we would do it in the most straightforward way, we would write a loop with an if inside. But we can get rid of explicit branching with either one of these two tricks: multiplication by a bool or a ternary operator. This general trick is called predication. In SIMD, we have to do the same. --- - `_mm256_cmpgt_epi32`: compare the integers in two vectors and produce a mask of all ones all `>` comparisons and a mask of zeros for all `<=` comparisons Note: To perform predication in SIMD, we can use these instructions. --- ```cpp // for (int i = 0; i < N; i++) // s += (a[i] < 50 ? a[i] : 0); const reg c = _mm256_set1_epi32(50); const reg z = _mm256_setzero_si256(); reg s = _mm256_setzero_si256(); for (int i = 0; i < N; i += 8) { reg x = _mm256_load_si256( (reg*) &a[i] ); // [23, 71, 90, 40] reg mask = _mm256_cmpgt_epi32(c, x); // [-1, 0, 0, -1] = [23<50, 71<50, 90<50, 50<50] x = _mm256_blendv_epi8(z, x, mask); // [23, 0, 0, 40] // x = _mm256_and_si256(x, mask); // we can use bitwise "and" instead of blending s = _mm256_add_epi32(s, x); // [s1, s2, s3, s4] += [23, 0, 0, 40] } ``` Note: This example is still simple enough to be auto-vectorized. The compiler is capable of producing these 4 lines on its own - so let's move on to a more complex one. In this particular case, we can use bitwise "and" to do the blending. This is slightly faster because on this particular CPU the vector "and" takes one cycle fewer than "blend". --- - movemask - `_mm256_permutevar8x32_epi32`: select elements from a vector with given indices Despite the name, it doesn't *permute* values but just *copies* them to form a new vector (duplicates in the result are allowed) Note: This intrinsic with an intimidatingly long name takes two vector operands. The first one is data, and the second one is a vector of integers between 0 and 15 inclusive. These integers select, for each position of the output, which elements to pick of the data vector. This instruction thus can be used for arbitrary permutations of the data vector. --- ### Filtering ```cpp int a[N], b[N]; int filter() { int k = 0; for (int i = 0; i < N; i++) if (a[i] < P) b[k++] = a[i]; return k; } ``` 0.7-1.5 GFLOPS depending on `P` Note: Our final example is filtering an array -- that is, we need to construct another array that only has the elements that satisfy a given predicate (in this case, a simple comparison) in their original order. This is used, for example, in quicksort. The performance of the scalar version here varies because branch prediction success rate depends on the probability of this predicate --- - Calculate the predicate on an 8-element vector of data - Use the `movemask` instruction to get a scalar 8-bit mask - Use this mask to index a **lookup table**, returning the permutation moving the required elements to the beginning of the vector (in their original order) - Use the `_mm256_permutevar8x32_epi32` intrinsic to permute the values - Write the whole permuted vector to the output array - Popcount the mask and move the output pointer by that amount Note: Vectorizing this involves some precomputation. We need to use the fact that the permutation instruction can take arbitrary vector as selection indices. --- ```cpp struct Precalc { alignas(64) int permutation[256][8]; constexpr Precalc() : permutation{} { for (int m = 0; m < 256; m++) { int k = 0; for (int i = 0; i < 8; i++) if (m >> i & 1) permutation[m][k++] = i; } } }; constexpr Precalc T; ``` Note: First, we need to precompute the permutations. For each possible 8-bit mask, we need a permutation that compresses this vector by removing the elements that we do not need and collapsing all elements satisfying the predicate to the beginning of the vector --- ```cpp const reg p = _mm256_set1_epi32(P); int filter() { int k = 0; for (int i = 0; i < N; i += 8) { reg x = _mm256_load_si256( (reg*) &a[i] ); reg m = _mm256_cmpgt_epi32(p, x); int mask = _mm256_movemask_ps((__m256) m); reg permutation = _mm256_load_si256( (reg*) &T.permutation[mask] ); x = _mm256_permutevar8x32_epi32(x, permutation); _mm256_storeu_si256((reg*) &b[k], x); k += __builtin_popcount(mask); } return k; } ``` TODO: might be good time to explain ports and throughput Note: …and then we execute the algorithm itself. ---  ![](https://en.algorithmica.org/hpc/simd/img/filter.svg =500x) 3-7x faster than the scalar one Much faster on AVX-512: it has a dedicated `compress` instruction  Note: Мы можем использовать это в квиксорте --- ### Back to parsing Parse one integer first (`atoi`) - Find how many digits - Null the rest - Sum it up Note: Сейчас мы из этих строительных блоков соберем парсер для интов. --- ```cpp reg convert(reg x) { // merge into 8 2-digit numbers, one byte each const reg mul10 = _mm_set1_epi16((1 << 8) + 10); // (1, 10, 1, 10...) x = _mm_maddubs_epi16(x, mul10); // merge into 4 4-digit numbers, two bytes each const reg mul100 = _mm_set1_epi32((1 << 16) + 100); // (1, 100, 1, 100...) x = _mm_madd_epi16(x, mul100); // pack numbers tighter (second half will be a copy of first; we don't need it anyway) x = _mm_packus_epi32(x, x); // merge into 2 8-digit numbers, 4 bytes each (and their copy on the right half) const reg mul10000 = _mm_set1_epi32((1 << 16) + 10000); // (1, 10000, 1, 10000...) x = _mm_madd_epi16(x, mul10000); return x; } ``` --- ```cpp int k = 0; // how many numbers we've parsed int pos = 0; // how many bytes we've parsed while (pos + 15 < fsize) { reg x = _mm_loadu_si128( (reg*) &input[pos] ); const reg zero = _mm_set1_epi8('0'); reg mask = _mm_cmplt_epi8(x, zero); int m = _mm_movemask_epi8(mask); // zeros correspond to digits int d = __builtin_ffs(m); // relative position of first separator (1-indexed) buf_pos += d; // convert ASCII chars to 2x8 numbers 0..9, one byte each x = _mm_subs_epu8(x, zero); x = _mm_slli_epi64(x, 8 * (9 - d)); // shift right by (9 - d) bytes a[k++] = convert(x); } ``` --- ![](https://hackmd.io/_uploads/H1vc9b292.png =600x) ~1.9 cycles per byte (if properly optimized) Performance is proportional to the number count, not file size (performs terribly when integers are small) --- - When we apply madubbs to a chunk of input, we get useful partial sums - We can multiply these partial sums by required powers of 10 and add the results to needed array cells using (a lot of) lookup tables and permutations - Complicated idea to implement, but this way we truly process 32 bytes at once --- ![](https://hackmd.io/_uploads/BkU5qZnqn.png =600x) ~1.55 cycles per byte (now independent of number sizes) --- - Scalar solution was considerably simpler, but we can't vectorize it directly because of dependencies between cycles - …unless we "transpose" the input array, so that we fetch 8 numbers at a time from different "partitions" ```cpp for (int i = 0; i < fsize; i++) { char c = input[i]; if (c >= '0') x = 10 * x + (c - '0'); else { a[k++] = x; x = 0; } } ``` --- ```cpp int a[16], b[16]; reg r1 = _mm_load_si128((reg*) &a[0]); reg r2 = _mm_load_si128((reg*) &a[4]); reg r3 = _mm_load_si128((reg*) &a[8]); reg r4 = _mm_load_si128((reg*) &a[12]); 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ``` ```cpp _mm_unpacklo_epi32(r1, r2) = 0 4 1 5 _mm_unpackhi_epi32(r1, r2) = 2 6 3 7 ``` --- ```cpp reg t1 = _mm_unpacklo_epi32(r1, r3); reg t2 = _mm_unpackhi_epi32(r1, r3); reg t3 = _mm_unpacklo_epi32(r2, r4); reg t4 = _mm_unpackhi_epi32(r2, r4); 0 8 1 9 2 10 3 11 4 12 5 13 6 14 7 15 ``` --- ```cpp reg q1 = _mm_unpacklo_epi32(t1, t3); reg q2 = _mm_unpackhi_epi32(t1, t3); reg q3 = _mm_unpacklo_epi32(t2, t4); reg q4 = _mm_unpackhi_epi32(t2, t4); 0 4 8 12 1 5 9 13 2 6 10 14 3 7 11 15 ``` Generalizing, we can transpose $n x n$ matrix in $O(n log n)$ SIMD instructions --- ```cpp void transpose(char *a, reg *b) { reg r1 = _mm256_load_si256((reg*) &a[0 * P]); reg r5 = _mm256_load_si256((reg*) &a[1 * P]); reg r3 = _mm256_load_si256((reg*) &a[2 * P]); reg r7 = _mm256_load_si256((reg*) &a[3 * P]); reg r2 = _mm256_load_si256((reg*) &a[4 * P]); reg r6 = _mm256_load_si256((reg*) &a[5 * P]); reg r4 = _mm256_load_si256((reg*) &a[6 * P]); reg r8 = _mm256_load_si256((reg*) &a[7 * P]); reg d1 = _mm256_unpacklo_epi8(r1, r2); reg d2 = _mm256_unpackhi_epi8(r1, r2); reg d3 = _mm256_unpacklo_epi8(r3, r4); reg d4 = _mm256_unpackhi_epi8(r3, r4); reg d5 = _mm256_unpacklo_epi8(r5, r6); reg d6 = _mm256_unpackhi_epi8(r5, r6); reg d7 = _mm256_unpacklo_epi8(r7, r8); reg d8 = _mm256_unpackhi_epi8(r7, r8); reg q1 = _mm256_unpacklo_epi8(d1, d3); reg q2 = _mm256_unpackhi_epi8(d1, d3); reg q3 = _mm256_unpacklo_epi8(d2, d4); reg q4 = _mm256_unpackhi_epi8(d2, d4); reg q5 = _mm256_unpacklo_epi8(d5, d7); reg q6 = _mm256_unpackhi_epi8(d5, d7); reg q7 = _mm256_unpacklo_epi8(d6, d8); reg q8 = _mm256_unpackhi_epi8(d6, d8); r1 = _mm256_unpacklo_epi8(q1, q5); r2 = _mm256_unpackhi_epi8(q1, q5); r3 = _mm256_unpacklo_epi8(q2, q6); r4 = _mm256_unpackhi_epi8(q2, q6); r5 = _mm256_unpacklo_epi8(q3, q7); r6 = _mm256_unpackhi_epi8(q3, q7); r7 = _mm256_unpacklo_epi8(q4, q8); r8 = _mm256_unpackhi_epi8(q4, q8); _mm256_store_si256(b + 0, r1); _mm256_store_si256(b + 1, r2); _mm256_store_si256(b + 2, r3); _mm256_store_si256(b + 3, r4); _mm256_store_si256(b + 4, r5); _mm256_store_si256(b + 5, r6); _mm256_store_si256(b + 6, r7); _mm256_store_si256(b + 7, r8); } ``` For 8x8, we have 3x8=24 unpacklo/unpackhi (ports 1 and 2) + 8 store (port 2) We need (24 + 8) / 2 cycles to process 256 bytes → 16/256 = 1/16 = 0.0625 bytes per cycle --- ```cpp for (int i = 0; i < 8; i++) { const reg zero_char = _mm256_set1_epi32('0'); const reg zero_int = _mm256_setzero_si256(); const reg ten = _mm256_set1_epi32(10); reg c = _mm256_cvtepu8_epi32( _mm_loadl_epi64((__m128i*) &in_buffer[i]) ); reg mask = _mm256_cmpgt_epi32(zero_char, c); // 1s are separators c = _mm256_sub_epi32(c, zero_char); reg y = _mm256_add_epi32(_mm256_mullo_epi32(x, ten), c); y = _mm256_blendv_epi8(y, zero_int, mask); _mm256_store_si256((reg*) out_buffer[i], y); } ``` --- transpose + process + transpose back (4x data) + filter results 0.0625 + 0.2 + 0.25 + 0.5 ≈ 1 ![](https://hackmd.io/_uploads/SkLc9-393.png =500x) --- We can use `maddubs` to 'compress" the input stream 2x (0.0625 + 0.**3** + 0.25 + 0.5) / 2 ≈ 0.55 --- ![](https://hackmd.io/_uploads/rJL55Z35h.png =600x) 7-8 faster than scalar parsing, 50-100 faster than standard libraries --- ### Negative numbers? - For one-by-one methods, just check for sign and negate the number, costs 3 instructions per number (fetch-cmp, neg and cmov) - For transpose-based solution, we can add a separate flag - We can also always correct it after parsing the integers --- ### Floating-point numbers? 3.14159265358979323846 6.0221408e+23 Complicated, but the computation-heavy part is still parsing digits --- ### Parsing CSV ```text respondent_id,talk_id,talk_score 1,1,7 1,2,9 2,1,6 2,2,5 ... ``` --- ```text respondent_id,talk_id,talk_score,feedback 1,1,, 1,2,10,"Great talk, learned a lot" 2,1,7.5,"讲得很好，学到了很多" 2,2,1,"This ""transpose"" algorithm is too data-specific.\nAnd who even cares about parsing anyway?" ... ``` - Find all quoted strings (bitmask of chars preceded by an odd number of quotes, except escaped ones; can be done with a lookup table) - Find all commas (separators) not inside quotes - Parse every field --- ### Parsing JSON ```json { "statuses": [ { "metadata": { "result_type": "recent", "iso_language_code": "ja" }, "created_at": "Sun Aug 31 00:29:15 +0000 2014", "id": 505874924095815681, "id_str": "505874924095815681", "text": "@aym0566x \n\n名前:前田あゆみ\n第一印象:なんか怖っ！😋✨✨", "source": "<a href=\"http://twitter.com/download/iphone\" rel=\"nofollow\">Twitter for iPhone</a>", "truncated": false, "in_reply_to_status_id": null, "in_reply_to_status_id_str": null, "in_reply_to_user_id": 866260188, ... ``` --- ```cpp void parse_int() { // ... } void parse_string() { // ... } void parse_list() { // ... } void parse_object() { while (next_symbol() != ']') { if (next_symbol() = '[') { parse_list(); } else if (next_symbol() = '{') { parse_object(); } else if (next_symbol() = '"') { parse_string(); } else if (is_digit(next_symbol())) { parse_int(); } } } ``` - We want a goto-based automata, with states and transitions - Use SIMD for a "Stage 1" preprocessing and converting tokens to transitions - No need to perform any comparisons in "stage 2" - There won't be branch mispredictions if JSON structure is predictable --- ### Further reading - *Algorithms for Modern Hardware*, Chapter 10 - Geoff Langdale ([@geofflangdale](https://twitter.com/geofflangdale), Wojciech Muła ([@pshufb](https://twitter.com/pshufb)), Daniel Lemire ([@lemire]((https://twitter.com/lemire))), and me ([@sergey_slotin]((https://twitter.com/sergey_slotin))) - en.algorithmica.org/hpc/algorithms/parsing-integers - en.algorithmica.org/hpc/algorithms/printing-integers? Note: ... --- ![](https://hackmd.io/_uploads/HkM536iqn.png =700x) --- ![](https://hackmd.io/_uploads/H1d53Toc2.png =700x) --- ![](https://hackmd.io/_uploads/SJJpE6i9n.png =350x) --- ![](https://hackmd.io/_uploads/HJgEBai53.png =400x) --- ![](https://hackmd.io/_uploads/rk8JXTic3.png =300x) ![](https://hackmd.io/_uploads/r1DxXajc3.png =300x) --- ![](https://hackmd.io/_uploads/S18VMpi93.png =450x) ![](https://hackmd.io/_uploads/H1TrQpsc3.png =450x)