Conjugate Gradient Solvers with High Accuracy and Bit-wise Reproducibility between CPU and GPU using Ozaki scheme

The Ozaki scheme is the error-free transformation of dot product/matrix multiplication. This subsection presents a brief overview of the scheme. For further details of the Ozaki scheme, see the original paper [17].

Here, we explain the case of a dot product, but this scheme can be naturally extended to any dot product-based operations such as matrix-vector multiplication and matrix-matrix multiplication. Algorithm 1 shows the entire procedure of the Ozaki scheme for the dot product of two vectors and , where is the set of floating-point numbers (in this study, FP64). Briefly, this method consists of the following three steps:

1. Element-wise splitting of the input vectors into several split vectors
   
2. Computation of the all-to-all products of those split vectors
   
3. Element-wise summation (reduction) of the above all-to-all-product results
   

This method can be understood as an extension of high-precision arithmetic with multiple components (e.g., double-double arithmetic [7]) into the vector level. Specifically, first, the input vectors are split element-wise into the summation of several vectors using Algorithm 2 as

(1)

(2)

1. If ${{\boldsymbol {x}_{\tt split}}^{(p)}}_i$ and ${{\boldsymbol {y}_{\tt split}}^{(q)}}_j$ are non-zero elements,
   $\left|{{\boldsymbol {x}_{\tt split}}^{(p)}}_i\right| \ge \left|{{\boldsymbol {x}_{\tt split}}^{(p+1)}}_i\right|$ and $\left|{{\boldsymbol {y}_{\tt split}}^{(q)}}_j\right| \ge \left|{{\boldsymbol {y}_{\tt split}}^{(q+1)}}_j\right|$.
   
2. $\left({\boldsymbol {x}_{\tt split}}^{(p)}\right)^{T}{\boldsymbol {y}_{\tt split}}^{(q)}={\tt fl} \left(\left({\boldsymbol {x}_{\tt split}}^{(p)}\right)^{T}{\boldsymbol {y}_{\tt split}}^{(q)}\right),\\ 1\le p\le s_x, 1\le q\le s_y$.
   

The former implies that the accuracy of the final result can be controlled by omitting some lower split vectors. The key point of the latter is that the inner products of the split vectors can be computed with the standard floating-point arithmetic: for FP64 data, the DDOT routine provided in BLAS such as MKL and cuBLAS can be used⁹. Since ${\tt fl} \left(\left({\boldsymbol {x}_{\tt split}}^{(p)}\right)^{T}{\boldsymbol {y}_{\tt split}}^{(q)}\right)$ has no round-off error, that is, it is error-free and reproducible, even if it is computed with a non-reproducible operation.

Subsequently, the accurate and reproducible result is obtained by the summation of the all-to-all inner products of the split vectors. In this study, we compute the summation by a correctly rounded method, NearSum [21], to observe the result obtained by completely eliminating the rounding-error introduced in inner products. Although the summation can also be computed using working precision floating-point operations, the following points must care:

• The summation must be computed using a reproducible method. Since the summation is performed element-wise, fixing the computational order is neither difficult nor costly.
  
• log2 in Algorithm 2 must be computed by a reproducible method on different platforms, because the accuracy is not standardized in IEEE and may differ on different platforms (e.g., the accuracy is different between x86 and NVIDIA GPUs). However, this is not a concern when using a correctly rounded summation.
  

The CG method solves Ax = b where A is a symmetric positive definite matrix. Algorithm 3 shows the typical algorithm and corresponding BLAS routines for each linear algebra computation. Among them, the factor that may disturb the reproducibility of computed solutions on many core processors is the operations that consist of the inner product operations, namely sparse matrix-vector multiplication (SpMV), dot product (DOT), and 2-norm (NRM2). In this study, these operations are computed using the Ozaki scheme. SpMV and DOT achieve the correctly rounded results. The NRM2 is implemented using DOT as $r=\sqrt {{\tt DOT}(\boldsymbol {x}, \boldsymbol {x})}$, and the square root is performed on the standard FP64 operation.

In AXPY, we need to ensure consistency in whether or not the FMA operation is used. In this study, we use the AXPY implementation that uses the FMA operation. In addition, to prevent fast and less accurate computations for mathematical functions (e.g., -fp-model fast on ICC), we disable any less accurate options.

All computations are implemented using the standard FP64 arithmetic and the standard FP64 BLAS routines. Each linear algebra operation is performed through the corresponding BLAS routine shown in Algorithm 3 by using Intel MKL on CPUs and NVIDIA cuSparse (for SpMV) and cuBLAS (for the others) on GPUs. On the SpMV routines, the symmetrical structure of matrices is not considered (i.e., symmetric matrices are given to the computation after being expanded into general matrices). The coefficient sparse matrix is stored using the common Compressed Sparse Row (CSR) format. We note here that the choice of format does not affect the performance discussion in this study, as this study aims to examine the relative difference. In the GPU implementations, the BLAS routines are performed on GPUs, whereas the scalar value computations are performed on CPUs.

