For solving systems of linear algebraic equations with blockvediagonal matrices arising in geoelectrics and di usion problems, the parallel matrix square root method, conjugate gradient method with preconditioner, conjugate gradient method with regularization, and parallel matrix sweep algorithm are proposed and some of them are implemented numerically on multi-core CPU Intel. Investigation of e ciency and optimization of parallel algorithms for solving the problem with quasi-model data are performed. The problem with quasi-model data is solved.
Systems of linear algebraic equations (SLAE) with block- vediagonal matrices arise when solving some mathematical modelling problems, in particular, geoelectrics and di usion problems. The geoelectrics problems are very important for investigating the crust heterogeneity. This problem is described by a differential Poissons equation. After using a nite-di erence approximation, the two- and three-dimensional lateral logging problems are reduced to solving illconditioned systems having large-scale block- ve diagonal matrix [1]. Another important problem is the multicomponent di usion problem when it is necessary to know concentration distribution of di using components at every moment of time. This problem is also reduced to solving a SLAE with block- vediagonal matrix after approximating systems of di erential equations [2]. These types of problems has the matrix with a special structure, we also describe.
In this paper, for solving SLAE with block- vediagonal matrices, direct and iterative parallel algorithms are designed: conjugate gradient method with preconditioner (PCGM), conjugate gradient method with regularization (PCGR), square root method (PSRM), and parallel matrix sweep algorithm (PMSA). ? This work was partially supported by Program of UrB RAS (project no. 15-7-1-13) and by the Russian Foundation for Basic Research (project no. 15-01-00629 a).
The PCGM, PCGR, and PSRM algorithms are implemented on multi-core Intel processor. The parallel matrix sweep algorithm will be implemented on multi-core Intel processor later.
Investigation of e ciency and optimization of parallel algorithms are performed. The problem with quasi-model data is solved. 2
Parallel Algorithms for Solving SLAE
>8 C0Y0 D0Y1 + E0Y2 = F0; >><> B1Y0 + C1Y1 D1Y2 + E1Y3 = F1;
AiYi 2 BiYi 1 + CiYi DiYi+1 + EiYi+2 = Fi; i = 2; :::; N >>> AN YN 2 BN YN 1 + CN YN DN YN+1 = FN ; >: AN+1YN 1 BN+1YN + CN+1YN+1 = FN+1; where Yi are unknown n{vectors, Fi are given n{vectors, Ai; Bi; Ci; Di; Ei are given n n { matrices.
As mentioned in introduction, after a nite-di erence approximation, the geoelectrics problems and di usion problems can be reduced to solving SLAE with block- vediagonal matrices with the structure presented in Fig. 1. This special case arising when we are solving di usion problems for higher order scheme and geoelectric problems.
Figure 1 show the special case of SLAE (
For solving SLAE with block- vediagonal matrices, parallel algorithms based on conjugate gradient method with preconditioner (PCGM), conjugate gradient method with regularization (PCGR), parallel square root method (PSRM), and parallel matrix sweep algorithm (PMSA) are proposed. 2.1
Parallel Conjugate Gradient Method with Regularization One of the fast iterative algorithms for solving a SLAE with a symmetric positive de nite matrix is the conjugate gradient method (CGM) [3].
We write the system (
Ax = b;
When we are solve geoelectric problems with non-uniformly scaled coe cients and when the system is ill-conditioned, we consider a regularization to prevent e ect of rounding errors. We add regularization with alpha parameter of regularization to ensure stability of the method with the following equation: Aex = b; Ae = A +
E; where is the parameter of regularization.
The CGM method has the following form: One of the fast iterative algorithms for solving a SLAE with a symmetric positive de nite matrix is the conjugate gradient method (CGM). The introduction of a preconditioner is applied to accelerate the convergence of the iterative process. Preconditioning consists in the fact that the initial system of equations Ax = b is replaced by the system
C 1Ax = C 1b;
for which the iterative method converges essentially faster. The condition for the
choice of the preconditioner C is the following:
cond(A~) << cond(A); cond(A~) = ~max ; cond(A) =
~min
max ;
min
where cond(A) and cond(A~) are the condition numbers of the matrices A and
A~; max; ~max and min; ~min are the largest and smallest eigenvalues of the
matrices A and A~, respectively.
(
For system of equations (
The condition for stopping the CGM iterative process with preconditioner is Azk b = kbk < ".
Parallelization of the PCGR and PCGM is based on splitting the matrix
A by horizontal lines into m blocks, and the solution vector z and the vector
of the right-hand part bof the SLAE are divided into m parts such that n =
m L, where n is the dimension of the system of equations, m is the number of
processors, and Lis the number of lines in the block (Fig. 2). At every iteration,
each processor computes its part of the solution vector. In the case of the
matrixvector multiplication (A by z), each processor multiplies its part of the lines of the
matrix A by the vector z. The host-CPU (master) is responsible for transferring
data and, also, calculates its part of the solution vector.
For solving SLAE (
ST y = b; Sz = y: (8) To solve systems (8), we use the recursion formulas (9)
The main idea of parallelization of the square root method for multiprocessor
computers with shared memory is based on parallel computing the elements
sij ;j = i; :::; n; of each i -th line of the matrix S : The i -th line is divided into m
parts so that n i = m Li, where i is the line number, n is the dimension of
the system, m is the number of processors, and Li is the number of line items
that are evaluated by each processor (Fig. 3) described in [4].
Let us consider system (
Yi
AiYi 2 BiYi 1 + CiYi DiYi+1 + EiYi+2 : Similarly to Yi, we de ne the operators Ui; Vi; Pi; Qi ; Wi .
A reduced system of equations with respect to the parameters YK; YK+1 is constructed as follows. In L subdomains de ned by the intervals (K; K +M +1), K = 0 ; M; :::; N M; we consider the following problems: 8 8 8 8 > 00 1 00 1 00 1 00 1 > > > >> Uin = 0; UKn = BB@:0:CCA ; UKn+1 = BB@:0:CCA ; UKn+M = BB@:0:CCA ; UKn+M+1 = BB@:0:CCA ; > > > > > :> 1 0 0 0 > > >> Vi1 = 0; VK1 = BB@:0:CAC ; VK1+1 = BB@:0:CAC ; VK1+M = BB@:0:CCA ; VK1+M+1 = BB@:0:CCA ; > > > > > >> 0 0 0 0 > <
: : : : : : : : : : : : : : > 00 1 00 1 00 1 00 1 > > > >> Vin = 0; VKn = BB@:0:CAC ; VKn+1 = BB@:0:CAC ; VKn+M = BB@:0:CCA ; VKn+M+1 = BB@:0:CCA ; > > > > > >: 0 1 0 0 > > >> Pi1 = 0; PK1 = BB@:0:CAC ; PK1 +1 = BB@:0:CAC ; PK1 +M = BB@:0:CCA ; PK1 +M+1 = BB@:0:CCA ; > > > > > >> 0 0 0 0 > <
: : : : : : : : : : : : : : > 00 1 00 1 00 1 00 1 > > > >> Pin = 0; PKn = BB@:0:CAC ; PKn+1 = BB@:0:CAC ; PKn+M = BB@:0:CCA ; PKn+M+1 = BB@:0:CCA ; > > > > > >: 0 0 1 0 > > >> Qi1 = 0; Q1K = BB@:0:CAC ; Q1K+1 = BB@:0:CAC ; Q1K+M = BB@:0:CCA ; Q1K+M+1 = BB@:0:CCA ; > > > > > >> 0 0 0 0 > < > 00 1 00 1 00 1 00 1 > > > >> Qin = 0; QnK = BB@:0:CAC ; QnK+1 = BB@:0:CAC ; QnK+M = BB@:0:CCA ; QnK+M+1 = BB@:0:CCA ; > > > > > >: 0 0 0 1 where i = K + 2; :::; K + M 1:
The following theorem is proved in [5]. (15) (16)
Theorem. If Ui1; :::; Uin are solutions of (10), Vi1; :::; Vin are solutions of (11),
Pi1; :::; Pin are solutions of (12), Qi1; :::; Qin are solutions of (13), Wi are solutions
of (14), and Yi are solutions of given problem (
Yi = Ui1Ui2:::Uin YK + Vi1Vi2:::Vin YK+1 + Pi1Pi2:::Pin YK+M + + Qi1Qi2:::Qin YK+M+1 + Wi :
Substituting relations (15) into the given system (
A~Y = F :
Reduced system (16) is one of vector equations with block-sevendiagonal matrices of coe cients; in each line, one of the seven vector elements YK being on the left or right from the main diagonal is zero. After nding vectorparametres YK ; YK+1, other required unknown vectors are expressed through vector-parameters and are found in each subdomain L independently by formula (16).
The direct parallel matrix sweep algorithm for solving vector system (
Problems (10){(14) and system of equations (16) can be solved by the matrix sweep algorithms (Gaussian elimination methods) for solving systems of equations with block- vediagonal and block-sevendiagonal matrices, respectively. The formulas of the matrix sweep algorithms for solving SLAEs with block- vediagonal and block-sevendiagonal matrices are deduced similarly to formulas of the corresponding scalar sweep algorithms [6].
The stable parallel matrix sweep algorithm for solving SLAE with block- vediagonal matrices can be e ectively implemented on parallel computing systems with distributed memory with L processors (the number is equal to the number of subdomains). Problems and unknown vector-parameters Yi ( i = K + 2; :::; K + M 1) inside of each subdomain L are computed on L processors independently. In addition, the parallel matrix sweep algorithm can be e ectively implemented on multi-core processor and graphic processors. 3
The parallel algorithms were implemented on a multi-core Intel processor using the technology of OpenMP, and the Windows API development tools.
With the help of the parallel preconditioned conjugate gradient method, and square root method, we solved the problem of nding a potential distribution in a conducting medium with known quasi-model solution.
The source data and quasi-model solution of the problem were provided by the Department of Borehole Geophysics, Institute of Geology and Geophysics, Siberian Branch of RAS (Novosibirsk).
After discretization the problem is reduced to solving a SLAE with an illconditioned symmetric positive de ned block- vediagonal matrix of dimension 134862 134862 with 247 square blocks.
The numerical solution of the problem is compared with the quasi-model solution by means of calculating the relative error =
Y M
Y N = Y M ; (18) where Y M is the quasi-model solution of the problem, Y N is the numerical solution of the problem.
The condition is chosen as a stopping criterion for the iterative PCGM.
A priori we nd the condition number of the original matrix: cond(A) = max min
This problem was solved by the parallel conjugate gradient method with preconditioner, parallel conjugate gradient method with regularization, and parallel square root method. The numerical solution of the problem coincides with the quasi-model solution with accuracy
For the quasi-model data the numerical solution of the problem is presented in Fig. 5.
The computation times of solving the SLAE in the potential distribution problem on the hybrid computing system are presented in Table 1. This system is installed in the Department of Ill-posed Problems of Analysis and Applications of the Institute of Mathematics and Mechanics UrB RAS. The computing system consists of 4-core processor Intel Core I5-750.
Note that the time for solving the problem by PCGM without preconditioner on one core Intel Core I5-750 for P CGM = 10 3 was 30 minutes.
Table 1 shows that the preconditioner decreases essentially the time of solving the problem.
At the beginning, for a multi-core processor with shared memory, the PCGM, PCGR, and PSRM parallelization was implemented by means of the operating system (OS) threads by the development tools Windows API. For parallel implementation of the computing block of the program, the threads were created. Each of these threads executed on the OS \logical processor" and computed the data portion. At the end of each computing block, the barrier synchronization of the threads was perfomed.
In order to optimize the program and reduce the computation time, the OpenMP technology was used. Automatic parallelization of loops was carried out by the OpenMP library using special compiler directives.
The interval of size L of the loop variable i was divided into m parts. Each thread of the process calculated its p-th part of the data, where p = L=m (Fig. 4).
So, the PSRM is the fastest method. The computation time for solving the SLAE on a 4-core CPU Intel is reduced to several seconds (in comparison with 30 min by CGM without preconditioner). 4
For solving SLAE with block- vediagonal matrices arising in geoelectrics and di usion problems, the parallel conjugate gradient method and parallel square root method with preconditioner are proposed and numerically implemented on a multi-core processor Intel. Investigation of e ciency and optimization of parallel algorithms for solving the problem with quasi-model data are performed.
The calculation results show that the use of parallel PCGR, PSRM, and PCGM with preconditioner allows us to solve rather e ectively problems with ill-conditioned matrices on multi-core CPU.
For the future work we want to compare parallel matrix sweep algorithm with others and to implement and optimize for GPU.