We have developed a library for fast, accurate, and reproducible fundamental linear algebra operations that are frequently used in scientific computations as the underlying layer. The library is called Exact BLAS (Basic Linear Algebra Subprograms). In Robust, our focus is on enhancing and expanding the ExBLAS library but even more on developing robust iterative linear solvers, which are frequently used in various scientific computations and consume significant amount of their execution time. Our first target was the Preconditioned Conjugate Gradient (PCG) method for solving symmetric sparse systems such as the ones from mashing surfaces of planes and cars. This time we have approached the problem from different perspective: instead of making all floating-point operations strictly reproducible, which can become expensive, we have decided to at first identify these parts in the PCG solver that are prone to non-reproducibility and lack of accuracy. While working on this task, we have gained precision knowledge on iterative solvers, architectures, and compiler optimization (e.g. random replacement of instructions in the favor of the fused multiply-add, fma, instruction). Therefore, we fix these issues by expanding and modifying the ExBLAS approach to tackle parallel computations of residual (dot product of two vectors), restricting compiler to the explicit user-defined usage of the fma instruction, etc. However, after analyzing the actual required precision of the PCG solver, we understand that there is a possibility to enhance our approach even further. Thus, we use a less expensive part of ExBLAS (a short vector of floating-point numbers to store both error and result) with additional optimization technique. Phenomenally, a vector with only three elements is enough to get accurate and reproducible result in practice due to very small and insignificant error that is propagated to the tail of this vector although this is difficult to prove using the rigorous and often pessimistic numerical analysis tools. This approach enhances the performance of robust PCG even further. Furthermore, we prove that both approaches deliver identical intermediate and final results in terms of residuals and number of iterations for various number of resources used but also cross-platform reproducibility, which is a rare case. We conduct our tests on synthetic matrices as well as on real cases from the SuiteSparse matrix collection. We extend this research to the Krylov Subspace method called Preconditioned BiCGStab that works on unsymmetric matrices as well as to its pipelined version, which means better suited for runs on large scale. With this, we provide reproducible version of iterative solvers that are capable to work on different types of matrices.
