## 2 Sequential Algorithm and Its Software Implementation
### Improved Implementation of the Conjugate Gradient Method  
In work [23], a modification of the above conjugate gradient method was proposed, which takes into account rounding errors accumulated during of calculations when making a decision to terminate the iterative process. We set \(p^{(0)}=0\), \(s=1\), and an arbitrary initial approximation \(x^{(1)}\). Then we repeatedly perform the following sequence of actions:  
\[r^{(s)} =\begin{cases}A^{T}\big{(}A\,x^{(s)}\,-\,b\big{)},&\text{if}\quad s =1,\\ r^{(s-1)}\,-\,\frac{q^{(s-1)}}{\big{(}p^{(s-1)},q^{(s-1)}\big{)}},&\text{if} \quad s\geqslant 2,\\ \sigma_{r^{(s)}}^{2}= \begin{cases}\big{(}A^{T}\big{)}^{\circ 2}\Big{(}A^{\circ 2 }x^{\circ 2}+b^{\circ 2}\Big{)},&\text{if}\quad s=1,\\ \sigma_{r^{(s-1)}}^{2}+\frac{\big{(}q^{(s-1)}\big{)}^{\circ 2}}{\big{(}p^{(s-1)},q^{(s-1 )}\big{)}^{2}},&\text{if}\quad s\geqslant 2,\\ \end{cases} \] \[\frac{\Delta^{2}\,\sum\limits_{n=1}^{N}\big{(}\sigma_{r^{(s)}}^{2} \big{)}_{n}}{\|r^{(s)}\|_{2}^{2}}\geqslant 1,\]  
\[p^{(s)} =p^{(s-1)}\,+\,\frac{r^{(s)}}{\big{(}r^{(s)},r^{(s)}\big{)}},\] \[q^{(s)} =A^{T}\big{(}A\,p^{(s)}\big{)},\] \[x^{(s+1)} =x^{(s)}\,-\,\frac{p^{(s)}}{\big{(}p^{(s)},q^{(s)}\big{)}},\] \[s =s+1.\]  
Here, \({}^{\circ 2}\) is Hadamard power, that is, element-wise raising of a vector/matrix to the second power, \(\Delta\)--relative rounding error. For half-precision calculation (float16) \(\Delta=10^{-3.3}\), for single precision calculation (float32) \(\Delta=10^{-7.6}\), for calculation with double precision (float64) \(\Delta=10^{-16.3}\), and for calculating with quadruple precision \(\Delta=10^{-34.0}\). This algorithm can be written as the following pseudocode (see Algorithm 2). In this algorithm, the lines highlighted in red correspond to the actions that must be performed in order to implement the improved criterion for terminating the iterative process. If these lines are removed and the condition \(True\) is changed to \(s\leqslant N\), then we get the \"classical\" implementation of the conjugate gradient method for solving System (1). The Python code for the function that implements the conjugate gradient method for solving System (1) with the improved iterative process termination criterion is presented at Listing 2. This software implementation contains the following features:  
1. The function returns (1) the array x that contains the solution of System (1) found by the improved Algorithm 2, (2) the number of iterations \(s\) that needed to be done by the algorithm to find an approximate solution, (3) the array x_classic, which contains the solution found using the \"classical\" Algorithm 1 (this array contains such a solution only if the number of iterations performed by the algorithm \(s\geqslant N+1\)). 2. This implementation (and the corresponding pseudocode) takes into account that the result of each dot product is used twice, so the result of the dot product is stored in a separate variable. 3. At the first iteration of the algorithm (\(s=1\)), the value of \(\sigma_{r}^{2}\) is calculated using the command sigma2_r = dot(A.T**2, dot(A**2, x**2) + b**2) A feature of Python is that under the arrays A**2 and A.T**2, which contain the matrices \(A^{\circ 2}\) and \(\big{(}A^{T}\big{)}^{\circ 2}\), a separate memory space is allocated. This is very bad, because it is assumed that problem (1) is being solved with a huge matrix \(A\). Therefore, when solving real \"large\" problems, an array containing matrix \(A\) can take up a large part of the memory. As a result, there may not be enough space in memory for additional arrays. This problem can be worked around in several ways. In the first case, it is possible to issue the specified command as a separate function in the C/C++/Fortran programming language, using for calculations only the elements of the array in which matrix \(A\) is stored.