2
$\begingroup$

I am solving the linear elasticity equation within a FEM library with a complex 3D geometry. The resulting linear system is solved with CG, preconditioned by AMG (Algebraic Multigrid). The computed numerical solution is correct. However, when running the program with different numbers of processors, say going from 5 to 10 MPI ranks, the of iterations count slightly changes from 300 to 303.

The solution looks identical of course, but I am quite puzzled by this fact. Is this something that can happen with AMG when running in parallel? I am using TrilinosML, so the AMG implementation is kind of a black box to me right now.

Improve this question
$\endgroup$
2
  • $\begingroup$ Yes that's okay $\endgroup$
    – CuteCompute
    Sep 28 at 6:33
  • $\begingroup$ Do you get exactly the same results every time you run the code with 10 MPI ranks? $\endgroup$
    – Brian Borchers
    Sep 28 at 22:55

2 Answers 2

Reset to default
12
$\begingroup$

This is something that can happen in almost any numerical algorithm running in parallel.

It's important to know that floating-point addition is not associative due to round-off errors. Thus you can't count on

$(a+b)+c=a+(b+c)$

Different orders of operations can then result in answers that are different and you can then see different numbers of iterations used in iterative methods.

Parallel algorithms that distribute load across processors typically do not ensure that all floating point arithmetic operations are performed in exactly the same order regardless of the number of processors used. In distributed-memory message-passing systems, there is often non-determinism in the time it takes for messages to pass from processor to processor, so you might not get exactly the same answer in repeated executions of the code on the same number of processors.

Improve this answer
$\endgroup$
3
  • $\begingroup$ I get the part about non-associativity of floating point; but how does non-determinism in the time it takes for messages to pass from processor to processor play a part? $\endgroup$
    – NNN
    Sep 28 at 7:38
  • 4
    $\begingroup$ @NNN if you read the MPI standard for (say) MPI_REDUCE it says 'The “canonical” evaluation order of a reduction is determined by the ranks of the processes in the group. However, the implementation can take advantage of associativity, or associativity and commutativity in order to change the order of evaluation. This may change the result of the reduction for operations that are not strictly associative and commutative, such as floating point addition." I.e. your implementation can run a summation at a single memory location based on order of arrival. $\endgroup$
    – origimbo
    Sep 28 at 15:09
  • 1
    $\begingroup$ Even on a shared memory multiprocessor system, it can be extremely hard to ensure that operations are always performed in precisely the same order in a repeated run of the same code. See for example the "Conditional Numerical Reproducibility" feature of the Intel MKL library. $\endgroup$
    – Brian Borchers
    Sep 28 at 22:54
4
$\begingroup$

@brianborchers has already given one answer. Here is another: When you throw a matrix into an AMG algorithm, it needs to figure out a coarsening strategy whereby it determines a smaller linear system that approximates the big linear system you are trying to solve. It will preferentially do that by merging degrees of freedom located on the same process, until the local problems become so small that merging needs to happen across processes.

It should be clear that this strategy introduces a dependency on the number of processes you have overall, and you see that in the number of iterations required to solve the problem.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.