Поэлементное произведение между вектором и матрицей с использованием подпрограмм GNU Blas

Я работаю над C, используя библиотеку GNU для научных вычислений. По сути, мне нужно сделать эквивалент следующего кода MATLAB:

x=x.*(A*x);

где x — gsl_vector, а A — gsl_matrix.

Мне удалось сделать (A*x) с помощью следующей команды:

gsl_blas_dgemv(CblasNoTrans, 1.0, A, x, 1.0, res);

где res — это другой gsl_vector, в котором хранится результат. Если матрица A имеет размер m * m, а вектор x имеет размер m * 1, то вектор res будет иметь размер m * 1.

Теперь осталось произвести поэлементное произведение векторов x и res (в результате должен получиться вектор). К сожалению, я застрял на этом и не могу найти функцию, которая это делает.

Если кто-нибудь может помочь мне в этом, я был бы очень благодарен. Кроме того, кто-нибудь знает, есть ли документация по GNU лучше, чем https://www.gnu.org/software/gsl/manual/html_node/GSL-BLAS-Interface.html#GSL-BLAS?-Interface, который до сих пор меня смущает .

Наконец, потеряю ли я время выполнения, если выполню этот шаг, просто используя цикл for (размер вектора составляет около 11000, и этот шаг будет повторяться 500-5000 раз)?

for (i = 0; i < m; i++)
    gsl_vector_set(res, i, gsl_vector_get(x, i) * gsl_vector_get(res, i));

Спасибо!


c vector gnu elementwise-operations blas
person TheRevanchist    schedule 14.07.2016    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Функция, которую вы хотите:

gsl_vector_mul(res, x)

Я использовал Intel MKL, и мне нравится документация на их веб-сайте для этих процедур BLAS.

person Bonan    schedule 14.07.2016
comment
Действительно, это, кажется, работает именно так, как я хотел. - person TheRevanchist; 16.07.2016

arrow_upward
0
arrow_downward

Цикл for подходит, если GSL хорошо спроектирован. Например, gsl_vector_set() и gsl_vector_get() могут быть встроенными. Вы можете сравнить время работы с gsl_blas_daxpy. Цикл for хорошо оптимизирован, если результат синхронизации аналогичен.

С другой стороны, вы можете попробовать гораздо лучшую библиотеку матриц Eigen, с помощью которого вы можете реализовать свою операцию с кодом, подобным этому

x = x.array() * (A * x).array();
person kangshiyin    schedule 14.07.2016