Optimizations method weighted

[CarlosPenaDePedro]

This pull request highlights changes that can improve the performance of the code. The first change involves iterating over sparse matrices using indexes instead of iterators, which allows for better compiler optimizations and is generally more suited for structures with a known size, such as vectors.

The second change introduces OpenMP parallelization for certain independent sections of the code, such as the treatment of rows, to further improve performance.

These changes are applied only to specific parts of the code relevant to our current interests. However, we understand that if these changes are accepted, a major refactor would be required to standardize their usage across all parts of the code that rely on this iterator structure. Such a refactor would represent a long-term development effort.

We also acknowledge that the iterator structure may have a specific purpose and might not be replaced. Our goal is to showcase the observed improvements on MN5 and explain how we achieved them.

We applied these changes to the ForceMissing calculation (Section 3) and the MissingValueTreatment, specifically the MissingIfHeaviestMissing case (Section 9).

Other relevant sections include:

1. Section 8, which corresponds to the matrix copy and represents the major bottleneck, as we could not modify it within MIR. Perhaps the eckit copy constructor could be revised to optimize this process.
2. Section 10, which covers the interpolation calculation itself. This section is already parallelized using OpenMP due to ECKIT_OMP and MIR_OMP being enabled.

This table presents the average burst time for the MIR Method Weighted sections across different servers during an output step of Tco1279-eORCA12 simulations using the DE340 output plans. The data focuses exclusively on the ocean servers, showcasing only eORCA12 data with missing values.

This is the base case.

Using indexes instead of iterators

• Section 3: 5.4x Average burst time speedup respect base
• Section 9: 2.85x Average burst time speedup respect base

Using indexes and 16 OMP threads

• Section 3: 29.42x Average burst time speedup respect base
• Section 9: 18.97x Average burst time speedup respect base

[FussyDuck]

All committers have signed the CLA.

[CarlosPenaDePedro]

This part corresponds to Section 3. If I am not mistaken, this calculation could be performed during matrix creation, as the original empty rows in the weight matrix are part of the geometry. A vector containing the indices of these empty rows could be stored as a private member and cached. This approach could save significant time, especially with larger matrices, where traversing all rows is costly. For example, in the base case profiling with eORCA12, each MIR API call takes approximately 40 ms.

[pmaciel]

The OMP header should look into the conditional pre-processor define (defined in #include "mir/api/mir_config.h")
#if mir_HAVE_OMP

And in case there's no OpenMP found the pragma, or defined symbol, should be doing nothing

[wdeconinck]

In atlas we have something just for that :

#if ATLAS_HAVE_OMP
#define ATLAS_OMP_STR(x) #x
#define ATLAS_OMP_STRINGIFY(x) ATLAS_OMP_STR(x)
#define atlas_omp_pragma(x) _Pragma(ATLAS_OMP_STRINGIFY(x))
#else
#define atlas_omp_pragma(x)
#endif

With something similar you can then use:

mir_omp_pragma( omp parallel for reduction(vec_merge_sorted:forceMissing) )

Or define the reduction:

#define mir_omp_parallel_for_reduction(x) mir_omp_pragma( omp parallel for reduction(x) )

and then use

mir_omp_parallel_for_reduction(vec_merge_sorted:forceMissing)
for( ... ) {...}