Mixed-precision improvements

[nrseman]

This PR provides several improvements to the mixed-precision BiCGSTAB implementation:

1. Hand-optimized SPMV and ILU0/DILU implementations for 4x4 blocks
2. Hand-optimized SPMV and ILU0/DILU implementations for 2x2 blocks
3. Templatized fall-back routines for SPMV and ILU0/DILU for nxn blocks with n > 4
4. Mixed-precision implementation of GhostLastMatrixAdapter
5. Mixed-precision implementation of WellModelGhostLastMatrixAdapter
6. Parallel scalar product that takes advantage of the fact that ghost cells are stored last

The impact of all these changes are:

1. support for arbitrary block-sizes
2. much better parallel scaling
3. support for second-stage ILU0/DILU smoothers in the CPR+AMG preconditioner

For more information on how to run the opm simulator with mixed-precision, please see opm/simulators/linalg/mixed/READAME.md

UPDATE:
You can now select mixed-precision cpr from the command line using the --linear-solver=mixed-cprw option

[blattms]

jenkins build this please

[nrseman]

You can now select mixed-precision cpr from the command line using the --linear-solver=mixed-cprw option

[nrseman]

Just rebased on master branch

[alfbr]

Anything happening with this one?

[nrseman]

Sorry for the late reply, @alfbr. I just came back from vacation. @SoilRos, I know you have built and run the code behind this PR. Did you also start reviewing the code?

[SoilRos]

@nrseman I just built it to have a sense of the speed up, but I did not have a detailed look besides skimming its contents to switch on/off the vector dot product.

I got the impression this is not ready for review because it is still not rebased to master (parts of its contents are already merged) and there are some comments that seem to be for debugging or forgotten. Could you have a look at that?

I noticed that there is duplication of some algorithms where their only change are the vector length specific calls. Do you think is there a way to unify those? That would make these functions significantly more maintainable.

In my testing, the specialization of the scalar vector dot product seems to be fine in parallel but it is slower in sequential mode, at least for NORNE. I mentioned this to you before holidays, did you have a look at why is that?

[nrseman]

Welcome back @SoilRos!

Here are my responses to your comments:

I got the impression this is not ready for review because it is still not rebased to master (parts of its contents are already merged) and there are some comments that seem to be for debugging or forgotten. Could you have a look at that?

I have rebased the code on master twice since I created the PR. OPM is a fast-moving project and I do not rebase on an ongoing basis unless I am aware of breaking changes. If you run into a situation where rebasing is required to review the code, just drop a comment in the PR. To avoid any confusion in the future, I never submit a PR unless I consider the code ready for review. I have just rebased for a third time, so right now everything should be up to date. I'll try to find and delete the comments you are alluding to.

I noticed that there is duplication of some algorithms where their only change are the vector length specific calls. Do you think is there a way to unify those? That would make these functions significantly maintainable.

Specialization is the name of the game to squeeze performance out of the hardware. Can you please be more specific on which algorithms you are referring to? Only then can I assess whether or not they can be unified.

In my testing, the specialization of the scalar vector dot product seems to be fine in parallel but it is slower in sequential mode, at least for NORNE. I mentioned this to you before holidays, did you have a look at why is that?

Do you mean the opposite? Both the plot you sent me and my own testing shows that in sequential mode NORNE is faster when using my dot product instead of the default. In parallel there are only tiny performance differences. Of course, NORNE runs more or less entirely in-cache and detailed profiling shows significant MPI overhead which dwarfs any improvements to the dot product.

[SoilRos]

To avoid any confusion in the future, I never submit a PR unless I consider the code ready for review.

Alright, no problem, I was not aware this is how you operate.

Can you please be more specific on which algorithms you are referring to? Only then can I assess whether or not they can be unified.

I am talking about the mixed BiCGStab, ILU and DILU algorithms. Their (non-trivial) logic is repeated several times with virtually no changes between them, aside of a different function call names. To unify them, templates are your friend. There are many ways to do it, this is just an example:

