Использование вложенных векторов против плоской векторной оболочки, странное поведение

Проблема

В течение долгого времени у меня создалось впечатление, что использование вложенного std::vector<std::vector...> для имитации N-мерного массива вообще плохое, поскольку память не является гарантией непрерывности и может иметь недостатки в кэше. Я подумал, что лучше использовать плоский вектор и карту из нескольких измерений в 1D и наоборот. Итак, я решил проверить его (код указан в конце). Это довольно просто, я приурочил чтение/запись к вложенному 3D-вектору против моей собственной 3D-обертки 1D-вектора. Я скомпилировал код с g++ и clang++, при включенной оптимизации -O3. Для каждого запуска я изменил размеры, поэтому я могу получить довольно хорошее представление о поведении. К моему удивлению, это те результаты, которые я получил на своей машине MacBook Pro (Retina, 13 дюймов, конец 2012 года), 2,5 ГГц i5, 8 ГБ оперативной памяти, OS X 10.10.5:

g++ 5.2

dimensions       nested   flat
X   Y   Z        (ms)     (ms) 

100 100 100  ->  16       24
150 150 150  ->  58       98
200 200 200  ->  136     308
250 250 250  ->  264     746
300 300 300  ->  440    1537

clang++ (LLVM 7.0.0)

dimensions       nested   flat
X   Y   Z        (ms)     (ms) 

100 100 100  ->  16       18
150 150 150  ->  53       61
200 200 200  ->  135     137
250 250 250  ->  255     271
300 300 300  ->  423     477


Как вы можете видеть, обертка "сгладить" никогда не избивает вложенную версию. Более того, реализация g++ libstdС++ выполняется довольно плохо по сравнению с реализацией libС++, например, для 300 x 300 x 300 версия сглаживания почти в 4 раза медленнее, чем вложенная версия. libС++, похоже, имеет равную производительность.

Мои вопросы:

  • Почему не сгладить версию быстрее? Разве это не должно быть? Я что-то пропустил в тестовом коде?
  • Кроме того, почему g++ libstdС++ так плохо работает при использовании сглаженных векторов? Опять же, лучше не работать?

Код, который я использовал:

#include <chrono>
#include <cstddef>
#include <iostream>
#include <memory>
#include <random>
#include <vector>

// Thin wrapper around flatten vector
template<typename T>
class Array3D
{
    std::size_t _X, _Y, _Z;
    std::vector<T> _vec;
public:
    Array3D(std::size_t X, std::size_t Y, std::size_t Z):
        _X(X), _Y(Y), _Z(Z), _vec(_X * _Y * _Z) {}
    T& operator()(std::size_t x, std::size_t y, std::size_t z)
    {
        return _vec[z * (_X * _Y) + y * _X + x];
    }
    const T& operator()(std::size_t x, std::size_t y, std::size_t z) const
    {
        return _vec[z * (_X * _Y) + y * _X + x];
    }
};

int main(int argc, char** argv)
{
    std::random_device rd{};
    std::mt19937 rng{rd()};
    std::uniform_real_distribution<double> urd(-1, 1);

    const std::size_t X = std::stol(argv[1]);
    const std::size_t Y = std::stol(argv[2]);
    const std::size_t Z = std::stol(argv[3]);


    // Standard library nested vector
    std::vector<std::vector<std::vector<double>>>
        vec3D(X, std::vector<std::vector<double>>(Y, std::vector<double>(Z)));

    // 3D wrapper around a 1D flat vector
    Array3D<double> vec1D(X, Y, Z);

    // TIMING nested vectors
    std::cout << "Timing nested vectors...\n";
    auto start = std::chrono::steady_clock::now();
    volatile double tmp1 = 0;
    for (std::size_t x = 0 ; x < X; ++x)
    {
        for (std::size_t y = 0 ; y < Y; ++y)
        {
            for (std::size_t z = 0 ; z < Z; ++z)
            {
                vec3D[x][y][z] = urd(rng);
                tmp1 += vec3D[x][y][z];
            }
        }
    }
    std::cout << "\tSum: " << tmp1 << std::endl; // we make sure the loops are not optimized out
    auto end = std::chrono::steady_clock::now();
    std::cout << "Took: ";
    auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
    std::cout << ms << " milliseconds\n";

    // TIMING flatten vector
    std::cout << "Timing flatten vector...\n";
    start = std::chrono::steady_clock::now();
    volatile double tmp2 = 0;
    for (std::size_t x = 0 ; x < X; ++x)
    {
        for (std::size_t y = 0 ; y < Y; ++y)
        {
            for (std::size_t z = 0 ; z < Z; ++z)
            {
                vec1D(x, y, z) = urd(rng);
                tmp2 += vec1D(x, y, z);
            }
        }
    }
    std::cout << "\tSum: " << tmp2 << std::endl; // we make sure the loops are not optimized out
    end = std::chrono::steady_clock::now();
    std::cout << "Took: ";
    ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
    std::cout << ms << " milliseconds\n";
}

