2
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • 3
    $\begingroup$ 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? $\endgroup$
    – cbcoutinho
    Dec 8, 2016 at 16:54
  • $\begingroup$ 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. $\endgroup$
    – DJames
    Dec 8, 2016 at 16:56
  • 7
    $\begingroup$ 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. $\endgroup$
    – BlaB
    Dec 8, 2016 at 17:45
  • 3
    $\begingroup$ Since you are doing element wise operation, you could check vectorization capabilities of your platform, like Intel AVX and use specialised instruction. $\endgroup$
    – Moonwalker
    Dec 8, 2016 at 20:45

3 Answers 3

Reset to default
7
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – nougako
    Oct 7, 2020 at 16:45
  • 1
    $\begingroup$ @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. $\endgroup$
    – David Ketcheson
    Oct 8, 2020 at 6:47
  • $\begingroup$ 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/…? $\endgroup$
    – nougako
    Oct 8, 2020 at 19:46
2
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.