template<int B>
static void prec_mapplyc(rec_t *P, double *x) {
  // switch specialized implementation
}

template<int B>
static void bsr_vmspmv(bsr_matrix *A, const double *x, double *y) {
  // switch specialized implementation
}

template<int B>
int bslv_pbicgstabm(bslv_memory *mem, bsr_matrix *A, const double *b, double *x)
{

    // ...
    mem->use_dilu ? prec_dilu_factorize<B>(P,A) : prec_ilu0_factorize<B>(P,A); // choose 
    // ...
    for(j=0;j<max_iter;j++)
    {
        prec_mapplyc<B>(P,q_j);                                          //q_j=P.q_j;
        bsr_vmspmv<B>(A,q_j,v_j);                                        //v_j= A.q_j

        // ...
    }
    prec_mapplyc<B>(P,x_j);                                       //x_j=P.x_j;
    // ...
}

Do you mean the opposite?

Yes, sorry. You are right, is the parallel version the slow one. Could you find a case where you parallel dot product is faster? Otherwise is hard to justify adding code that does the same thing.

Could you add some tests that touch this code during testing? For example, by running a system/regression test that already works, but with these solvers instead.

[SoilRos]

There seem to be a lot of comments from my side, but worry not, most of them are minor issues. All in all, the PR looks good and there is a good case to merge it. It is also nice to see more integration with the linear solver classes.

[SoilRos]

This looks like a bit set, right? If so, please use std::vector<bool> or std::vector<int>. Note that std::vector<bool> is a funny type and will allocate N bits instead of N sizeof(bool), thus saving storage, but it will give you proxy objects to modify them (i.e., you do not get a proper pointer for the actual data). In either case, that should be the right tool for bit fields in C++.

[nrseman]

Let's hold off on addressing this until we've had a discussion on index sets. Depending on the solution, this variable may disappear altogether.

[SoilRos]

Similar comments to this class, block_sizes may be private, and virtual keyword is superfluous.

[SoilRos]

@blattms may I ask you to review this class? I do not have enough experience with the parallel index sets to assess what is going on here with confidence.

[nrseman]

Fixed! Thanks!

[nrseman]

Regarding GhostLastScalarProduct:

Yes, sorry. You are right, is the parallel version the slow one. Could you find a case where you parallel dot product is faster? Otherwise is hard to justify adding code that does the same thing.

My testing shows that my implementation consistently beats the default implementation. Below are a few results for norne run with mixed-ilu0 and a tolerance of 1e-3. The differences are not big, but as we have already discussed, for norne the linear solver runs entirely in cache for ranks > 2 and there is a significant mpi overhead.

Results using GhostLastScalarProduct

ranks steps newton linear solve
-------------------------------
16    366   1488   36363  19.12
 8    363   1471   35870  27.11
 4    364   1500   35926  42.97
 2    363   1485   35531  93.29
-------------------------------

Results using default scalar product

ranks steps newton linear solve
-------------------------------
16    358   1468   36212  19.81
 8    362   1478   35865  28.36
 4    365   1489   35540  43.73
 2    361   1480   35382  96.25
-------------------------------

[nrseman]

Regarding unifying algorithms

I am talking about the mixed BiCGStab, ILU and DILU algorithms. Their (non-trivial) logic is repeated several times with virtually no changes between them, aside of a different function call names.

I'm not a big fan of using templates (or C++ for that matter) for performance critical code. The reason why is that it makes studying the generated assembly code a small nightmare. Instead, I would prefer to use other techniques to unify the algorithms:

• For BiCGStab, function pointers set at compile time should do the trick. Note that the current PR bypasses the C-implementation of BiCGStab with the ISTL version. The legacy C-implementation is there for reference, since it is still faster for sequential runs than ISTL.
• For ILU/DILU it is relatively straight forward to unify the ILU and DILU factorzations. Unifying implementations for different vector lengths while retaining performance and requires some more thought. Note that I used to have a generic implementation for this, I deemed it too slow.

I suggest we leave this for a future PR to avoid scope creep.

[nrseman]