ИЗМЕНИТЬ

Изменение Array3D<T>::operator() возвращается к

return _vec[(x * _Y + y) * _Z + z];

в соответствии с предложением @1201ProgramAlarm действительно избавляется от "странного" поведения g++ в том смысле, что плоские и вложенные версии берут сейчас примерно одно и то же время, Однако это все еще интригует. Я думал, что вложенное будет намного хуже из-за проблем с кешем. Могу ли я просто повезти и все выделены в память?

Ответы

Ответ 1

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

Я поставил исходный код с производительностью bugfixed на Godbolt, чтобы мы могли посмотреть на asm внутреннего цикла, скомпилированного g++ 5.2, с -O3. (Яблочная вилка clang может быть похожа на clang3.7, но я просто посмотрю на версию gcc.) Там много кода из функций С++, но вы можете щелкнуть правой кнопкой мыши по исходной строке, чтобы прокрутить окна asm до код для этой строки. Кроме того, наведите указатель на строку источника, выделенную жирным шрифтом asm, который реализует эту строку, или наоборот.

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

## outer-most loop not shown

.L213:  ## middle loop (over `y`)
    test    rbp, rbp        # Z
    je      .L127           # inner loop runs zero times if Z==0
    mov     rax, QWORD PTR [rsp+80]   # MEM[(struct vector * *)&vec3D], MEM[(struct vector * *)&vec3D]
    xor     r15d, r15d        # z = 0
    mov     rax, QWORD PTR [rax+r12]  # MEM[(struct vector * *)_195], MEM[(struct vector * *)_195]
    mov     rdx, QWORD PTR [rax+rbx]  # D.103857, MEM[(double * *)_38]

## Top of inner-most loop.
.L128:
    lea     rdi, [rsp+5328]   # tmp511,   ## function arg: pointer to the RNG object, which is a local on the stack.
    lea     r14, [rdx+r15*8]  # D.103851,  ## r14 = &(vec3D[x][y][z])
    call    double std::generate_canonical<double, 53ul, std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >(std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul>&)  #
    addsd   xmm0, xmm0    # D.103853, D.103853  ## return val *= 2.0: [0.0, 2.0]
    mov     rdx, QWORD PTR [rsp+80]   # MEM[(struct vector * *)&vec3D], MEM[(struct vector * *)&vec3D]   ## redo the pointer-chasing from vec3D.data()
    mov     rdx, QWORD PTR [rdx+r12]  # MEM[(struct vector * *)_150], MEM[(struct vector * *)_150]
    subsd   xmm0, QWORD PTR .LC6[rip]     # D.103859, ## and subtract 1.0:  [-1.0, 1.0]
    mov     rdx, QWORD PTR [rdx+rbx]  # D.103857, MEM[(double * *)_27]
    movsd   QWORD PTR [r14], xmm0 # *_155, D.103859        # store into vec3D[x][y][z]
    movsd   xmm0, QWORD PTR [rsp+64]      # D.103853, tmp1  # reload volatile tmp1
    addsd   xmm0, QWORD PTR [rdx+r15*8]   # D.103853, *_62  # add the value just stored into the array (r14 = rdx+r15*8 because nothing else modifies the pointers in the outer vectors)
    add     r15, 1    # z,
    cmp     rbp, r15  # Z, z
    movsd   QWORD PTR [rsp+64], xmm0      # tmp1, D.103853  # spill tmp1
    jne     .L128     #,
 #End of inner-most loop

.L127:  ## middle-loop
    add     r13, 1    # y,
    add     rbx, 24           # sizeof(std::vector<> == 24) == the size of 3 pointers.
    cmp     QWORD PTR [rsp+8], r13    # %sfp, y
    jne     .L213     #,

 ## outer loop not shown.

И для плоского контура:

 ## outer not shown.
.L214:
    test    rbp, rbp        # Z
    je      .L135       #,
    mov     rax, QWORD PTR [rsp+280]  # D.103849, vec1D._Y
    mov     rdi, QWORD PTR [rsp+288]  # D.103849, vec1D._Z
    xor     r15d, r15d        # z
    mov     rsi, QWORD PTR [rsp+296]  # D.103857, MEM[(double * *)&vec1D + 24B]

