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

Question

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

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

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

Функция kron, описанная ниже, передает векторы A и B длины vector_size и вычисляет их произведение Кронекера, которое хранится в 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, но поиск в Google принес мне немного радости для этого запроса.

Answer 1

Поскольку ваши тела цикла полностью независимы, есть способ ускорить это. Проще всего уже воспользоваться преимуществами нескольких ядер, прежде чем думать о 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];
    }
}

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

Answer 2

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

Answer 3

Решение найдено (благодаря @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 являются неинициализированной (обнуленной) матрицей / вектором правильной размерности.

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

Answer 4

Если ваш компилятор поддерживает 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 компилятор может предположить, что этого не произойдет.

Answer 5

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

Эта стратегия может быть обобщена путем использования оператора 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], кэшируя его в локальной переменной, что может дать незначительный прирост скорости. (Обратите внимание, что эта версия не распараллеливается, поскольку она изменяет значение указателей, но все равно будет работать с развертыванием цикла.)

Answer 6

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

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

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

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

Answer 7

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

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

Удачи!