  The Preconditioned Conjugate Gradient method is often employed for the
solution of linear systems of equations arising in numerical simulations of
physical phenomena. While being widely used, the solver is also known for its
lack of accuracy while computing the residual. In this article, we propose two
algorithmic solutions that originate from the ExBLAS project to enhance the
accuracy of the solver as well as to ensure its reproducibility in a hybrid MPI
+ OpenMP tasks programming environment. One is based on ExBLAS and preserves
every bit of information until the final rounding, while the other relies upon
floating-point expansions and, hence, expands the intermediate precision.
Instead of converting the entire solver into its ExBLAS-related implementation,
we identify those parts that violate reproducibility/non-associativity, secure
them, and combine this with the sequential executions. These algorithmic
strategies are reinforced with programmability suggestions to assure
deterministic executions. Finally, we verify these approaches on two modern HPC
systems: both versions deliver reproducible number of iterations, residuals,
direct errors, and vector-solutions for the overhead of less than 37.7 % on 768
cores.