.L136:  ## inner-most loop
    imul    rax, r12        # D.103849, x
    lea     rax, [rax+rbx]    # D.103849,
    imul    rax, rdi        # D.103849, D.103849
    lea     rdi, [rsp+5328]   # tmp520,
    add     rax, r15  # D.103849, z
    lea     r14, [rsi+rax*8]  # D.103851,       # &vec1D(x,y,z)
    call    double std::generate_canonical<double, 53ul, std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >(std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul>&)  #
    mov     rax, QWORD PTR [rsp+280]  # D.103849, vec1D._Y
    addsd   xmm0, xmm0    # D.103853, D.103853
    mov     rdi, QWORD PTR [rsp+288]  # D.103849, vec1D._Z
    mov     rsi, QWORD PTR [rsp+296]  # D.103857, MEM[(double * *)&vec1D + 24B]
    mov     rdx, rax  # D.103849, D.103849
    imul    rdx, r12        # D.103849, x       # redo address calculation a 2nd time per iteration
    subsd   xmm0, QWORD PTR .LC6[rip]     # D.103859,
    add     rdx, rbx  # D.103849, y
    imul    rdx, rdi        # D.103849, D.103849
    movsd   QWORD PTR [r14], xmm0 # MEM[(double &)_181], D.103859  # store into the address calculated earlier
    movsd   xmm0, QWORD PTR [rsp+72]      # D.103853, tmp2
    add     rdx, r15  # tmp374, z
    add     r15, 1    # z,
    addsd   xmm0, QWORD PTR [rsi+rdx*8]   # D.103853, MEM[(double &)_170]   # tmp2 += vec1D(x,y,z).  rsi+rdx*8 == r14, so this is a reload of the store this iteration.
    cmp     rbp, r15  # Z, z
    movsd   QWORD PTR [rsp+72], xmm0      # tmp2, D.103853
    jne     .L136     #,

.L135:  ## middle loop: increment y
    add     rbx, 1    # y,
    cmp     r13, rbx  # Y, y
    jne     .L214     #,

 ## outer loop not shown.

У вашего MacBook Pro (Late 2012) есть Intel IvyBridge, поэтому я использую номера для этой микроархитектуры от Таблицы инструкций Agner Fog и руководство микроархива. Вещи должны быть в основном одинаковыми на других процессорах Intel/AMD.

Единственный 2,5 ГГц мобильный IvB i5 - i5-3210M, поэтому ваш процессор имеет 3MiB кеша L3. Это означает, что даже ваш самый маленький тестовый пример (100 ^ 3 * 8B на double ~ = 7.63MiB) больше, чем ваш кеш последнего уровня, поэтому ни один из ваших тестовых случаев не подходит в кеше вообще. Вероятно, это хорошо, потому что вы выделяете и дефолт инициализируете как вложенные, так и плоские перед тем, как тестировать их. Тем не менее, вы тестируете в том же порядке, который вы назначили, поэтому, если вложенный массив по-прежнему остается кешем после обнуления плоского массива, плоский массив может оставаться горячим в кэше L3 после цикла синхронизации по вложенному массиву.

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


