2

Consider the following subroutine, which calculates a sparse matrix product in a Poisson equation solver.

SUBROUTINE mut_A_sparse(n, w, v)
    INTEGER, INTENT(IN) :: n
    REAL, INTENT(IN)    :: w(n*n)
    REAL, INTENT(OUT)   :: v(n*n)

    v = -EOSHIFT(w, -n) - EOSHIFT(w, -1) + 4*w - EOSHIFT(w, 1) - EOSHIFT(w, n)
    v((n + 1):n*n:n) = v((n + 1):n*n:n) + w(n:(n*n - 1):n)
    v(n:(n*n - 1):n) = v(n:(n*n - 1):n) + w((n + 1):n*n:n)
END SUBROUTINE mut_A_sparse

I'm not sure what this compiles to. Will the compiler store the intermediate result of -EOSHIFT(w, -n) to some temporary array, or calculate it on the fly for each index of v? The TA told me it is the latter case, but as far as I know, EOSHIFT is not ELEMENTAL, so how does the compiler know calculating it on the fly will not cause accidentally quadratic complexity?

Additionally, what about array slices? I wonder if the following creates an array temporary. I am asking because passing sliced arrays to functions/procedures will induce a copy during argument binding, but what about primitive operators like + - * /?

z(::2) = a(::2) + b(::2) + c(::2)

I guess I should read the assembly, but the assembly output is a complete mess beyond my comprehension.

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
nalzok
  • 14,965
  • 21
  • 72
  • 139
  • It is easier to look at optimized GIMPLE. – Vladimir F Героям слава Apr 29 '22 at 18:38
  • @VladimirFГероямслава Thanks for the tip, but I am using `mpiifort` which does not seem to support GIMPLE output? – nalzok Apr 29 '22 at 18:42
  • 3
    The standard won't say how this is done, only define the what the result of the calculation is - as long as the answer is correct the compiler can do what it wants. So without knowing exactly which implementation (compiler) you are using we can't answer this. Even then it may choose more than one method depending upon a number of circumstances. – Ian Bush Apr 29 '22 at 19:03
  • No, that is only for gfortran. The language-lawyer tag is probably not necessary as this will depend more on what indvidual compilers are able to do. – Vladimir F Героям слава Apr 29 '22 at 20:27

0 Answers0