# What's the fastest implementation of elementwise vector multiplication in Fortran?

My fortran code contains lines like the following

integer, parameter :: dbp = kind(1.0d0)
integer, parameter :: n = 1 000 000

real(dbp) :: x(n), y(n), z(n)

y(:) = x(:) * z(:)


I would like to take advantage (if possible) of some optimised maths libraries to carry out this operation. I have found a lapack routine dgbmv which multiplies a matrix by a vector. This would suit my needs if I create a diagonal matrix such that

$$\left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) = \left( \begin{array}{ccc} x_1 & & & & \\ & x_2 & & & \\ & & \ddots & & \\ & & & x_n & \\ \end{array} \right)\left( \begin{array}{c} z_1 \\ z_2 \\ \vdots \\ z_n\end{array} \right)$$

But I don't know if this is the best way to go about calculating x(:)*z(:). Is there a more appropriate way?

• It looks like you're just doing an element-wise vector multiplication, which would be y = x*z' - no external library required. If you already have x as a vector instead of a matrix, just use it directly. Why do you need lapack for this? – cbcoutinho Dec 8 '16 at 16:54
• My motivation is optimisation. Currently I just use standard fortran, but I am wondering if using lapack (or some other library) would make this operation faster. – DJames Dec 8 '16 at 16:56
• No, it cannot be made faster by using a library. What you are doing is simply an element-wise vector multiplication, which is already a very fast instruction (and a very basic one). The only things that can make it faster is to use compiler options or parallelism. Compiling with -Ofast should help increase the speed of the multiplication. Otherwise, you could transform your element-wise vector multiplication into a loop and render that loop parallel with OpenMP. There is no case where creating a matrix to then apply the multiplication could be in any way faster. – BlaB Dec 8 '16 at 17:45
• Since you are doing element wise operation, you could check vectorization capabilities of your platform, like Intel AVX and use specialised instruction. – Moonwalker Dec 8 '16 at 20:45

The cost of the multiplication is almost insignificant compared to the cost of loading the data from memory (and writing it back). If you're worried about performance, you should be thinking about data locality. Perform more flops with each of your $x_i$ and $z_i$ values (if that's possible) when you load them from memory and before proceeding on.

Setting up a diagonal matrix will at best make no difference, but is more likely to disastrously degrade performance.

• I don't understand "Perform more flops with each of your 𝑥𝑖 and 𝑧𝑖 values (if that's possible) when you load them from memory and before proceeding on.". First can you explain how to perform more flops in a scalar product x_i * z_i if x_i and z_i` are fixed? And then, I don't see how increasing flops can increase speed. – nougako Oct 7 at 16:45
• @nougako It's possible that there are other calculations you need to perform (perhaps a bit later) with the same data. In that case, you can speed up the overall execution time by doing all the calculations on a given set of data when it is loaded into memory, instead of loading it twice. – David Ketcheson Oct 8 at 6:47
• I am currently suspecting if data (un)locality in my program is the reason it runs slow and I just posted a question on this stack. Would you please take a look there scicomp.stackexchange.com/questions/36065/…? – nougako Oct 8 at 19:46

You are doing an element wise operation between two vectors, so if is possible is better use functions designed for vectors.

The library Lapack is for linear algebra, but it is made for operations or methods with an high level respect this basic operation.

For an optimized form of your code you can try to use some compiler option as BlaisB suggested.

An other way can be to use more basic level library, similar Blas (note the Lapack is builded over Blas).

For example MKL has got v?mul that performs a vector-vector element wise multiplication. See also this question in MKL forum.

Rephrasing what other people have said, you can treat diagonal matrices as vectors, which will reduce the memory required to store them, and the number of computations unless you are using sparse routines. You can follow similar approaches with other structured matrices like tridiagonal matrices.