Вы делаете здесь несколько вещей, которые являются сверхъестественными и делают это настолько медленным, что выполнение вне порядка может скрыть дополнительную латентность изменения y, даже если ваши внутренние векторы z не являются полностью смежными.

  • Вы запускаете медленный PRNG внутри синхронизированного цикла. std::uniform_real_distribution<double> urd(-1, 1); - дополнительные накладные расходы поверх std::mt19937 rng{rd()};, которые уже медленны по сравнению с задержкой FP-добавления (3 цикла), или по сравнению с пропускной способностью кэша L1D по 2 за цикл. Все это дополнительное время, выполняемое PRNG, приводит к выполнению команды вне очереди, чтобы запустить инструкции по индексированию массива, чтобы конечный адрес был готов к моменту, когда данные были. Если у вас много промахов в кеше, вы в основном просто измеряете скорость PRNG, потому что она дает результаты намного медленнее, чем 1 за такт.

    g++ 5.2 не полностью встраивает код urd(rng), а соглашение о вызове System 86 x86-64 не имеет регистров XMM с сохранением номера. Поэтому tmp1/tmp2 необходимо проливать/перезагружать для каждого элемента, даже если они не были volatile.

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

    clang (с libС++) полностью встраивает PRNG, поэтому переход к следующему z - это просто add reg, 8, чтобы увеличить указатель как в плоских, так и вложенных версиях. Вы можете получить такое же поведение из gcc, получив итератор за пределами внутреннего цикла, или получив ссылку на внутренний вектор, вместо повторного использования operator[] и надеясь, что компилятор поднимет его для вас.

    Intel/AMD FP add/sub/mul пропускная способность/задержка не зависят от данных, за исключением денормалов. (x87 также замедляется для NaN и, возможно, бесконечности, но SSE этого не делает. 64-разрядный код использует SSE даже для скалярного float/double.) Таким образом, вы могли бы просто инициализировать ваш массив нулями или PRNG outisde цикл синхронизации. (Или оставил их обнуленными, так как конструктор vector<double> делает это для вас, и на самом деле он берет дополнительный код, чтобы он не был в тех случаях, когда вы собираетесь писать что-то еще.) Производительность деления и sqrt зависит от данных на некоторые процессоры и намного медленнее, чем add/sub/mul.

  • Вы пишете каждый элемент прямо перед его чтением, внутри внутреннего цикла. В источнике это выглядит как store/reload. То, что gcc на самом деле делает, к сожалению, но clang с libС++ (который встраивает PRNG), преобразует тело цикла:

           // original
           vec3D[x][y][z] = urd(rng);
           tmp1 += vec3D[x][y][z];
    
           // what clang asm really does
           double xmm7 = urd(rng);  
           vec3D[x][y][z] = xmm7;
           tmp1 += xmm7;
    

    В clang asm:

                                   # do { ...
    
    addsd   xmm7, xmm4             # last instruction of the PRNG
    movsd   qword ptr [r8], xmm7   # store it into the Z vector
    addsd   xmm7, qword ptr [rsp + 88]
    add     r8, 8                  # pointer-increment to walk along the Z vector
    dec     r13                    # i--
    movsd   qword ptr [rsp + 88], xmm7
    jne     .LBB0_74               # }while(i != 0);
    

    Это разрешено, потому что vec3D не volatile или atomic<>, поэтому было бы undefined поведение для любого другого потока записывать эту память одновременно. Это означает, что он может оптимизировать хранение/перезагрузку объектов в памяти только в хранилище (и просто использовать значение, которое он хранит, без перезагрузки). Или полностью оптимизируйте хранилище, если он может доказать, что это мертвый магазин (хранилище, которое ничего не может прочитать, например, неиспользуемой переменной static).

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

    Однако даже фактическое сохранение/перезагрузка в asm все еще менее чувствительно к ошибкам кэша, чем простая загрузка!

    Store- > load forwarding по-прежнему работает, даже если хранилище пропущено в кеше. Таким образом, кэш-промах в векторе Z напрямую не задерживает критический путь. Это только замедляет вас, если выполнение вне порядка не может скрыть латентность промаха в кеше. (Магазин может уйти в отставку, как только данные будут записаны в буфер-хранилище (и все предыдущие инструкции удалены). Я не уверен, что загрузка может уйти в отставку до того, как кэш-строка даже сделает ее L1D, если она его данные из пересылки хранилища.Она может быть в состоянии, поскольку x86 позволяет переупорядочивать StoreLoad (магазины могут стать глобально видимыми после загрузки). В этом случае хранилище/перезагрузка добавляет только 6 циклов задержки для результата PRNG ( от критического пути от одного состояния PRNG до следующего состояния PRNG). И это только пропускная способность узкого места, если он кэширует настолько, что буфер хранилища заполняется и предотвращает выполнение новых операций хранения, что, в свою очередь, в конечном итоге предотвращает новые ошибки выдается в ядро ​​внешнего порядка, когда станция резервирования или ROB заполняется не выполненными или не удаленными (соответственно) uops.

    При обратном индексировании (исходная версия плоского кода), вероятно, основным узким местом было разбросанные магазины. IDK, почему clang сделал намного лучше, чем gcc. Может быть, clang удалось инвертировать цикл и перемещать память в последовательном порядке. (Поскольку он полностью ввел PRNG, не было никаких вызовов функций, которые требовали бы, чтобы память соответствовала программному порядку.)

    Перемещение каждого вектора Z в порядке означает, что промахи кэша относительно далеко друг от друга (даже если каждый вектор Z не соприкасается с предыдущим), что дает много времени для выполнения хранилищами. Или даже если перенаправленная с хранения загрузка не может фактически уйти в отставку, пока кэш L1D фактически не владеет линией кэша (в измененном состоянии протокола MESI), спекулятивное выполнение имеет правильные данные и не нужно ждать латентности кэш-промах. Окно инструкции вне порядка, вероятно, достаточно велико, чтобы критический путь, вероятно, был остановлен, прежде чем загрузка может выйти на пенсию. (Кэш-пропускные нагрузки обычно очень плохие, потому что зависимые инструкции не могут выполняться без данных для их работы. Таким образом, они намного легче создают пузырьки в конвейере. С полным отсутствием кэша от DRAM с задержкой более 300 циклов, а окно вне порядка - 168 мкп на IvB, оно не может скрыть всю задержку для выполнения кода с точностью до 1 мкп (приблизительно 1 команда) за такт.) Для чистых магазинов, окно заказа выходит за пределы размера ROB, потому что им не нужно фиксировать L1D для выхода на пенсию. Фактически, они не могут совершать действия до тех пор, пока они не уйдут на пенсию, потому что тот момент, когда они, как известно, не являются спекулятивными. (Так что сделать их глобально видимыми раньше, чем это предотвратит откатывание при обнаружении исключения или неправильной спекуляции.)

    У меня нет libc++, установленного на моем рабочем столе, поэтому я не могу сравнить эту версию с g++. С g++ 5.4 я нахожу Nested: 225 миллисекунд и Flat: 239 миллисекунд. Я подозреваю, что дополнительные умножения на массивы являются проблемой и конкурируют с инструкциями ALU, которые использует PRNG. Напротив, вложенная версия, переработавшая кучу слежения за указателем, попадающая в кеш L1D, может происходить параллельно. Мой рабочий стол - Skylake i7-6700k на 4.4 ГГц. SKL имеет размер ROB (ReOrder Buffer) 224 uops и RS из 97 uops, поэтому окно вне порядка очень велико. Он также имеет FP-добавление латентности в 4 цикла (в отличие от предыдущих uarches, где это было 3).

  • volatile double tmp1 = 0; Аккумулятор volatile, который заставляет компилятор хранить/перезагружать его каждую итерацию внутреннего цикла. Полная латентность цепи зависимостей, связанной с циклом, в внутренний цикл равен 9 циклам: 3 для addsd и 6 для хранения-пересылки из movsd в хранилище до перезагрузки movsd. (clang сбрасывает перезагрузку в операнд памяти с помощью addsd xmm7, qword ptr [rsp + 88], но с той же разницей. ([rsp+88] находится в стеке, где хранятся переменные с автоматическим хранением, если они должны быть выписаны из регистров.)

    Как отмечалось выше, вызов non-inline-функции для gcc также приведет к разрыву/перезагрузке в соглашении на вызов x86-64 System V (используемом всеми, кроме Windows). Но умный компилятор мог бы сделать, например, 4 вызова PRNG, а затем 4 массива массивов. (Если бы вы использовали итератор, чтобы убедиться, что gcc знал, что векторы, содержащие другие векторы, не изменяются.)

    Использование -ffast-math позволило бы авто-векторизовать компилятор (если бы не PRNG и volatile). Это позволило бы вам быстро пройти через массивы, чтобы отсутствие реальной локальности между различными Z-векторами могло быть реальной проблемой. Это также позволит компиляторам развернуть несколько аккумуляторов, чтобы скрыть задержку FP-добавления. например они могут (и clang) сделать asm эквивалентом:

    float t0=0, t1=0, t2=0, t3=0;
    for () {
       t0 += a[i + 0];
       t1 += a[i + 1];
       t2 += a[i + 2];
       t3 += a[i + 3];
    }
    t0 = (t0 + t1) + (t2 + t3);
    

    У этого есть 4 отдельные цепи зависимостей, таким образом он может держать 4 FP добавляет в полете. Поскольку IvB имеет задержку в 3 цикла, по одной на единицу времени для addsd, нам нужно только сохранить 4 в полете, чтобы насытить его пропускную способность. (Skylake имеет 4c латентность, 2 за пропускную способность тактовой частоты, что и mul или FMA, поэтому вам нужно 8 аккумуляторов, чтобы избежать узких мест с задержкой. На самом деле, еще лучше. как выяснил этот вопрос, Хасуэлл справился с еще большим количеством аккумуляторов при приближении к максимальной нагрузке.)

    Что-то вроде этого было бы гораздо лучшим испытанием того, насколько эффективным является цикл над Array3D. Если вы хотите полностью отключить цикл от целиком, используйте результат. Проверьте свой микрообъект, чтобы убедиться, что увеличение размера проблемы масштабирует время; если нет, то что-то оптимизировалось, или вы не тестируете то, что, по вашему мнению, испытываете. Не делайте временную временную привязку volatile!!

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


Могу ли я просто повезти и все выделены в память?

Да, это, вероятно, происходит для многих небольших распределений, выполненных в порядке, когда вы ничего не выделяли и не освобождали, прежде чем делать это. Если они были достаточно большими (обычно одна страница 4kiB или больше), glibc malloc переключится на использование mmap(MAP_ANONYMOUS), а затем ядро ​​выберет рандомизированные виртуальные адреса (ASLR). Таким образом, с большим Z, вы можете ожидать, что местность будет ухудшаться. Но, с другой стороны, большие Z-векторы означают, что вы тратите больше времени на один непрерывный вектор, поэтому пропустить кеш при изменении yx) становится относительно менее важным.

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

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


Различные компиляторы/библиотеки могут сильно повлиять на этот странный тест. В моей системе (Arch Linux, i7-6700k Skylake с максимальным турбонаддувом 4,4 ГГц) лучшее из 4 работает в 300 300 300 для g++ 5.4 -O3:

Timing nested vectors...
        Sum: 579.78
Took: 225 milliseconds
Timing flatten vector...
        Sum: 579.78
Took: 239 milliseconds

 Performance counter stats for './array3D-gcc54 300 300 300':

    532.066374      task-clock (msec)         #    1.000 CPUs utilized          
             2      context-switches          #    0.004 K/sec                  
             0      cpu-migrations            #    0.000 K/sec                  
        54,523      page-faults               #    0.102 M/sec                  
 2,330,334,633      cycles                    #    4.380 GHz                    
 7,162,855,480      instructions              #    3.07  insn per cycle         
   632,509,527      branches                  # 1188.779 M/sec                  
       756,486      branch-misses             #    0.12% of all branches        

   0.532233632 seconds time elapsed

против. g++ 7.1 -O3 (который, по-видимому, решил разветкиться на то, что g++ 5.4 не сделал)

Timing nested vectors...
        Sum: 932.159
Took: 363 milliseconds
Timing flatten vector...
        Sum: 932.159
Took: 378 milliseconds

 Performance counter stats for './array3D-gcc71 300 300 300':

    810.911200      task-clock (msec)         #    1.000 CPUs utilized          
             0      context-switches          #    0.000 K/sec                  
             0      cpu-migrations            #    0.000 K/sec                  
        54,523      page-faults               #    0.067 M/sec                  
 3,546,467,563      cycles                    #    4.373 GHz                    
 7,107,511,057      instructions              #    2.00  insn per cycle         
   794,124,850      branches                  #  979.299 M/sec                  
    55,074,134      branch-misses             #    6.94% of all branches        

   0.811067686 seconds time elapsed

против. clang4.0 -O3 (с gcc libstdС++, а не libС++)

perf stat ./array3D-clang40-libstdc++ 300 300 300
Timing nested vectors...
        Sum: -349.786
Took: 1657 milliseconds
Timing flatten vector...
        Sum: -349.786
Took: 1631 milliseconds

 Performance counter stats for './array3D-clang40-libstdc++ 300 300 300':

   3358.297093      task-clock (msec)         #    1.000 CPUs utilized          
             9      context-switches          #    0.003 K/sec                  
             0      cpu-migrations            #    0.000 K/sec                  
        54,521      page-faults               #    0.016 M/sec                  
14,679,919,916      cycles                    #    4.371 GHz                    
12,917,363,173      instructions              #    0.88  insn per cycle         
 1,658,618,144      branches                  #  493.887 M/sec                  
       916,195      branch-misses             #    0.06% of all branches        

   3.358518335 seconds time elapsed

Я не копал в том, что clang сделал не так, или попробуйте с -ffast-math и/или -march=native. (Тем не менее, вы ничего не сделаете, если не удалите volatile.)

perf stat -d не отображает больше промахов в кеше (L1 или последний уровень) для clang, чем gcc. Но он показывает, что clang делает больше, чем вдвое больше нагрузок L1D.

Я попытался использовать неквадратный массив. Это почти точно то же самое время, когда общее количество элементов подсчитывается одинаково, но изменяя окончательный размер до 5 или 6.

Даже незначительное изменение в C помогает и "сглаживается" быстрее, чем вложенный с gcc (от 240 мс до 220 мс для 300 ^ 3, но едва ли имеет значение для вложенных.):

 // vec1D(x, y, z) = urd(rng);
    double res = urd(rng);
    vec1D(x, y, z) = res;   // indexing calculation only done once, after the function call
    tmp2 += vec1D(x, y, z);
    // using iterators would still avoid redoing it at all.

Ответ 2

Это из-за того, как вы заказываете свои индексы в 3D-классе. Поскольку ваш внутренний цикл меняет z, это самая большая часть вашего индекса, поэтому вы получаете много промахов в кеше. Переупорядочивайте индексирование до

_vec[(x * _Y + y) * _Z + z]

и вы должны увидеть лучшую производительность.

Ответ 3

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

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

В основном есть две вещи, которые делают кеширование особенно эффективным:

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

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

Чтобы оценить влияние эффектов косвенности и кеша на эту проблему, предположим, что у нас есть X = Y = Z = 1024

Судя по этому вопросу, одна строка кеша (L1, L2 или L3) имеет длину 64 байта, что означает 8 двойных значений. Предположим, что кеш L1 имеет 32 кбайта (4096 удвоений), кэш L2 имеет 256 кбайт (32 тыс. Удваивается), а кэш L3 - 8 МБ (1 М удваивается).

Это означает, что - если кеш заполнен никакими другими данными (это смелое предположение, я знаю) - в сплющенном случае только каждое четвертое значение y приводит к пропуску кеша L1 (латентность кэша L2 вероятно, составляет около 10-20 циклов), только каждое 32-е значение y приводит к пропуску кэша L2 (латентность кэша L3 - это некоторое значение, меньшее 100 циклов), и только в случае промаха в кеше L3 нам действительно нужно получить доступ основная память. Я не хочу раскрывать здесь весь расчет, так как учет всей иерархии кэша делает его немного сложнее, но позвольте сказать, что почти все обращения к памяти могут быть кэшированы в сплющенном случае.

В исходной формулировке этого вопроса сглаженный индекс вычислялся по-разному (z * (_X * _Y) + y * _X + x), увеличение значения, которое изменяется в самой внутренней петле (z), всегда означает скачок _X * _Y * 64 bit, что приводит к гораздо больше нелокальной памяти, что увеличило количество ошибок кэша на большие суммы.

Во вложенном случае ответ сильно зависит от значения Z:

  • Если Z довольно большой, большинство обращений к памяти являются смежными, поскольку они относятся к элементам одного вектора во вложенном vector<vector<vector>>>, которые выкладываются смежно. Только когда значение y или x увеличивается, нам нужно фактически использовать косвенную функцию для извлечения указателя начала следующего самого внутреннего вектора.
  • Если Z довольно мал, нам нужно часто менять наш доступ к памяти "базовый указатель", что может привести к промахам в кеше, если области хранения самых внутренних векторов помещаются довольно случайным образом в память. Однако, если они более или менее смежны, мы можем наблюдать минимальные различия в производительности кэша.

Поскольку возник вопрос о выходе сборки, позвольте мне дать краткий обзор:

Если вы сравниваете сборку вложенного и сплющенного массива, вы заметите много сходства: существует три эквивалентных вложенных цикла и переменные счета x, y и z хранятся в регистрах. Единственное реальное отличие - кроме того, что вложенная версия использует два счетчика для каждого внешнего индекса, чтобы избежать умножения на 24 при каждом вычислении адресов, а сглаженная версия делает то же самое для самого внутреннего цикла и умножения на 8 - может быть найдена в y, где вместо того, чтобы просто увеличивать y и вычислять сплющенный индекс, нам нужно сделать три взаимозависимые загрузки памяти, чтобы определить базовый указатель для нашего внутреннего цикла:

mov rax, QWORD PTR [rdi]
mov rax, QWORD PTR [rax+r11]
mov rax, QWORD PTR [rax+r10]

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

Ответ 4

(На самом деле это не отвечает на вопрос. Я думаю, что сначала прочитал его назад, полагая, что ОП только что нашел то, что я ожидал, что вложенные векторы медленнее, чем плоские.)


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

Линия кэша - только 64B. Это 8 double s. Локальность на уровне страницы имеет значение из-за ограниченных записей TLB, а предварительная выборка требует последовательного доступа, но вы все равно получите (достаточно близко) вложенные векторы, которые распределяются сразу с большинством реализаций malloc. (Это тривиальный микрообъект, который ничего не делает, прежде чем выделять его vector s. В реальной программе, которая выделяет и освобождает некоторую память, прежде чем делать много небольших распределений, некоторые из них могут быть разбросаны по всему.)


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

Ссылка/указатель на std::vector просто указывает на блок фиксированного размера, который содержит текущий размер, выделенное пространство и указатель на буфер. IDK, если какие-либо реализации помещают буфер сразу после данных управления как часть одного и того же блока malloced, но, вероятно, это невозможно, потому что sizeof(std::vector<int>) должен быть постоянным, чтобы вы могли иметь вектор векторов. Посмотрите asm on godbolt: функция, которая просто возвращает v[10], принимает один груз с массивом arg, но две нагрузки с std::vector arg.

В реализации вложенного вектора загрузка v[x][y][z] требует 4 шагов (при условии, что указатель или ссылка на v уже находится в регистре).

  • load v.buffer_pointer или v.bp или то, что его вызывает реализация. (Указатель на массив std::vector<std::vector<double>>)
  • load v.bp[x].bp (указатель на массив std::vector<double>)
  • load v.bp[x].bp[y].bp (указатель на массив double)
  • load v.bp[x].bp[y].bp[z] (double мы хотим)

Собственный трехмерный массив, смоделированный с помощью одного std::vector, просто выполняет:

  • load v.bp (указатель на массив double)
  • load v.bp[(x*w + y)*h + z] (double мы хотим)

Множественный доступ к одному и тому же симулятору 3D-массива с разными x и y требует вычисления нового индекса, но v.bp останется в регистре. Итак, вместо 3 промахов кэша, мы получаем только один.

Перемещение 3D-массива по порядку скрывает штраф за реализацию вложенного вектора, потому что существует цикл над всеми значениями в самом внутреннем векторе, скрывающим накладные расходы при изменении x и y. Здесь предваряется предварительная выборка смежных указателей во внешних векторах, а Z достаточно мала в тестировании, что цикл над одним внутренним вектором не будет вызывать указатель для следующего значения y.


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

Ответ 5

Могу ли я просто повезти и все выделены в память?

Вероятно, да. Я немного изменил ваш образец, поэтому у нас есть эталон, который больше концентрируется на различиях между двумя подходами:
  • заполнение массива выполняется в отдельном проходе, поэтому исключается случайная скорость генератора
  • удален volatile
  • исправлена ​​небольшая ошибка (tmp1 была напечатана вместо tmp2)
  • добавлена ​​часть #if 1, которая может использоваться для рандомизации размещения vec3D в памяти. Как оказалось, это имеет огромное значение для моей машины.

Без рандомизации (я использовал параметры: 300 300 300):

Timing nested vectors...
    Sum: -131835
Took: 2122 milliseconds
Timing flatten vector...
    Sum: -131835
Took: 2085 milliseconds

Итак, есть небольшая разница, сглаженная версия немного быстрее. (Я провел несколько тестов и поставил здесь минимальное время).

При рандомизации:

Timing nested vectors...
    Sum: -117685
Took: 3014 milliseconds
Timing flatten vector...
    Sum: -117685
Took: 2085 milliseconds

Таким образом, эффект кэша можно увидеть здесь: вложенная версия на 50% медленнее. Это связано с тем, что CPU не может предсказать, какая область памяти будет использоваться, поэтому его префикс кэша неэффективен.


Здесь измененный код:

#include <chrono>
#include <cstddef>
#include <iostream>
#include <memory>
#include <random>
#include <vector>

template<typename T>
class Array3D
{
    std::size_t _X, _Y, _Z;
    std::vector<T> _vec;
    public:
    Array3D(std::size_t X, std::size_t Y, std::size_t Z):
        _X(X), _Y(Y), _Z(Z), _vec(_X * _Y * _Z) {}
    T& operator()(std::size_t x, std::size_t y, std::size_t z)
    {
        return _vec[(x * _Y + y) * _Z + z];

    }
    const T& operator()(std::size_t x, std::size_t y, std::size_t z) const
    {
        return _vec[(x * _Y + y) * _Z + z];
    }
};

double nested(std::vector<std::vector<std::vector<double>>> &vec3D, std::size_t X, std::size_t Y, std::size_t Z) {
    double tmp1 = 0;
    for (int iter=0; iter<100; iter++)
        for (std::size_t x = 0 ; x < X; ++x)
        {
            for (std::size_t y = 0 ; y < Y; ++y)
            {
                for (std::size_t z = 0 ; z < Z; ++z)
                {
                    tmp1 += vec3D[x][y][z];
                }
            }
        }
    return tmp1;
}

double flatten(Array3D<double> &vec1D, std::size_t X, std::size_t Y, std::size_t Z) {
    double tmp2 = 0;
    for (int iter=0; iter<100; iter++)
        for (std::size_t x = 0 ; x < X; ++x)
        {
            for (std::size_t y = 0 ; y < Y; ++y)
            {
                for (std::size_t z = 0 ; z < Z; ++z)
                {
                    tmp2 += vec1D(x, y, z);
                }
            }
        }
    return tmp2;
}

int main(int argc, char** argv)
{
    std::random_device rd{};
    std::mt19937 rng{rd()};
    std::uniform_real_distribution<double> urd(-1, 1);

    const std::size_t X = std::stol(argv[1]);
    const std::size_t Y = std::stol(argv[2]);
    const std::size_t Z = std::stol(argv[3]);

    std::vector<std::vector<std::vector<double>>>
        vec3D(X, std::vector<std::vector<double>>(Y, std::vector<double>(Z)));

#if 1
    for (std::size_t i = 0 ; i < X*Y; i++)
    {
        std::size_t xa = rand()%X;
        std::size_t ya = rand()%Y;
        std::size_t xb = rand()%X;
        std::size_t yb = rand()%Y;
        std::swap(vec3D[xa][ya], vec3D[xb][yb]);
    }
#endif

    // 3D wrapper around a 1D flat vector
    Array3D<double> vec1D(X, Y, Z);

    for (std::size_t x = 0 ; x < X; ++x)
    {
        for (std::size_t y = 0 ; y < Y; ++y)
        {
            for (std::size_t z = 0 ; z < Z; ++z)
            {
                vec3D[x][y][z] = vec1D(x, y, z) = urd(rng);
            }
        }
    }

    std::cout << "Timing nested vectors...\n";
    auto start = std::chrono::steady_clock::now();
    double tmp1 = nested(vec3D, X, Y, Z);
    auto end = std::chrono::steady_clock::now();
    std::cout << "\tSum: " << tmp1 << std::endl; // we make sure the loops are not optimized out
    std::cout << "Took: ";
    auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
    std::cout << ms << " milliseconds\n";

    std::cout << "Timing flatten vector...\n";
    start = std::chrono::steady_clock::now();
    double tmp2 = flatten(vec1D, X, Y, Z);
    end = std::chrono::steady_clock::now();
    std::cout << "\tSum: " << tmp2 << std::endl; // we make sure the loops are not optimized out
    std::cout << "Took: ";
    ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
    std::cout << ms << " milliseconds\n";
}