Эффективное вычисление продуктов кронекера в C

Я довольно новичок в C, но не нуждаюсь в чем-то быстрее, чем python для большинства моих исследований. Однако выясняется, что недавняя работа, которую я выполнял, требовала вычисления довольно больших векторов / матриц, и поэтому решение C + MPI могло бы быть в порядке.

Математически говоря, задача очень проста. У меня много векторов размерности ~ 40k и вы хотите вычислить произведение Кронекера выбранных пар этих векторов, а затем суммировать эти продукты кронекера.

Вопрос в том, как это сделать эффективно? Есть ли что-то неправильное в следующей структуре кода, используя для циклов или получить эффект?

Описанная ниже функция kron передает векторы A и B длин vector_size и вычисляет их произведение kronecker, которое он хранит в C , vector_size*vector_size .

 void kron(int *A, int *B, int *C, int vector_size) { int i,j; for(i = 0; i < vector_size; i++) { for (j = 0; j < vector_size; j++) { C[i*vector_size+j] = A[i] * B[j]; } } return; } 

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

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

    Для векторов с двойной точностью (одноточечная и сложная аналогичны), вы можете использовать BLAS-процедуру DGER (обновление первого ранга) или аналогично делать продукты по одному, так как они все на векторах. Сколько векторов вы умножаете? Помните, что добавление пучка векторных внешних продуктов (которые вы можете рассматривать как продукты Kronecker as) заканчивается как умножение матрицы-матрицы, которое DGEMM BLAS может эффективно обрабатывать. Возможно, вам придется писать свои собственные подпрограммы, если вам действительно нужны целые операции.

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

     #pragma omp parallel for for(int i = 0; i < vector_size; i++) { for (int j = 0; j < vector_size; j++) { C[i][j] = A[i] * B[j]; } } 

    В настоящее время это поддерживается многими компиляторами.

    Вы также можете попытаться вытащить некоторые распространенные выражения из внутреннего цикла, но достойные компиляторы, например, gcc, icc или clang, должны сделать это довольно хорошо самостоятельно:

     #pragma omp parallel for for(int i = 0; i < vector_size; ++i) { int const x = A[i]; int * vec = &C[i][0]; for (int j = 0; j < vector_size; ++j) { vec[j] = x * B[j]; } } 

    BTW, индексирование с помощью int обычно не подходит. size_t - это правильный typedef для всего, что связано с индексацией и размерами объектов.

    Если ваш компилятор поддерживает C99 (и вы никогда не пропускаете тот же вектор, что и A и B ), подумайте о компиляции в поддерживающем C99 режиме и изменении вашей сигнатуры функции:

     void kron(int * restrict A, int * restrict B, int * restrict C, int vector_size); 

    Ключевое слово restrict обещает компилятору, что массивы, на которые указывают A , B и C , не являются псевдонимами (перекрываются). При написании вашего кода, компилятор должен перезагрузить A[i] при каждом выполнении внутреннего цикла, потому что он должен быть консервативным и предположить, что ваши магазины в C[] могут изменять значения в A[] . Под restrict , компилятор может предположить, что этого не произойдет.

    Решение найдено (спасибо @Jeremiah Willcock): привязки BLAS от GSL, похоже, прекрасно делают трюк. Если мы постепенно выбираем пары векторов A и B и добавляем их к некоторому «текущему полному» вектору / матрице C , следующая модифицированная версия указанной функции kron

     void kronadd(int *A, int *B, int *C, int vector_size, int alpha) { int i,j; for(i = 0; i < vector_size; i++) { for (j = 0; j < vector_size; j++) { C[i*vector_size+j] = alpha * A[i] * B[j]; } } return; } 

    точно соответствует функции BLAS DGER (доступной как gsl_blas_dger ), функционально говоря. Начальная функция kron - это DGER с alpha = 0 а C - неинициализированная (обнуленная) matrix / вектор правильной размерности.

    Оказывается, в конечном итоге было бы проще просто использовать привязки python для этих библиотек. Однако, я думаю, я многому научился, пытаясь понять это. В других ответах есть несколько полезных советов, проверьте их, есть ли у вас такая же проблема. Спасибо всем!

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

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

    Все вышеперечисленное (вероятно) будет быстрее, чем код, написанный строго на python. Если вам нужна более высокая скорость, я бы предложил пару вещей:

    1. Посмотрите на использование Fortran вместо C. Компиляторы Fortran, как правило, лучше оптимизируют числовые вычисления (одним исключением было бы, если вы используете gcc, поскольку оба его компилятора C и Fortran используют один и тот же бэкэнд).
    2. Рассмотрим распараллеливание вашего алгоритма. Есть варианты Fortran, которые я знаю, которые имеют инструкции параллельного цикла. Я думаю, что вокруг есть некоторые C-дополнения, которые делают то же самое. Если вы используете ПК (и одноточность), вы также можете использовать графический процессор вашей видеокарты, который по сути является действительно дешевым процессором массива.

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

    Эту страtagsю можно обобщить, используя оператор switch во внешнем цикле, при этом размеры массивов делятся на два, три, четыре и пять, или что-то другое. Это может дать большой выигрыш в производительности и совместим с предложениями 1 и 3 для дальнейшей оптимизации / параллелизации. Хороший компилятор может даже сделать что-то подобное для вас (он же разворачивается).

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

     int i, j; for(i = 0; i < vector_size; i++) { int d = *A++; int *e = B; for (j = 0; j < vector_size; j++) { *C++ = *e++ * d; } } 

    Это также позволяет избежать доступа к значению A [i] несколько раз, кэшируя его в локальной переменной, что может дать вам небольшое ускорение скорости. (Обратите внимание, что эта версия не является параллельной, поскольку она изменяет значение указателей, но все равно будет работать с разворачиванием цикла.)

    Чтобы решить вашу проблему, я думаю, вы должны попробовать использовать Eigen 3, это библиотека C ++, которая использует все функции матрицы!

    Если у вас есть время, перейдите к документации! знак равно

    Удачи !