Multi-D pointer with constant stride and Rank R-1 into array of Rank R

1

Hi there,
I have an equation of state table that is a 4D table of the kind stateArray(x,rho,T,Ye),
where x is the state variable, e.g., entropy, potential,…, while (rho,T,Ye) are the grid variables on which the values are tabulated.
In order to use SIMD to interpolate all state variables at once given a particular (rho,T,Ye) value I have them as the fastest running index.

For convenience however I also want to have direct access to e.g. all entropy values as a function of (rho,T,Ye) so I also have a 3D-Array of Entropy(rho,T,Ye), etc…

Now in order to save space I use 3D pointers of constant stride into the 4D array, i.e.,

Entropy(1:nro,1:ntt,1:nye) => stateArray(1,1:nro,1:ntt,1:nye)
Pressure(1:nro,1:ntt,1:nye) => stateArray(2,1:nro,1:ntt,1:nye)

All state values should have a constant stride in the stateArray because each state value is just nStateVariable away from the previous one.

Intel Fortran is able to handle it no problem, however GFortran complains that:
“Error: Rank remapping target must be rank 1 or simply contiguous at (1)”

How can I rewrite this pointer so that it works in both compilers, or is this just a case of Intel Fortran being less strict about standards?

2 
implicit none
real, target :: a(3,5,5,5)
real, pointer :: b(:,:,:) => null()
b => a(1,1:5,1:5,1:5)
end

This seems to work: https://godbolt.org/z/4dWWscssv

1 Like
3

Thanks Ivan,
I might have wrongly remembered that syntax not working before.

4

Welcome to the Discourse @RobertB !

@ivanpribec 's answer should solve the issue. However, it appears that Gfortran, contrarily to ifort/ifx, performs a copy before the pointer association. So, be careful that even though you would get to a compiling version, you still have to check it is translated in what you seek.

Why don’t you split the stateArray into each one of the state variables as you also show, and then you construct when needed the array of state variables given a particular (rho,T,Ye).

Do you do it for all possible values of (rho,T,Ye)?

5

Hi @mEm,
I typically trilinearly interpolate all the state variables on an array of nr*(rho(i),T(i),Ye(i)) state coordinates.
There are about 32 state variables that I can then read as a contiguous dataset of 32 doubles at each 8 of the rho(i),T(i),Ye(i) locations. Where i is the radial index in an r,theta,phi coordinate system that is parallelized over theta and phi.
The point is to do only 8*nr random lookups, compared to 32*8*nr random lookups, with enough contiguous data to do SIMD operations on.

Hmm, what would it mean for it do a copy before associating the pointer. Wouldn’t the temporary copy be lost once the routine is exited and the stack is cleared? Or does it secretly allocate contiguous memory for the pointer and copy into that, which is of course not what i want.
Do you see the copy in the compiler assembly output?
I only see a bunch of mov commands that seem to move data into the pointer array descriptor.

6

It’s surprising, if it happens it’s a bug… Do you have some evidence of that (some sample code on godbolt for instance)?

Just pointing to the higher dimensions won’t magically enable SIMD… What is important is that the vector on which you would like SIMD to operate, have a stride equal to 1 (i.e. the elements are contiguous in memory). Pointing to non-contiguous memory doesn’t make it contiguous, the internal descriptor of the pointer will have a stride >1 for the first dimension…

I’ve read that some SIMD instructions are able to work on strided data, but I’d rather not cound on it in the general case. Note that if even for such instructions, the pointer trick would not be necessary…

If possible you should build your array with x in the last dimension stateArray(rho,T,Ye,x)

1 Like
7

Just for clarification, I mostly use the 4D stateArray for interpolation of the form (:,rho,T,Ye) which should have contiguous data in the first coordinate at all 8 interpolation grid points. The lookup for the particular combination of rho,T,Ye would be a mostly random access in memory space.
The 3D arrays are only used sporadically when I only want to interpolate a single quantity or for backward compatibility to the older implementation.

8

Upon checking, I do admit my previous statement is wrong ! So please ignore it.
Since I don’t want to hide the mistake, but at the same time I don’t want next readers to get misled by it, I recall somewhere (here I guess) having seen crossed-out sentences, to remark they are faulty.
If anyone here could suggest how to do it, I will happily implement it.

2 Likes