0

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!

TheRevanchist
  • 331
  • 1
  • 4
  • 12

2 Answers2

2

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.

Bonan
  • 96
  • 2
  • 5
0

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();
kangshiyin
  • 9,681
  • 1
  • 17
  • 29