Regarding adding tests

Could you add some tests that touch this code during testing? For example, by running a system/regression test that already works, but with these solvers instead.

Yes, I absolutely want to do this. In fact, it is something I have mentioned to Markus before. However, my preference is to do this through a dedicated PR.

[SoilRos]

Regarding tests, I think we can delay them given that your first implementation of mixed-precision already set a precedent. However, please do not make this a habit for future contributions. The intended workflow for collaboration is that PRs need to be accompanied with their respective tests right away, not later.

I strongly suggest that we address the unification of algorithms within the scope of this PR. Your suggestion to delay this is exactly how technical debt is accumulated and I am not comfortable with that here. Using templates was just a suggestion as we already use C++, but other options are also welcomed!

[SoilRos]

My testing shows that my implementation consistently beats the default implementation.

I got my results with a tolerance of 1e-4, but I do not think that explains our different numbers. In fact, we are getting contradictory results and, as you said, this mostly demonstrate why norne is not the best benchmark for this feature. What I am trying to say is: if the premise of this class is that it delivers faster results than the current approach, ideally, I would like to see a case where I can assert that very premise. I thought that maybe you already tried this in a bigger case, like SPE10, but if not, I will give it a try sometime this week.

[SoilRos]

I'd also like to check that the two new block sizes at least run on my side. What cases are you using to check those?

[nrseman]

I'd also like to check that the two new block sizes at least run on my side. What cases are you using to check those?

I use the following test cases from opm-tests which Markus brought to my attention on May 19:

• blocksize 5: SPE1CASE1_BRINE_THERMAL.DATA
• blocksize 4: spe9_solvent/SPE9_CP_SOLVENT_CO2.DATA
• blocksize 2: h2store/H2STORE.DATA

Matthew also provided scaled up version of the above cases for our joint ECMOR paper.

[nrseman]

Regarding tests

Regarding tests, I think we can delay them given that your first implementation of mixed-precision already set a precedent. However, please do not make this a habit for future contributions. The intended workflow for collaboration is that PRs need to be accompanied with their respective tests right away, not later.

The questions on tests came up during the previous mixed-precision PR. I outlined the following plan in an email to Markus and yourself on July 8:

Hi Markus,

I have not added any tests yet. The current PR is just a partial progress towards the current goal of absorbing mixed-precision into the ISTL-framework. There are at least two more PRs pending after we are done with this one. I will add a test in the last PR, but I still need to learn how to do it.

SINTEF expressed some concern last time that my PR was too big and Equinor recommended that I submit more frequent PRs to speed-up the process. My current approach is an attempt to follow their advice.

Thanks,
KJ

The current PR is the first of the two additional PRs referred to in the email. Given that there were no objections at the time, I am simply following the plan. I appreciate that you are not flipping the tables on me now.

[nrseman]

Regarding GhostLastScalarProduct:

I got my results with a tolerance of 1e-4, but I do not think that explains our different numbers. In fact, we are getting contradictory results and, as you said, this mostly demonstrate why norne is not the best benchmark for this feature. What I am trying to say is: if the premise of this class is that it delivers faster results than the current approach, ideally, I would like to see a case where I can assert that very premise. I thought that maybe you already tried this in a bigger case, like SPE10, but if not, I will give it a try sometime this week.

I did test with spe10, but I did not save the results. I reran the 8-way comparison last night for mixed-ilu0 at a tolerance of 1e-4 and with maximium number of linear iterations set to 1024. The results below show an 8.3% speed-up, but that also accounts for the slight reduction in linear iterations. The speed-up per linear iteration is 4.4%. There is of course no guarantee that the perturbations caused by altering the order in which the summation is performed will lead to a reduction in linear iterations. Thus, I am always looking at the per-linear-iteration speed-up.

dot        ranks  steps  newton  linear  solve
------------------------------------------------
ghostLast  8      75     509     166738  5491.22
default    8      77     535     172968  5948.88
------------------------------------------------

