How to achieve CPU performance

[Try2Code]

hi!

Background: I am working with a large Fortran code base and the team I am working with will replace Fortran-OpenACC with cpp/kokkos. so we try to solve a similar technical problem like EAMxx did. Couple of years ago @PeterCaldwell gave a talk on what was called scream at that time during a WarmWorld kickoff meeting in Bonn, Germany. Besides the many social and strategic aspects, he mentioned special calls where needed to achieve decent performance with Kokkos on CPUs using OpenMP.

Since we now started to replace fortran with cpp we face very similar problems I think. My current questions are:

• Is CPU performance still relevant?
• If yes: What is currently implemented to have fortran-like CPU performance?
• What was the strategy in EAMxx in the code transistion regarding the memory layout going from fortran to cpp/kokkos? classic cpp loops would need a reverse loop order compared to fortran, but kokkos offeres something like Kokkos::LayoutLeft/Kokkos::LayoutRight

I would be very happy to get some feedback. If needed I can set up video calls for a more interactive form of discussion.

my target model is the ICON model (www.icon-model.org)

[bartgol]

A few thoughts:

• We were able to keep the same layout (namely, LayoutRight, or "C-style") for all archs, since in most of EAMxx we use hierarchical parallelism over a few dimensions. As such, when it's time to dispatch the inner for over a warp, we are already dealing with 1d views. In our case, views are 2d (ncol,nlev), and in the outer layer of the nested for we subview them at a particular col index, leaving just (nlev) views, and for 1d views, layout left or right is the same. This is different from a flat parallel for (or, in kokkos lingo, a parallel for with RangePolicy exec policy). In the "range" for case, a 2d array should have LayoutLeft on GPU in order to get coalesced access, and LayoutRight on CPU for cached access. In eamxx, we do have a few of these kernels, but the vast majority of the code uses hierarchical parallelism, which means layout right everywhere is still best.
• CPU performance IS important. We don't know what tomorrow will bring, but even today, one could run low-resolution ensembles (e.g., for ML purposes) on CPU. So running efficiently on CPU can still be important. Also, our application is mostly memory bound, so a CPU machine with fast memory (HBM) can still be a reasonable option.
• On CPU, the best way to get FLOPS/W (flops per watt) is to use vector instructions. Fortran is great at vectorization, as it's built for array operations (e.g., stuff like a(:,:)=b(:,:,k) works well), but C/C++ have a harder time, especially in the case of multiple layers of nested inline functions, where the compiler may give up and not vectorize. Several apps are coming up with their version of a "simd-like" data structure. You can prob recycle what kokkos has, or write your own (it's not a ton of code). In eamxx, we wrote our own, since we also wanted masked vector operation (where the op is only applied for entries where the corresponding entry of a "mask" vector is 1), which kokkos' simd struct did not support at the time. We did some study and noticed that a pack size of N or 2*N (where N is the vector length supported by the cpu) offers the best speed up. E.g., on the old KNL we had a speedup of about 3.5 when using N=8, IIRC. All that been said, there are algorithms that are not very well proned to vectorization, so you may a) do away with trying to use simd types or b) try to reformulate the algorithm (or its impl) so that you can.

[Try2Code]

thx @bartgol for considering such an extensive answer!

ICON uses an nproma-based approach, where 1d-horizontal fields are split up into blocks of nproma size chunks. By setting nproma at runtime, the binary can be tuned for different cache sizes, OpenMP-threads and for different grids. For GPUs we use large values so that at best only a single block exists. But the 3d-fields are declared as fields(1:nproma,1:nlev,1:nblock). Hence the loops are written according to this specific memory layout. this can cause issues with Fortran + OpenACC, which is our current production solution for GPUs. With some tweaks the code works fine on vector machines, though.

We face the same problem when it comes to hardware choices: CPU performance is still important but there is no clear winner for the future. That's why Fortran/OpenACC needs to be replaced by a more portable solution. we cannot change the memory layout all at once because main parts of the model will stay in Fortran for a while. Can you maybe point me to the code, that shows how you the interaction with cpp and fortran code is done? Do you use a dedicated memory manager for both Fortran and cpp?

thx in advance

[rljacob]

Everything was converted to C++. Physics, dynamics, radiation, everything. Then the flexibility to go between CPU and GPU friendly memory layouts comes from the C++ framework.

[bartgol]

The dycore and the surface coupling are the 2 places where we still use some F90 arrays: in the dycore, the init phase (which involves A LOT of code) is still done in F90, so we share arrays, and in the SC we must use f90 cause that's the lang of the coupler. For the cpl, you may look at this file, which imports data from the cpl f90 arrays into the component cxx ones. It's not an easy read though...

But the bottom line is that we do keep everything on device as much as possible.

[bartgol]

Ah, there are other places where we interact with f90 or with code with different layouts:

• in the cosp code we grab "host mirrors" of our fields, and pass them to f90. But the layout is the same
• in the radiation code we do transpose some arrays, so that the rad package can work (thought the pkg itself is still in C++)