CODEWITHSUNDEEP
cvectorgnublaselementwise-operations
TheRevanchist
TheRevanchist

Reputation: 331

Elementwise product between a vector and a matrix using GNU Blas subroutines

I am working on C, using GNU library for scientific computing. Essentially, I need to do the equivalent of the following MATLAB code:

x=x.*(A*x);

where x is a gsl_vector, and A is a gsl_matrix.

I managed to do (A*x) with the following command:

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

where res is an another gsl_vector, which stores the result. If the matrix A has size m * m, and vector x has size m * 1, then vector res will have size m * 1.

Now, what remains to be done is the elementwise product of vectors x and res (the results should be a vector). Unfortunately, I am stuck on this and cannot find the function which does that.

If anyone can help me on that, I would be very grateful. In addition, does anyone know if there is some better documentation of GNU rather than https://www.gnu.org/software/gsl/manual/html_node/GSL-BLAS-Interface.html#GSL-BLAS-Interface which so far is confusing me.

Finally, would I lose in time performance if I do this step by simply using a for loop (the size of the vector is around 11000 and this step will be repeated 500-5000 times)? 

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

Thanks!

Upvotes: 0

Views: 343

Answers (2)

Bonan
Bonan

Reputation: 96

The function you want is:

gsl_vector_mul(res, x)

I have used Intel's MKL, and I like the documentation on their website for these BLAS routines.

Upvotes: 2

kangshiyin
kangshiyin

Reputation: 9781

The for-loop is ok if GSL is well designed. For example gsl_vector_set() and gsl_vector_get() can be inlined. You could compare the running time with gsl_blas_daxpy. The for-loop is well optimized if the timing result is similar.

On the other hand, you may want to try a much better matrix library Eigen, with which you can implement your operation with the code similar to this

x = x.array() * (A * x).array();

Upvotes: 0

Related Questions