The Ozaki scheme is installed into DOT, NRM2, and CSRMV in the aforementioned FP64 implementation. In the Ozaki scheme, the internal computations are performed using MKL on CPUs and cuSparse and cuBLAS on GPUs. The splitting and summation portions are parallelized using OpenMP on CPUs and CUDA on GPUs. In Algorithm 2, to perform lines 9 and 10 correctly, the order of expression evaluation must be honored by using the compiler option “-fprotect-parens” on CPUs and the intrinsic for arithmetic on GPUs. 2^τ in lines 7 and 8 is computed using NextPowTwo [18]. AXPY is implemented using FMA.

In addition, we employ several techniques for speedup.

1. For SpMV, as the coefficient matrix is not changed during the iterations, the splitting of the matrix is needed only once before the iteration starts. This contributes to a nontrivial performance increase, because matrix splitting becomes the major cost in the Ozaki scheme on memory-bound operations.
   
2. The inner products of the split vectors can be performed using a dense matrix multiplication (GEMM) by combining multiple split vectors into a single matrix, as shown in Figure 1. This strategy contributes to the speedup of memory-bound operations in the Ozaki scheme by reducing memory access. The same concept can be applied to SpMV: the computation can be performed using a sparse matrix - dense matrix multiplication routine (SpMM), which is available on MKL as mkl_dcsrmm and on cuSparse as cusparseDcsrmm. On CPUs, however, we do not use this technique because performance is degraded with mkl_dcsrmm¹⁰.
   
3. In SpMV, we use asymmetric splitting [18]. σ at line 8 in Algorithm 2 determines how many bits are stored in each element of the split matrices/vectors; smaller σ increases the number of bits that can be held. σ is determined not to cause an overflow in the computation of the product of the split matrices and vectors but is chosen to be as small as possible in the powers of 2. We can bias the ρ at line 2 to reduce the number of split data on either the matrix or vector by increasing ρ on one side and decreasing ρ on another side by the same amount. If SpMM is used in the computation, the reduction of the number of split matrices increases the chance of a speedup in general. Our GPU implementation decreases the ρ on the matrix side by the minimum amount that decreases the number of split matrices. On the other hand, the CPU implementation, which does not use SpMM in its computation, increases the ρ on the matrix side by the maximum amount that does not change the number of split matrices. This strategy increases the chance of reduction in the number of split vectors. The optimal ρ is determined by trial and error by performing the matrix splitting with the ρ reduced step by step. As this determination is performed only once before starting the iterations of the CG method, the cost is not high.
   

Table 2: (a) Relative true residual and (b) number of iterations on different platforms. Note: in all cases, true residual was computed on the correctly rounded operations using the Ozaki scheme. FP64Oz-CR got identical result on all platforms.

Table 2 shows both the relative true residual (||b − Ax_i||/||b||) when ||r_i||/||b|| ≤ 10^{− 16} and the number of iterations across four platforms. In most cases, the results (both the solution and the number of iterations) of FP64 are non-identical among different platforms. However, the results of FP64Oz-CR across the four platforms were identical; that is, reproducibility was ensured by the Ozaki scheme. In addition, some cases converged with fewer iterations by FP64Oz-CR with accurate computation when compared with FP64 implementations.

Figure 2 shows the residual plots for every 10 iterations on four platforms. Solid lines show the relative residual ||r_i||/||b|| and dotted lines show the relative true residual ||b − Ax_i||/||b||¹². Significant differences can be observed in the convergence plots among the five cases, but in all cases, the solution converges to a similar value on the same order (however, most of them are non-identical, as shown in Table 2). In terms of the residual r_i, FP64Oz-CR often converged with fewer iterations than FP64, but in those cases, the stopping criterion of 10^{− 16} was too small.

Before presenting the experimental results, we discuss the expected performance. When SpMV with the Ozaki scheme is computed using SpMM (i.e., our GPU implementation), the execution time overhead of FP64Oz-CR against FP64 per iteration can be roughly estimated. Ideally, if the matrix is sufficiently dense, the overhead becomes close to d times, where d is the number of split matrices, which can be obtained by performing the splitting of the coefficient matrix once. When SpMV is computed using SpMM, it requires 4d times execution time overhead against the standard floating-point operation in terms of memory access to the matrix (we assume that the cost to the vector can be ignored as it is small enough when compared with the matrix). The 3/4 of the 4d times overhead arises in the splitting process (the accesses to vector x at lines 9 and 10 in Algorithm 3); however, it is eliminated on CG methods since the splitting is performed only once before iteration. In CG methods, the SpMV cost is usually dominant, as matrix-vector multiplication is an O(n²) operation, whereas the others are O(n). Therefore, only the d times overhead appears. However, if the matrix is highly sparse, and the execution efficiency becomes low owing to the frequent random memory access, the cost of other operations than SpMV becomes nonnegligible, and the overhead becomes unpredictable. As DOT requires 4d times overhead (assuming cache hits all split vectors), the cost of DOT (and NRM2) may become dominant instead of SpMV. Moreover, the value of d is unpredictable in CG methods because the vectors are updated during iteration. In addition, as d increases, the chance of cache-missing increases.

