The paper presents a parallel calculation scheme for distributed memory systems for a modification of the AMLI preconditioner that performs incomplete inversions of matrix blocks by a technique that produces approximate inverse of a matrix by an orthogonalization procedure, referred to as the incomplete basis matrix method. Theoretical estimates of algorithm's performance and the results of its experimental testing are presented. The results of experimental analysis of preconditioner parameters influence on the convergence of the iterative algorithm are given. The characteristics of the proposed algorithm are compared with those of known algorithms implemented in the hypre software package. The tests on the matrices obtained from finite element discretization of soil stress-strain state modelling problem showed that the proposed algorithm has better performance for significantly ill-conditioned linear systems when solving them on a small number of computational resources.
The need to solve linear algebraic systems arises in mathematical modelling of most physical processes, in particular when modelling soil stress-strain state. In many cases, linear systems have ill-conditioned sparse matrices of large size. In these situations, iterative methods, such as conjugate gradient method (CG) or biconjugate gradient stabilized method (BiCGstab), are most commonly used [1].
The main method used to improve the convergence and the accuracy while solving ill-conditioned linear systems is the use of preconditioners [1, 2] - matrices, whose multiplication on linear system matrix leads to the decrease of condition number. The most commonly used methods for constructing preconditioners are incomplete matrix decomposition (e.g. ILU decomposition [1]), incomplete matrix inversion (e.g. polynomial preconditioners [3]), Geometric and Algebraic Multigrid (AMG) [4] methods, the Algebraic Multilevel Iteration method (AMLI) [5,6], and wavelet preconditioners [7,8].
The so-called basis matrix method [9, 10] is one of the methods developed to analyse ill-conditioned and rectangular linear systems. This method constructs iteratively the inverse of a given matrix using some orthogonalization procedure. Based on this method, an incomplete inversion algorithm is created, that can be used as a preconditioner. Further, it is combined [11] with the AMLI preconditioner and applied for solving elasticity problems.
The speed of the AMLI preconditioner, similarly to the multigrid ones, in particular the AMG, is low comparing to the incomplete decomposition and the incomplete inversion preconditioners and its application takes most of the time spent to obtain solutions of linear systems. This makes crucial the development of parallel schemes for these algorithms. While the AMLI is shown to have good scalability on vector and parallel-vector machines [12, 13], its parallel implementations on cluster systems are poorly studied. 2 Preconditioner on the base of the AMLI and the basis matrix method The AMLI preconditioner [5] is based on the representation of an input matrix as having 2x2 blocks with its subsequent LU factorization:
A11 A = A21
A12 I
= A22 A21 A111 A11 I
A12
S where S = A22 A21 A111 A12 is the Schur complement of A .
After approximating the Schur complement with a polynomial and the recursive
application of the expansion (
1
A of the size dimA = [K * L, K * L] that is divided into L * L blocks Aij , i, j = 1,..., L of the size dimAij = [K , K ] we have
A1(1k ) = ALk ,Lk , dimA1(1k ) = [K , K ],
A1(2k ) = ( ALk ,Lk 1,..., ALk ,L ), dimA1(2k ) = [K , kK ],
A2(1k ) = ( ALk 1,Lk ,..., AL,Lk ), dimA2(1k ) = [kK , K ],
(
AL,Lk 1 ... ... ...
ALk 1,L ... , dimA2(2k ) = [kK , kK ],
ALL dimS (k 1) = dimA2(2k ) , S (0) = ALL1, dimM (k ) = [(k 1)K , (k 1)K ], k = 1,..., L 1 and M (L1) is an approximate inverse of A .
The AMLI algorithm can be optimally stabilized by choosing the degree of an
approximating polynomial [12]. However, the applicability of this theory while solving
stress-strain problems described by hyperbolic differential equations wasn't proved.
As such problems are of our specific interest, we consider the simplest case of
Pv (t) = P1 (t) = 1 t . Then, the Schur complement approximation becomes
S (k 1) M (k 1)
and the recursive scheme (
I M (k ) = A2(1k ) A1(1k ) 1
xN = A1(1N ) wN , xi = wi A1(1i) A1(2i) ( xi1 ,..., xN ), i = 0,..., N 1.
The ``reverse`` part consists in the sequential computation of the resulting vector A1(1k ) 1
In the scheme (
Summarizing, the sequential computational scheme for the preconditioning
procedure (
11 A = A2(01)1 ... A2(01)N
A(0) 121 ... ... ...
...
...
A(N 1) 11 A2(1PN1)
A1(20)N ... A1(2NN1) ,
A(N ) 11
A1(2i) = A1(2i,)i1,..., A1(2i)N , A2(1i) = A2(1i),i1,..., A2(1i)N .
Considering the application of the preconditioner (
On the ``direct`` part of the procedure, values of the following vectors are
sequentially computed:
(
11 are the approximates to the corresponding inverse matrices obtained by the incomplete basis matrix method [16].
As the initial matrix is sparse, we use the compressed row storage (CSR) scheme for all the matrices.
The approximation accuracy is controlled on the iterations of the algorithm by calculating ii for all the processed rows. When 1 becomes bigger than a given threshold value, we perform complete inversion. As the resulting matrices aren't always dense, we also use the CSR storage scheme doing complete inversion. 3
Parallel implementation of the preconditioner
To make the AMLI algorithm combined with matrix blocks incomplete inversions
using the basis matrix method run in parallel on distributed memory systems, the
following data partitioning scheme can be applied. Each of i row blocks
A2(01i) ,..., A2(1ii1) , A(i) , A(i)
11 12,i1,..., A1(2i)N
is distributed on a system of P processes by blocks of rows. Thus, each process
stores blocks of rows of the submatrices A(i) , A(i) , A2(1i)j and can perform
multiplica11 12
tions that arise while applying the preconditioner (
The computational scheme that we use for multiplying a sparse matrix A of size M M by a vector v has the following form:
Bp = [rp0,rp1 ], iP=0Bi = [0,..., M ], Bi Bj = , i, j of the matrix A rows and the corresponding elements of the vector v ; 2. Denote a set of columns of the matrix A in which non-zero elements are present in the row block Bp as Cp ; 3. The process p sends values of the elements of the vector v from the index set Ci Bp to the processes i, Ci Bp and receives the values of the elements from the index set Cp Bj from the processes j, C p Bj ; 4. After the synchronization is performed on the step 3, each of the processes independently multiplies the block of rows of the matrix A by the vector v . The result of this multiplication has the same partitioning as the vector v has.
Performing the ``direct`` part (
The matrices A (i)
11 can be obtained in parallel the following way: The process p handles the sequence k p2 > k k p1 of the matrices A(k) 11 . The sequences are constructed to make the processes evenly loaded as far as it is possible. This is ensured by approximately the same number of non-zero elements in the matrices A1(1k) that are processed by each of the processes. The processes exchange the values of the elements of the matrices A1(1i) to make the matrices A1(1k) , k p2 > k k p1 completely stored in the memory of the process p . The process p sends row blocks Bp of the matrices A1(1k) , ki2 > k ki1 to the processes i, i p and receives row blocks Bi of the matrices A1(1k) , k p2 > k k p1 from them. Each of the processes independently performs a complete or incomplete inversion of the matrices. The processes exchange the values of the elements of the matrices A (i) 11 to restore the initial data partitioning. The process p sends row blocks Bi of the matrices A1(1k) , k p2 > k k p1 to the processes i, i p and receives row blocks Bp of
Since the data partitioning of the blocks of rows within the matrices A1(2i) , A2(1ji) may differ for the different values of i , an additional synchronization should be performed when vectors are multiplied by these matrices. Denote a block of rows from the submatrices A1(2i) , A2(1ji) that is handled by the process p as Bpi . Then, before multiplying a vector v with the partitioning Bi = (B0i ,..., BPi ) by the matrix A2(1i) , we need to synchronize its elements to make each of the processes store the block Bpi ...BpN . To do so, the process p sends the block Bpi Brj of the elements of the vector to the process r if Bpi Brj , j i . Then, the process p receives a corresponding block from the process r if Bpj Bri , j i .
Similar procedures should be performed when a vector v = (vi1,..., vN ) is multiplied by the matrix A(i)
12 for each of its components. When the non-changeable
component x j is multiplied by the matrices A1(2i) during the ``reverse`` part (
When the preconditioner is used with an iterative method for solving a linear alge braic system, we perform a calculation of the matrices A11 and the sets of the elements that need synchronization once on the first iteration. The execution of these operations will be further denoted as ``initialization procedure''. 4
Estimation of algorithms performance We estimate the performance of the considered algorithms on the base of the following assumptions: the matrices A1(1i) have the same fill-in factor k11 , and the matrices A1(2i)j , A2(1i)j - the same fill-in factor k2 which does not dependent on the block size; A percentage k f of the matrices A(i)
11 for which the full inversion procedure is applied does not depend on the block size; The inverse matrices of A1(1i) are completely filled; The number of elements that each of the processes must synchronize before performing a parallel multiplication is considered proportional (with the coefficient of proportionality equal to ks ) to the number of non-zero elements in the block of matrix rows. We assume that each of the processes performs only one pair of data exchange operations doing the synchronization; The submatrices are considered evenly distributed across the system.
Under these assumptions, the time sequential algorithm spend to perform
calculations on (
2 n 2(k f (1 k f )k11 ) (N 1)k2 4n n(N 1)]
T1(n, N ) = k p[ N where k p is the performance coefficient.
The minimal time is achieved here for the number of blocks
N = Nm = 2 k f (1 k f )k11 k2 n.
The time spent on computation of the matrices with an assumption that the computational complexity of one step of the incomplete basis matrix method is proportional to the maximal number of non-zero elements in the rows of the matrices A(k) 11 can be estimated as follows: A (k) 11 n3 n2 Ti1(n, N ) = Nk p2 k f N 3 k11 N 2
where k p2 is the performance coefficient.
(
Taking into account (
TP (n, N ) = P1 T1(n, N ) 2 n ks ((1 k f )k11 ke P
N 1
N k2 )
P 1
P
k f (2N (2 k f P 2k f ) 2)ka where P is the number of processes.
The minimal time spent on data exchanges for a fixed N is achieved when the
number of processes is equal to
(
We evaluate the computation time of the parallel algorithm using the linear model of a pairwise data exchange operation:
Ts (n) = n/ke ka Pm (N ) = 2nks ((1 k f )k11 2kekak f N
N 1
N k2 ) .
When doing parallel computations, the time spent to calculate the matrices A (k) 11 can be estimated with the assumption that each of the processes performs incomplete
N
or complete inversion of P submatrices (we consider this number to be integer) as
TiP (n, N ) = P1 Ti1(n, N ) NP (P 1)( Nk121Pnk2e ka (1 k f )k11 k f N 2nP2 ke ka ).
(
Experimental analysis of the algorithm’s performance We apply the proposed parallel modification of the preconditioner to solve linear systems arising while modelling stress-strain state of the soil under a dam using the model of dynamic consolidation of saturated soils as described in [11].
The two-dimensional model has the form of the hyperbolic differential equations' system with respect to 8 unknown function - the horizontal and the vertical components uw , vw , usk , vsk of water and soil skeleton displacement vectors and their first derivatives. The used boundary and the initial conditions, solution domain, and parameter values are given in [11]. The problem was discretised by the Galerkin finite element method on a rectangular mesh with respect to the space variables and by the Crank-Nicolson scheme with respect to the time variable.
While forming linear systems, their rows were ordered to form blocks related to a single unknown function. To minimize a number of synchronizations we use the Cuthill-McKee ordering of mesh nodes. Performing the diagonal scaling of rows we obtain matrices that have the following properties: non-symmetricity; loosely coupled block structure; high fill-in factor; high condition number due to different speeds of the simulated processes. As in the considered model of soil stress-strain state we have 8 unknown functions, the matrices can be easily represented as 2x2 block matrices. Such properties of matrices make them an appropriate testing case for the AMLI preconditioner that is built upon the idea of block matrix subdivision and thus can be efficient working with matrices that have loosely coupled diagonal blocks in their structure. The usage of the incomplete basis matrix method here is beneficial because matrices contain significantly ill-conditioned blocks and other methods like incomplete LU decomposition can be unstable in this case.
The performance and the convergence of the considered algorithm were experimentally studied on the SCIT-3 cluster of the Institute of Cybernetics of National Academy of Sciences of Ukraine. The eight-core nodes based on Intel Xeon 5345 processors with 2GB RAM per core were used.
The algorithm was tested on the matrix which arises on the first time step with a length = 2hours . This matrix of this linear system, further denoted as SD, has the size equal to 107856 107856 with an average number of non-zero elements in a row equal to 23.42 . Its condition number estimated by the LSMR algorithm [17] has an order of 1012 .
Linear system was solved using the BiCGStab [1] algorithm with the proposed modification of the AMLI preconditioner. The BiCGStab was used because of nonsymmetric nature of the matrix. The submatrices A1(1k) were completely inverted when 1 < 0.05 . This resulted in the complete inversion of : 50% of the submatrices. The acceptable accuracy of the calculations (here and further - the value below which the residual has to be reduced on BiCGStab iterations) was equal to 109 .
The use of the BiCGStab algorithm itself, or with the procedure of a diagonal scaling of rows, did not allow obtaining solutions of the required accuracy. The computation time when the BiCGStab was used without the usage of any preconditioner was equal to 140 ms per iteration.
To experimentally assess the accuracy of the estimates (
The performance estimate (
A (k)
The time spent on the calculations of the matrices 11 (initialization procedure)
depending on the block size in the case of the sequential algorithm was measured for
the block size ranging from 125 to 4000. Then, the time was assessed by the estimate
(
Further, we studied an influence of the size of blocks on which the matrix is split on the iterative method’s convergence. The number of iterations required to achieve 109 order of accuracy was measured along with the overall time spent by the BiCGstab algorithm and the proposed preconditioner running on a single processor core. Here, an increase of the block size leads to a decrease of the fill-in factor of the matrices A(k)
11 . As the incomplete basis matrix method performs better in terms of accuracy on highly filled sparse matrices, the decrease of the fill-in factor results in the increase of the number of completely inverted matrices. This accelerates the convergence of the iterative algorithm (Fig. 2) along with the significant increase of its running time. The block size for which the required accuracy is reached in the minimal time was in the range from 75 to 125 elements.
A(k)
The decrease of the number of fully inverted submatrices 11 reduces the time
spent on preconditioner initialization and thus significantly decreases the running time
the parallel algorithm in the case when a small number of computational resources are
used. On the other hand, the speed-up of the initialization operation is higher for the
k
higher values of f , so the time spent on initialization becomes close to equal for the
k
different values of f when the number of processes increase. The time spent on the
parallel calculations on (
The efficiency of the proposed parallel preconditioner was compared with the efficiency of known algorithms implemented in the hypre software package v.2.10.0b [18]. The only combination of an iterative method and a preconditioner implemented in the hypre which allowed obtaining solutions of 109 accuracy order was the Gmres [1] and the AMG preconditioner [3] (we will further denote this pair as the Gmres+AMG). The total time needed to obtain the solution of the linear system SD with 109 order of accuracy is given in Fig. 3. Here the proposed parallel preconditioner was used with the block size equal to 250 and 500. It was faster than the Gmres+AMG when the number of processes was less than 20 with the block size equal to 500. However, the scalability of the proposed parallel algorithm was lower comparing with the Gmres+AMG. The parallel computation scheme was constructed for the AMLI preconditioner with incomplete inversion of matrix blocks using the so-called basis matrix method and the theoretical estimates of its performance were obtained and experimentally verified.
For the linear system obtained while modelling stress-strain state of soil which has the size equal to 107856 107856 , the condition number of 1012 order, the maximal achieved speed-up was 24% when scaling the algorithm on 8 processes for the block size equal to 500 elements. The low scalability of the algorithm here can be explained by the large number of data exchange operations that must be performed before each matrix-vector multiplication.
We compared the characteristics of the proposed preconditioner with the AMG used as preconditioner in the Gmres algorithm implemented in the hypre software package. The scalability of the proposed algorithm was generally worse, but it had better performance for the considered significantly ill-conditioned linear system when solving it on less than 20 processor cores.