looping on flattened fortran matrices

Question

I have a bit of code that looks like this:

DO I=0,500
    arg1((I*54+1):(I*54+54)) = premultz*sinphi(I+1)
ENDDO

In short, I have an array premultz of dimension 54. I have an array sinphi of dimension 501. I want to take the first value of of sinphi times all the entries of premultz and store it in the first 54 entries of arg1, then the second value of of sinphi times all the entries of premultz and store it in the second54 entries of arg1, and so on.

These are flattened matrices. I have flattened them in the interest of speed, as one of the primary goals of this project is very fast code.

My question is this: is there a more efficient way of coding this sort of calculation in Fortran90? I know that Fortran has a lot of nifty array operations that can be done that I'm not fully aware of.

Thanks in advance.

Answer 1

This expression, if I've got things right, ought to create arg1 in one statement

arg1 = reshape(spread(premultz,dim=2,ncopies=501)*&
              &spread(sinphi,dim=1,ncopies=54),[1,54*501])

I've hardwired the dimensions here, that may or may not suit your purposes. The inner expression generates the outer product of premultz and sinphi, which is then reshaped into a vector. You may find you need to reshape the transpose of the outer product, I haven't checked things very carefully.

However, based on my experience with this sort of clever use of Fortran's array intrinsics I doubt that this, or most other clever uses of Fortran's array intrinsics, will outperform the straightforward loop implementation you already have. For many of these operations the compiler is going to generate copies of arrays, and copying data is relatively expensive. Of course, this is an assertion you may want to test.

I'll leave it to you to decide if the one-liner is more comprehensible than the loops. Sometimes the expressivity of the array syntax comes at an acceptable cost in performance, sometimes it doesn't.