Tables 3 and 4 show the number of split matrices required to achieve correct-rounding (d), the execution time until convergence, and the execution time overhead (times) of FP64Oz-CR against FP64 on CPUs and GPUs, respectively. Note that the required number of splits becomes minus-one on GPUs when compared with that on CPUs as we use the asymmetric-splitting technique described in Section 4. The FP64 implementations, which are the baselines for the comparison, can be seen as the ideal performance as they are constructed on vendor-implemented routines. As the CPU implementation does not use SpMM, we consider that the GPU implementation shows more desirable results. The number of split matrices corresponds to the expected minimum overhead, as explained before. When compared with the number of split matrices, approximately 1.3 – 7.2 times additional overhead was observed on GPUs, and we can see that the smaller nnz/n it has, the greater the overhead tends to take, following the previous discussion.

Figure 3 shows the breakdown of the execution time of 100 iterations on the tmt_sym and nd24k matrices on CPU1 and GPU1. Both matrices have the smallest and biggest sparsities, respectively. tmt_sym shows a high cost for level-1 BLAS operations, while nd24k shows a high cost for SpMV (SpMM). FP64Oz-dn shows the result executed with a specified d (the number of split matrices/vectors). When d is specified, the asymmetric-splitting technique is not used; only with FP64Oz-CR does the matrix splitting cost (SplitMat) include the tuning cost for the asymmetric splitting. However, it is not so large within 100 iterations.

Iakymchuk et al. have already proposed accurate and reproducible CG solvers [8, 9] based on the ExBLAS approach [10]. These CG solvers are parallelized with the flat MPI as well as MPI and OpenMP tasks but support only CPUs. We evaluate their performance and compare it against that of the proposed method. As with the proposed method, the ExBLAS-based implementation ensures correct rounding in all dot product operations in CG. The ExBLAS approach efficiently combines the Kulisch long accumulator [11] and floating-point expansions (FPEs). While the long accumulator is robust and designed for severe (ill-conditioned) cases, keeping every bit of information until the final rounding, FPEs are unevaluated sums (arrays of FP64) to target a limited range of numbers (e.g., nonsevere dynamic ranges and/or condition numbers). Hence, ExBLAS aims to use FPEs as much as possible owing to their speed and only occasionally long accumulators (e.g., when the accuracy of FPEs is insufficient, or at the final rounding to FP64). One clear advantage of the ExBLAS-based implementation against our proposed approach is its low memory consumption: it uses only 2097 bits for a long accumulator and 192-512 bits for an FPE per MPI process, while the proposed method consumes a large amount of memory for storing split matrices (i.e., in proportion to the number of split data d).

The ExBLAS-based implementations have two versions: the MPI-OpenMP hybrid parallel [8] and the flat MPI version [9]. However, the flat MPI version was faster than the hybrid version on both CPU1 and CPU2. Therefore, the following evaluation was conducted using only the flat MPI version. We conducted these experiments with the same conditions as the other evaluations except for the following points: The code¹³ was compiled using GCC 8.3.1 on CPU1 and GCC 7.5.0 on CPU2 with -mavx -fabi-version=0 -fopenmp. On both CPUs, we executed the code by mapping one MPI process per core. The ExBLAS-based implementation supports preconditioning, but it was disabled in this evaluation.

Table 5 shows the results obtained by the ExBLAS-based implementation (FP64Ex-CR) on CPU1 and CPU2. We note that the results of FP64Oz-CR and FP64Ex-CR may differ because the implementations of the CG algorithm itself are different and the strategies for reproducibility are clearly not the same. Since FP64Oz-CR and FP64Ex-CR adopt different strategies for accurate and reproducible computations, their overhead compared to the baseline FP64-based implementation depends on both the matrix at hand and the processors used. The FP64Ex-CR achieved better (lower) overhead on five out of eight matrices on CPU1. On the other hand, the proposed method FP64Oz-CR performs better on six out of eight matrices on CPU2. Even though FP64Ex-CR does not support GPUs, the best performance of FP64Oz-CR is on GPUs, where the potential of FP64Oz-CR flourishes with the use of SpMM and the reduction of the number of split matrices.