Как я могу ускорить этот код (MWE!), Например, используя ограничение

Question

Можно ли как-нибудь ускорить эту функцию:

void task(int I, int J, int K, int *L, int **ids, double *bar){ 
    double *foo[K];
    for (int k=0;k<K;k++)
        foo[k] = new double[I*L[k]];
        // I am filling these arrays somehow
        // This is not a bottleneck, hence omitted here        
    for (int i=0;i<I;i++)
        for (int j=0;j<J;j++){
            double tmp = 1.;
            for (int k=0;k<K;k++)
                tmp *= foo[k][i*L[k]+ids[j][k]]; //ids[j][k]<L[k]
            bar[i*J+j] = tmp;
        }
}

Типичные значения: I = 100,000, J = 10,000, K=3, L=[50,20,60].

Я читал, что ключевое слово / расширение __restrict__ может помочь, но я не уверен, как его применить здесь. Например, пытаясь поместить это в определение foo[k] = new double[...], я получаю error: '__restrict_ qualifiers cannot be applied to double. Кроме того, я не знаю, должен ли я / как я мог бы объявить ids и ids[j], 1<= j<= J как ограниченные.

Как примечание, в моем реальном коде я выполняю такие задачи параллельно в таком количестве потоков, в котором у моего ЦП есть ядра.

Я пишу в основном C-совместимый C ++, поэтому приветствуются решения на обоих языках.

Answer 1

restrict не имеет значения в этом случае.

Ваш алгоритм является мусором и не позволяет использовать длинные векторные операции (поэтому микрооптимизации здесь вообще не помогут).

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

Сначала пересмотрите алгоритм - преждевременная оптимизация не поможет, если алгоритм крайне неэффективен

Редактировать

После комментария ОП я просто хочу показать ему, какая разница между «наивным» и более эффективным (менее наивным, но труднее понять)

Рассмотрим четность 32-битного значения без знака. Наивный подход:

int very_naive_parity(const uint32_t val)
{
    unsigned parity = 0;

    for(unsigned bit = 0; bit < 32; bit++)
    {
        if(val & (1U << bit))
        {
            parity = !parity;
        }
    }
    return parity;
}

Это очень легко писать и понимать, но это крайне неэффективно. Для расчета этого соотношения будет выполнено не менее 288 инструкций.

эффективнее:

int parity(const uint32_t val)
{
    uint32_t tmp = val;

    tmp ^= tmp >> 16;
    tmp ^= tmp >> 8;
    tmp ^= tmp >> 4;
    return (0b110100110010110 >> (tmp & 0x0f)) & 1;
}

будет выполнено в 9 инструкциях (оба вычисления без функциональных прологов и эпилогов) Это сложнее понять? - Определенно да. Но, как я уже писал, эффективность, как правило, означает, что людям не так легко.

Answer 2

https://en.cppreference.com/w/c/language/restrict утверждает, что вы можете объявить массив из restrict указателей для удвоения, как в C99 / C11:

typedef double *array_t[10];
restrict array_t foo;        // the type of a is double *restrict[10]

Но только gcc принимает это. Я думаю, что это GCC-изм, не действительный ISO C11. (GCC также принимает
array_t restrict foo_r; но другие компиляторы этого тоже не принимают.)

ICC предупреждает "restrict" is not allowed, Clang отклоняет его с помощью

<source>:16:5: error: restrict requires a pointer or reference ('array_t' (aka 'double *[10]') is invalid)
    restrict array_t foo_r;
    ^

MSVC отклоняет его с error C2219: syntax error: type qualifier must be after '*'

Мы получаем, по сути, то же поведение в C ++ от этих компиляторов с __restrict, которое они принимают как расширение C ++ с той же семантикой, что и C99 restrict.

---