I also ran norne with mixed-ilu0 at a tolerance of 1e-4 to see if I can reproduce your findings. As you can see from the results below, all simulations are faster with the new scalar product.

Results using GhostLastScalarProduct

ranks steps newton linear solve
--------------------------------
16    347   1204   45494   21.82
 8    345   1211   45447   32.24
 4    349   1214   44573   49.50
 2    347   1217   44939  111.86
--------------------------------

Results using default scalar product

ranks steps newton linear solve
--------------------------------
16    348   1209   45426   23.08
 8    346   1209   45360   33.78
 4    347   1217   44610   51.95
 2    348   1200   44385  113.68
--------------------------------

[nrseman]

Regarding unifying algorithms:

I strongly suggest that we address the unification of algorithms within the scope of this PR. Your suggestion to delay this is exactly how technical debt is accumulated and I am not comfortable with that here. Using templates was just a suggestion as we already use C++, but other options are also welcomed!

I just pushed a commit that unifies the ILU0/DILU factorization algorithms.

Technical debt is to a large extent subjective and I have yet to work on a project that does not allow some degree of code duplication if it enhances performance or simplifies the mental map required to work in the code. OPM has relatively rich set of preconditioners and linear operators, many of which could have been unified, but are not, so I think there is precedence for my take on this. As the most likely maintainer of this code, I assume my perception of technical debt counts too?

PS! Note that the c-implementations of bicgstab are intended to be completely replaced by ISTL once the performance of its sequential implementation is competitive. As such I do not see a need to recreate what ISTL provides in c.

[nrseman]

Regarding unifying algorithms:

I just pushed a commit that uses function pointers to provide a generic implementation of the legacy bicgstab

[SoilRos]

Thanks for all the nice updates and clean ups, that's very much appreciated!

I got the different blocks working for me locally and while compiling I found a warning that should be easy to fix (see comment below: #7224 (comment)). So that looks good. I am currently trying to set up the SPE10 but I am having troubles with a keyword that is unrelated to this PR. Your results already look good though. You are running the cartesian model 2 from opm-test, right?

[nrseman]

I am currently trying to set up the SPE10 but I am having troubles with a keyword that is unrelated to this PR. Your results already look good though. You are running the cartesian model 2 from opm-test, right?

I am primarily running model2 from opm-data, but I have run the same case from opm-tests as well. The results I shared are from the former.

You may run into the same negative endpoint issue I did. Seems like a bug. My ugly hack is captured by the diff below. That got the case running.

diff --git a/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp b/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
index df45ad042..c07db44e7 100644
--- a/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
+++ b/opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
@@ -31,7 +31,7 @@ template <typename Scalar>
 void Opm::Satfunc::PhaseChecks::Oil::SOcr_GO<Scalar>::
 testImpl(const EclEpsScalingPointsInfo<Scalar>& endPoints)
 {
-    this->sogcr_ = endPoints.Sogcr;
+    this->sogcr_ = 0;//endPoints.Sogcr;
 
     if (! std::isfinite(this->sogcr_)) {
         this->setViolated();

[nrseman]

@SoilRos, I believe all issues have been addressed. Do you mind kicking off the CI? @blattms, would you be available to review the use of index sets in the MixedGhostLast linear operators and the GhostLastScalarProduct. I think that is the last item missing...

[SoilRos]

jenkins build this please

[SoilRos]

@nrseman The CI seems to be passing with no problem.

Just for reference. @nrseman and I met and walked through the differences in how we made our measurements. It turns out that our BIOS settings were different, which really affected the --mab-by numa command option for the MPI call. In my settings, I only have one numa node for the entire socket whereas @nrseman enables the AMD subdivision of CCDs to be their own NUMA domain while also disabling the hyper-threads. Then, mapping by numa domains in @nrseman case exhausts memory bandwidth and last level caches of each CCD as the program scales, whereas in my case processes are arbitrarily placed, most likely oversubscribing the available CCD resources. In the end, we achieved essentially the same results as @nrseman by manually pinning processes to each CCD with the --map-by l3cache:hwtcpus --bind-to hwthread option.