## 4 Parallel Version of the Improved Algorithm and Its Software Implementation
### Taking into Account the Capabilities of the Standard MPI-3
```
Data:\(A_{part}\), \(b_{part}\), \(x_{part}\equiv x_{part}^{(1)}\), comm_row, comm_col, \(N\) Result:\(x_{part}\) \(s\gets 1\) \(p_{part}\gets 0\) whileTrue do if\(s=1\)then \(Ax_{part\_temp}\gets A_{part}x_{part}\) \(Ax_{part}\leftarrow\) comm_row.Allreduce(\(Ax_{part\_temp}\)) \(b_{part}\gets Ax_{part}-b_{part}\) \(r_{part\_temp}\gets A_{part}^{T}b_{part}\) \(r_{part}\leftarrow\) comm_col.Allreduce(\(r_{part\_temp}\)) \(\sigma_{r\_part}^{2}\gets 0\) else \(r_{part}\gets r_{part}-\frac{q_{part}}{\text{scalar\_product\_pq}}\) \(\sigma_{r\_part}^{2}\leftarrow\sigma_{r\_part}^{2}+\frac{q_{part}^{2}}{(\text{ scalar\_product\_pq})^{2}}\) end scalar_product_temp\leftarrow(r_{part},r_{part}) \(\text{scalar\_product\_rr}\leftarrow\) comm_row.Allreduce(scalar_product_temp) \(\text{criterion\_temp}\leftarrow\frac{\Delta^{2}\sum\limits_{n}(\sigma_{r\_ part}^{2})_{n}}{\text{scalar\_product\_rr}}\) criterion\(\leftarrow\) comm_row.Allreduce(criterion_temp) if\(\text{criterion\ }\geqslant 1\)then return\(x_{part}\) end \(p_{part}\gets p_{part}+\frac{r_{part}}{\text{scalar\_product\_rr}}\) \(Ap_{part\_temp}\gets A_{part}p_{part}\) \(Ap_{part\_temp}\gets A_{part}^{T}(Ap_{part})\) \(q_{part\_temp}\leftarrow\) comm_col.Allreduce(\(q_{part\_temp}\)) \(\text{scalar\_product\_temp}\leftarrow(p_{part},q_{part})\) \(\text{scalar\_product\_pq}\leftarrow\) comm_row.Allreduce(scalar_product_temp) \(x_{part}\gets x_{part}-\frac{p_{part}}{\text{scalar\_product\_pq}}\) \(s\gets s+1\)
```
**Algorithm 3**Pseudocode for a (\"naive\") parallel version of the improved algorithm
```
defconjugate_gradient_method(A_part,b_part,x_part,comm_row,comm_col,N): N_part=size(x_part);M_part=size(b_part) x_classic_part=None delta=info(float64).eps S=1
```
**Algorithm 4**Pseudocode for a (\"naive\") parallel version of the improved algorithm
```P_part = zeros(N_part,dtype=float64) while True :  if s == 1 :  Ax_part_temp = dot(A_part, x_part)  Ax_part = empty(M_part,dtype=float64)  com_row.Allreduce([Ax_part_temp, M_part, MPI.DOUBLE],  [Ax_part, M_part, MPI.DOUBLE],  op=MPI.SUM)  b_part = Ax_part - b_part  r_part_temp = dot(A_part.T, b_part)  r_part = empty(N_part,dtype=float64)  comm_col.Allreduce([r_part_temp, N_part, MPI.DOUBLE],  [r_part, N_part, MPI.DOUBLE],  op=MPI.SUM)  sigma2_r_part = zeros(N_part,dtype=float64)  else :  r_part = r_part - q_part/scalar_product_pq  sigma2_r_part = sigma2_r_part + \  q_part**2/scalar_product_pq**2  scalar_product_temp = array(dot(r_part, r_part),  dtype=float64)  scalar_product_rr = array(0,dtype=float64)  comm_row.Allreduce([scalar_product_temp, 1, MPI.DOUBLE],  [scalar_product_rr, 1, MPI.DOUBLE],  op=MPI.SUM) criterion_temp = array(delta**2*sum(sigma2_r_part)/  scalar_product_rr,dtype=float64) criterion = array(0,dtype=float64)  comm_row.Allreduce([criterion_temp, 1, MPI.DOUBLE],  [criterion, 1, MPI.DOUBLE],  op=MPI.SUM) if criterion >= 1 :  return x_part, s, x_classic_part  p_part = p_part + r_part/scalar_product_rr  Ap_part_temp = dot(A_part, p_part)  Ap_part = empty(M_part,dtype=float64)  comm_row.Allreduce([Ap_part_temp, M_part, MPI.DOUBLE],  [Ap_part, M_part, MPI.DOUBLE],  op=MPI.SUM)  q_part_temp = dot(A_part.T, Ap_part)  q_part = empty(N_part,dtype=float64)  comm_col.Allreduce([q_part_temp, N_part, MPI.DOUBLE],  [q_part, N_part, MPI.DOUBLE],  op=MPI.SUM)  scalar_product_temp = array(dot(p_part, q_part),  dtype=float64)  scalar_product_pq = array(0,dtype=float64)  comm_row.Allreduce([scalar_product_temp, 1, MPI.DOUBLE],  [scalar_product_pq, 1, MPI.DOUBLE],  op=MPI.SUM)  x_part = x_part - p_part/scalar_product_pqs = s + 1 if s == N + 1 : x_classic_part = x_part.copy() ```
**Algorithm 4**Pseudocode for an efficient parallel version of the improved algorithm
``` Data:\(A_{part}\), \(p_{part}\), \(x_{part}\equiv x_{part}^{(1)}\), comm_row, comm_col, \(N\) Result:\(x_{part}\) \(s\gets 1\) \(p_{part}\gets 0\) whileTruedo if\(s=1\)then \(Ax_{part\_temp}\gets A_{part}x_{part}\) \(Ax_{part}\leftarrow\) comm_row.Allreduce(\(Ax_{part\_temp}\)) \(b_{part}\gets Ax_{part}-b_{part}\) \(r_{part\_temp}\gets A_{part}^{T}b_{part}\) \(r_{part}\leftarrow\) comm_col.Allreduce(\(r_{part\_temp}\)) \(\sigma_{r}^{2}\gets 0\) else \(r_{part}\gets r_{part}-\frac{q_{part}}{\text{scalar\_product\_pq}}\)  end if scalar_product_temp \(\leftarrow\) (r_{part},r_{part}\)  scalar_product_rr \(\leftarrow\) comm_row.Allreduce(scalar_product_temp) \(p_{part}\gets p_{part}+\frac{r_{part}}{\text{scalar\_product\_rr}}\) \(Ap_{part\_temp}\gets A_{part}p_{part}\) \(Ap_{part}\leftarrow\) comm_row.Allreduce(\(Ap_{part\_temp}\)) if\(s\geqslant 2\)then \(\sigma_{r}^{2}\leftarrow\sigma_{part}^{2}+\frac{q_{part}{}^{\odot 2}}{(\text{ scalar\_product\_pq})^{2}}\)  end if criterion_temp \(\leftarrow\)\(\frac{\Delta^{2}\sum\limits_{n}(\sigma_{r}^{2}\,\sigma_{r}^{2}\,)_{n}}{\text{ scalar\_product\_rr}}\) criterion \(\leftarrow\) comm_row.Allreduce(criterion_temp)  Wait(\(Ap_{part}\)) \(q_{part\_temp}\gets A_{part}^{T}(Ap_{part})\) \(q_{part}\leftarrow\) comm_col.Allreduce(\(q_{part\_temp}\))  scalar_product_temp \(\leftarrow\) (p_{part},q_{part}\)  scalar_product_pq \(\leftarrow\) comm_row.Allreduce(scalar_product_temp) Wait(criterion) if criterion \(\geqslant 1\)then return x_part  end if \(x_{part}\gets x_{part}-\frac{p_{part}}{\text{scalar\_product\_pq}}\) \(s\gets s+1\)  end if ```
**Algorithm 5**Pseudocode for an efficient parallel version of the improved algorithm
``` Data:\(A



```
defconjugate_gradient_method(A_part,b_part,x_part,comm_row,comm_col,N): N_part=size(x_part);M_part=size(b_part) x_classic_part=None delta=info(float64).eps requests=[MPI.Request()foriinrange(2)] s=1 P_part=zeros(N_part,dtype=float64) whileTrue: ifs==1: Ax_part_temp=dot(A_part,x_part) Ax_part=empty(M_part,dtype=float64) comm_row.Allreduce([Ax_part_temp,M_part,MPI.DOUBLE],[Ax_part,M_part,MPI.DOUBLE],op=MPI.SUM) b_part=Ax_part-b_part r_part=dot(A_part.T,b_part) r_part=empty(N_part,dtype=float64) comm_col.Allreduce([r_part_temp,N_part,MPI.DOUBLE],op=MPI.SUM) sigma2_r_part=zeros(N_part,dtype=float64) else: r_part=r_part-q_part/scalar_product_pq scalar_product_temp=array(dot(r_part,r_part),dtype=float64) scalar_product_rr=array(0,dtype=float64) comm_row.Allreduce([scalar_product_temp,1,MPI.DOUBLE],[scalar_product_rr,1,MPI.DOUBLE],op=MPI.SUM) Ap_part_temp=dot(A_part,P_part) Ap_part=empty(M_part,dtype=float64) requests[0]=comm_row.Allreduce([Ap_part_temp,M_part,MPI.DOUBLE],[Ap_part,M_part,MPI.DOUBLE],op=MPI.SUM) ifs>=2: sigma2_r_part=sigma2_r_part+\ q_part**2/scalar_product_pq**2 criterion_temp=array(delta**2*sum(sigma2_r_part)/ scalar_product_rr,dtype=float64) criterion=array(0,dtype=float64) requests[1]=comm_row.Allreduce([criterion_temp,1,MPI.DOUBLE],[criterion,1,MPI.DOUBLE],op=MPI.SUM)
```

**Listing 4** The Python code for an efficient parallel version of the improved algorithm

```
defconjugate_gradient_method(A_part,b_part,x_part,comm_row,comm_col,N): N_part=size(x_part);M_part=size(b_part) x_classic_part=None delta=info(float64).eps requests=[MPI.Request()foriinrange(2)] s=1 P_part=zeros(N_part,dtype=float64) whileTrue: ifs==1: Ax_part_temp=dot(A_part,x_part) Ax_part=empty(M_part,dtype=float64) comm_row.Allreduce([Ax_part_temp,M_part,MPI.DOUBLE],op=MPI.SUM) b_part=Ax_part-b_part r_part=dot(A_part.T,b_part) r_part=empty(N_part,dtype=float64) comm_col.Allreduce([r_part_temp,N_part,MPI.DOUBLE],op=MPI.SUM) sigma2_r_part=zeros(N_part,dtype=float64) else: r_part=r_part-q_part/scalar_product_rp q_part=q_part/scalar_product_rp q_part=q_part/scalar_product_rp q_part=empty(M_part,dtype=float64) request[0]=comm_row.Allreduce([x_part_temp,M_part,MPI.DOUBLE],op=MPI.SUM)
```

**Listing 5** The Python code for the efficient parallel version of the improved algorithm

```
defconjugate_gradient_method(A_part,b_part,x_part,comm_row,comm_col,N): N_part=size(x_part);M_part=size(b_part) x_classic_part=None delta=info(float64).eps requests=[MPI.Request()foriinrange(2)] s=1 P_part=zeros(N_part,dtype=float64) whileTrue: ifs==1: Ax_part_temp=dot(A_part,x_part) Ax_part=empty(M_part,dtype=float64) comm_row.Allreduce([Ax_part_temp,M_part,MPI.DOUBLE],op=MPI.SUM) b_part=Ax_part-b_part r_part=empty(M_part,dtype=float64) r_part=empty(N_part,dtype=float64) comm_col.Allreduce([r_part_temp,N_part,MPI.DOUBLE],op=MPI.SUM) sigma2_r_part=zeros(N_part,dtypeMPI.Request.Wait(requests[0],status=None) q_part_temp = dot(A_part.T,Ap_part) q_part = empty(N_part,dtype=float64) comm_col.Allreduce([q_part_temp,N_part,MPI.DOUBLE], [q_part,N_part,MPI.DOUBLE], op=MPI.SUM)  
scalar_product_temp = array(dot(p_part,q_part), dtype=float64) scalar_product_pq = array(0,dtype=float64) comm_row.Allreduce([scalar_product_temp,1,MPI.DOUBLE], [scalar_product_pq,1,MPI.DOUBLE], op=MPI.SUM)  
MPI.Request.Wait(requests[1],status=None) if criterion >= 1 : returnx_part,s,x_classic_part  
x_part = x_part - p_part/scalar_product_pq  
s = s + 1 ifs == N + 1 : x_classic_part = x_part.copy()