В качестве обходного пути вы можете вместо этого использовать foo вместо квалифицированного временного указателя вместо f[k][stuff]. Я думаю, это обещает, что память, на которую вы ссылаетесь через fk, не является к той же памяти, к которой вы обращаетесь через любые другие указатели в блоке, где объявлено fk.

double *__restrict fk = foo[k];
tmp *= fk[ stuff ];

Я не знаю, как пообещать компилятору, что ни один из f[0..K-1] не указывает псевдонимы друг друга. Я не думаю, что это делает это.

---

Вам не нужно __ограничение здесь.

Я добавил __restrict ко всем объявлениям указателей, например, int *__restrict *__restrict ids, и это не меняет asm вообще, согласно панели diff в проводнике компилятора Godbolt: https://godbolt.org/z/4YjlDA. As мы ожидаем, потому что псевдонимы на основе типов позволяют компилятору предполагать, что хранилище double в bar[] не изменяет ни один из int * элементов int *ids[]. Как говорили люди в комментариях, нет алиасинг здесь, что компилятор уже не может разобраться. И на практике кажется, что сортирует его без каких-либо дополнительных перезагрузок указателей.

Он также не может иметь псевдоним *foo[k], потому что мы получили указатели с new внутри этой функции. Они не могут указывать внутрь bar[].

(Все основные компиляторы x86 C ++ (GCC, clang, ICC, MSVC) поддерживают __restrict в C ++ с тем же поведением, что и C99 restrict: обещание компилятору, которое хранится через этот указатель, не изменяет объекты на которые указывает другой указатель. Я бы порекомендовал __restrict вместо __restrict__, по крайней мере, если вы в основном хотите переносимости между компиляторами x86. Я не уверен насчет этого.)

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

---

Первая версия вопроса содержала ошибку во внутреннем цикле: в ней было K++ вместо k++, так что это было чисто неопределенное поведение, и компиляторы стали странными. Asm не имеет никакого смысла (например, нет инструкции умножения FP, даже когда foo[] была функцией arg). Вот почему для измерения массива рекомендуется использовать имя типа klen вместо K.

После исправления этого в связи с Godbolt все равно нет разницы в ассемблере с / без __restrict для всего, но это намного более вменяемое.

Кстати, сделав double *foo[] функцией arg, мы посмотрим на asm только для основного цикла. И вам на самом деле понадобится __restrict, потому что хранилище до bar[] может изменить элемент foo[][]. Это не происходит в вашей функции, потому что компилятор знает, что память new не указана никакими существующими указателями , но он не знал бы, что если бы foo была функцией arg.

---

Есть небольшая часть работы внутри цикла - расширение 32-битных int результатов с расширением знака до их использования в качестве индексов массива с 64-битными указателями. Это добавляет цикл задержки где-то там, но не цепочку мультипликативных зависимостей FP, переносимых циклом, так что это может не иметь значения. Вы можете избавиться от одной инструкции во внутреннем цикле на x86-64, используя size_t k=0; в качестве самого внутреннего счетчика цикла. L[] - это 32-битный массив, поэтому i*L[k] должен быть расширен до знака внутри цикла. Нулевое расширение от 32 до 64-разрядных происходит бесплатно на x86-64, поэтому i * (unsigned)L[k] сохраняет инструкцию movsx в цепочке депонирования указателей. Тогда внутренний цикл, который выполняет gcc8.2, - это вся необходимая работа, требуемая вашей мерзкой структурой данных / компоновкой. https://godbolt.org/z/bzVSZ7

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

Он также не может автоматически векторизовать, потому что данные не являются смежными. Вы не можете получить непрерывные исходные данные из цикла j или i. По крайней мере, i будет простым шагом без необходимости повторять ids[j][k].

Если вы генерируете foo[k][...] и bar[...] транспонированными, то есть индексируете с помощью foo[k][ i + L[k] * ids[j][k] ], тогда у вас будет непрерывная память в src и dst, чтобы вы (или компилятор) могли использовать умножения SIMD.