Parallel Substructuring Method With Memory Cost Limits?

Abstract

In decomposition methods, the memory costs for solving the interface problem are increasing signi cantly with the number of subdomains increasing. Using of the substructuring method allows you to reduce the number of iterations when the system is being solved. At the same time, the global boundary sti ness matrix takes up more memory in comparison with global sti ness matrix. This imposes restrictions on the maximum size of the problem for which you can apply this method. Di erent approaches to reduce the costs and limitations of memory on stage of the construction and solving of the interface system exist. A layer-by-layer approach to the partitioning of a triangulated multiply connected domain into connected subdomains without branching of inner boundaries was presented. This makes it possible to avoid con icts with concurrent operations of the assembly type without using synchronization and critical sections. Parallel algorithm of the construction global boundary sti ness matrix with distributed storage of the matrix are considered when implementing using OpenMP and MPI technologies. This approach allows not only to reduce the limits on the maximum size of the solved problem, but also to resolve con icts of shared memory access by increasing the number of independent parallel tasks.

Keywords

• Parallel algorithms
• Decomposition method
• Substructur- ing method
• OpenMP
• MPI

-

The substructuring method is the non-overlapping domain decomposition [1]. However, the global boundary sti ness matrix has more nonzero elements with smaller dimension in comparison with global sti ness matrix [2]. Due to these properties, even the use of sparse matrix storage formats requires considerable memory costs. Block [3] and hierarchical [4] methods with low-rank approximation have received great development due to the current trends in reducing the amount of memory per processor core. In this paper, parallel substructuring method with distributed storage of the matrix is considered. The distributed ? The reported study was partially funded by RFBR according to the research project 17-01-00402-a. construction makes it possible to apply this approach not only in the classical substructuring method, but also in the hierarchical or block one, or in the case of constructing a preconditioner on the based on the substructuring method.

The second memory limit is the access speed. To accelerate calculations, it is necessary the number of simultaneous accesses to the same memory area both during reading and writing. To increase the parallelism of the method, it is required to arrange the grid and elements in such a way that independent threads simultaneously work with their own sections of memory without overlapping. In this case, We propose using a layer-by-layer partitioning method. In the course of this method, any unstructured mesh can be arranged.

The paper presents the results of parallel version of the substructuring method realized with the help of OpenMP and MPI technology. The potential for parallelization with the help of CUDA technology for computation on graphic accelerators is considered 2

Layer-by-layer Partitioning

We generalize the ordering of the 3D unstructured mesh T on a subdomains without branching of internal boundaries using layer-by-layer partitioning. We de ne in T the some set of the cells T T , assuming L1 = T . We construct the layers Lj = f 2 T j Lj T Lc 6= ? when j = c 1g, j 2 [2; m 1], according to the Algorithm 1. Here m { number of the layers, - cell of the mesh T .

Algorithm 1: Partitioning unstructured mesh T on layers.

Input data: T | unstructured mesh of cells

f g T .

Result: fLj g | set of the layers of cells. 1 De ne the cells corresponding to the vertices of T ; 2 8 2 T nd cells having a common vertex;

/* Form the layers of mesh cells 3 Let j = 1 and layer Lj = T ;

while Pjc=1 jLcj < jT j do 4 Construct layer Lj+1 of cells T ; 5 Let j j + 1; 2 T , T | given set */

The result of the Algorithm 1 is an ordered set of m layers Lj , in which the layers adjacent to Lj have numbers j 1 and j + 1. The resulting layers are discrete analogs of a domain of R3 of type (g; 2) : g > 0. All the vertices of the cells lie on the surface bounding the layer. In general, layer Lj consist of a set of sublayers fL(c)g, where each sublayer L(c) is the set of cells in Lj having at least j j one intersection with the cell in Lj . If some layer Lj consists of one sublayer L(jc), we assume that the layer Lj is connected.

Thus, the problem of forming subdomains of an unstructured mesh fTi : Si Ti = T g based on obtained partitioning on layers Lj is as follows:

nd the union of the layers Lj in the set of subdomains fTig satisfying certain restrictions. One of the conditions is the connectivity of the resulting subdomains Ti. The number of layers included in the subdomain is bounded below by two layers, which is a su cient condition for the existence of vertices in Ti that do not belong to the boundary Ti.

In order to construct connected union of layers in the subdomain of the mesh, it is necessary to nd sublayers and their connections. The algorithm 2 provides a search for sublayers and the de nition of their connections using a dual graph Gd(V; E) of connected cells in Lj along the faces. Further, the connected components Gd(V; E) are found. The set of the mesh cells corresponding the connected components Gd(V; E) will be denoted as sublayers L(jc); c = 1; 2; : : : ; j , where j = jfL(jc)gj | the number of the sublayers in layer Lj , = max j .

j Algorithm 2: The union of the layers fLig in linked subdomains fTlg.

Input data: fLj g | set of the layers of cells.

Result: fTig | set of the mesh T subdomains.

/* Find sublayers in each layer Lj */ 1 for Lj ; 8j 2 [1; m] do 2 Construct dual graph Gd(V; E) of the links of the cells of Lj along the faces; 3 4

De ne connected components Gd(V; E);

Assign sublayers L(jc) to the connected components Gd(V; E); 8 9 10 11

5 Construct the graph of sublayers G(V; E);

/* Form subgraphs Gl of th graph G(V; E): 6 i = 1, m = 1, mark(Lj ) = 0; 8 j 2 [1; m]; 7 while m < m do /* Find an unmarked layer with j = while > 1 do

Adjoin adjacent layers to the initial layer of the domain; De ne | the number of the connected components Gi; j j + 1; 12 i

i + 1; if v 2 VG ^ deg(v) = 1 ^ v , L(c) then

j

Adjoin fvg , Lj to Gj , Lj 1, j 2 [1; i].

Partitioning/ordered mesh T , assuming Ti , Gi. */ */

The Algorithm 2 extends the partitioning on layers [5] on multi-connected domains. Partitioning without branching on the basis of the neighborhood relation reduces to the use of the Algorithm 1 for the layer-by-layer partitioning of the triangulation T and the union layers Lj ; 8j 2 [1; m] by the Algorithm 2 in th subdomains Ti.

Substructuring Method

k( 2 ), and k6np = S

k=1 ules np such that Let the mesh domain be split into n non-overlapping subdomains = i6Sn i( 1 ), where i( 1 ) T j( 1 ) = ;, when i 6= j. Formally unite the subi=1 domains, preserving the structure of the decomposition (interior nodes remain interior, boundary { boundary). We introduce the decomposition of the second level with number of the subdomains equal number of the computational modi6(k+1) nnp

S i=k nnp

For decomposition the rst level, the system of equations is formed in such a way that the unknowns corresponding to the interior and boundary nodes take the form: ( 2 ) = k i( 1 ). 0 1(A(I1I))

0

A(B1B) 1(A(I1B)) 1 0 1u( 1 ) 1 0 1f ( 1 ) 1

I I 2(A(I1B)) CCCA BB@ 2u(I1) CC = BB 2fI( 1 ) CC :

A @ A uB fB

In this system, the unknowns stand for the boundary nodes are found from the solution of the system S( 1 )uB = f~B. Here f~B = i6n X (fB i=1

i(A(B1I))i(A(I1I))( 1 )ifI( 1 )) is the right-hand-side vector and S( 1 ) is the global boundary sti ness matrix (also known as the Schur complement matrix [2]) and is the sum of the local boundary sti ness matrices Si( 1 ) = A(B1B) i(A(B1I))i(A(I1I))( 1 )i(A(I1B)) such that i6n S( 1 ) = P Si(i), where index i corresponds to subdomains of the rst level.

i=1 4

Construct Global Boundary Sti ness Matrix

Global boundary sti ness matrix S( 2 ) corresponding to these subdomains are k np supplemented by zeros up to size global matrix S( 1 ) in this S( 1 ) = P Sk( 2 ). Note k=1 that, Matrices S( 2 ) do not depend on each other during the formation stage. This k excludes con icts when reading elements data and writing results of summing of local matrices Si( 1 ).

For reduce computational costs of RAM, each computational modules forms one row l of the matrix S( 1 ) at the same time. The order of the construction i of rows is given by subdomains. Parallel algorithm of the construction for np computational modules using OpenMP technology on the example of two subdomains would be as follows (note that np 6 n ):

Algorithm 3: Parallel algorithm of the global boundary sti ness matrix construction S( 1 ) = Pnp

k=1 Sk( 2 ) for OpenMP 1 2 3 4 5 // Forming of S( 2 )

1 for i = 1 : : : n =np do

Function of the inverse matrix i(A(I1I)) 1 ; // Parallelize loop by

OpenMP for l = 1 : : : nBi do v = iA(B1I)(l; ) ; v = v i(A(I1I)) 1 ; S( 1 )(l; ) = iA(B1B)(l; ) v iA(I1B) ;

i Record the S( 1 )(l; ) elements

i in the corresponding row of the matrix S( 2 ) ;

1 // Forming of S( 2 )

2 for i = n =np : : : n do

Function of the inverse matrix i(A(I1I)) 1 ; // Parallelize loop by

OpenMP for l = 1 : : : nBi do v = iA(B1I)(l; ) ; v = v i(A(I1I)) 1 ; S( 1 )(l; ) = iA(B1B)(l; ) v iA(I1B) ;

i Record the S( 1 )(l; ) elements

i in the corresponding row of the matrix S( 2 ) ;

2

Based on the computing experiment (See Table 1), run-time costs on each subdomain have the following order (top-down): 1. Inverse of the matrix i(A(I1I)) (Step 1 of Algorithm 3); 2. Record the row l of the local matrix S( 1 )(l; ) in the corresponding row of i the matrix S12 (Step 5 of Algorithm 3); 3. Operations of taking the row from the matrix, the matrix-vectors multiplications and the di erence of the vectors (steps 2-4 of Algorithm 3, denote by MatVec).

Taking the row from the matrix, the matrix-vector multiplications and the di erence of the vectors (see steps 2-4 of Algorithm 3) are parallelized using OpenMP technology, that for these operations gave an acceleration close to linear. Also, in the future, it is planned to parallelize these operations for use on graphics accelerators using CUDA technology.

Recording of the row l of the local matrix Si( 1 )(l; ) to the corresponding row of the matrix S( 2 ) (see step 5 of Algorithm 3) largely depends on the matrix k storage formats used. Matrices S( 2 ) is formed in a compressed CSR format, k because of its large size. Using arrays of vectors (C++ container std::vector) allows to reduce memory costs to a minimum, but the complexity of adding new elements to such storage formats is O(nB), where nB is the number of equations corresponding to the boundary nodes To reduce the costs of adding new elements of the row, we use sorted associative containers (C++ container std::map) and not sorted (C++ container std::unordered_map). In the sorted container, the complexity of nding and adding a new element is O(ln nB) (due to its red-black tree structure), but it signi cantly increase memory costs to store to store additional two pointers. The non-sorted container, for these operations has complexity form O( 1 ) to O(nB) (depending on the hash function) and request smaller memory sizes, compared to the container std::map, which is optimal for use in the constructing of matrices Sk( 2 ).

For a given non-singular matrix A there exists a decomposition on upper triangular matrix U and a lower triangular matrix L. If matrices A, U , L are invertible, then A 1 = U 1L 1. In the construction L and U , products with nonzero elements of the parent matrix were summed, since the factorization war carried out for matrices stored in the CSR format. Only sequential implementations of the LU-factorization on CPU was considered. This approach showed the best results for small matrices. Features of this algorithm do not allow to e ective use all the capabilities of the massively parallel GPU architecture.

A parallel version of the Sherman-Morrison method implemented using OpenMP technology was considered.

Algorithm 4: Parallel algorithm of Sherman-Morrison method for (A(I1I)) 1 1 x; y - vectors of dimension n ; for k = 1 : : : n do

// OpenMP 2 3 xT = A(I1I)(k; ) ek; // OpenMP y = xT Ainv; // Loop of OpenMP on i for i = 1 : : : n do for j = 1 : : : n do

Ainv(j; i) = Ainv(j; i) 4 (A(I1I)) 1 = Ainv ; (yi

Ainv(j; k)=(yk + 1)) ;

LU-factorization is faster than the Sherman-Morrison method in sequential implementations. In spite of it, LU-factorization does not parallelize well. Using parallel version of the Sherman-Morrison method provided a parallel acceleration about two times. To reduce the number of operations and the computational load, matrix i(A(I1I)) 1 is stored in its entirety before the completion of the constructing of S( 1 ) in the current subdomains, since it is not known how many i non-zero elements contain an inverse matrix. To solve the system of equations on the inner subdomains nodes iA(I1I)iuI = ifI iA(I1B)uB Krylov subspace methods [7] are used.

In the future, parallelization on the GPU is planned. To increase the parallelism, the order of the formation of the local boundary sti ness matrices will be given by layer-by-layer partitioning. Thus, a block of parallel GPU threads will process one row of the matrix S( 1 ). In total it will be possible to generate i i blocks that are able to write in di erent places of the same array without overlap and con icts due to layer-by-layer ordering. 5

~

Solving the S( 1 )uB = fB System Implementation of conjugate gradients method for solving the system S( 1 )uB = f~B on one of np computing node is shown in Algorithm 5. Here the np system matrix is represented in the form S( 1 ) = P Sk( 2 ). Index k corresponding k=1 to the number of computational module is omitted. Algorithm 5: Algorithm of the preconditioned conjugate gradient method for distributed matrix 1 u; r; p; q; z 2 RnB ; 2 i = 0; 10 11 12 13 14 15 16 17 3 r0 = f~B ; 4 u0 0 ; 5 z0 = M r0 ; // M (S( 1 )) 1 -- preconditioner 6 p0 = z0 ; 7 0 = (r0; z0) ; // sync point

while jjrijj2=jjf jj2 > " do 8 Assembly pi ; // Using MPI::Allreduce 9 qi = S( 2 )pi ; // k independent operations on each MPI process i = (ri; zi)=(qi; pi) ; // MPI::Allreduce lead to sync point ui+1 = ui + ipi; ri+1 = ri iqi; zi+1 = M ri+1 ; // M (S( 1 )) 1 -- preconditioner i+1 = (ri+1; zi+1) ; // MPI::Allreduce lead to sync point i+1 = i+1= i; pi+1 = ri+1 + i+1pi; i = i + 1 ;

In a matrix-vector product (See steps 8-9 of Algorithm 5) operation qi = S( 1 )pi is performed in two steps. First, the assembly of the vector pi from all MPI processes, then np independent operations (qk)i = S( 2 )pi. The result of knp products is the set of vectors qk such that qk 2 RnB and q = P qk. Summed of k=1 vector q are not required. Vector is stored by parts on various compute nodes.

Parallelization of the operations with vectors is realized as follows: operands are parts of vectors consisting of components corresponding to mesh nodes of the subdomain k( 2 ) for each parallel MPI process. The result vector assembly are not required, each MPI process jobs with its local vectors, and synchronization occurs when the scalar products (r; z), (q; p), (r; r) are computed. Scalar products are performed in two steps. First, local scalar products for the part vectors in parallel MPI processes, then summation of the local sums (on this step, implicit synchronization occurs). Inside the MPI process, vectors operations are parallelized using OpenMP technology with number of the thread equal number of the CPU cores on computational nodes.

To parallelize the GPU, it is planned to use an approach similar to that presented in paper [2].

Results

Numerical experiments were carried out on the test problem of the theory of elasticity. The computational domain had the form of a parallelepiped. Mesh consisted of 539000 hexagonal cells. 8 threads of the OpenMP are running inside each MPI process. In Table 2, the result of experiment are presented for dimension on 512 and 1024 subdomains with using 2, 4, 6 MPI processes. The presented results are obtained on 6 computational nodes, each of which contains two processors Intel Xeon E5-2609 and 64GB RAM. For n = 512, this problem took approximately 100GB of RAM, for n = 1024, problem took approximately 80GB of RAM.

The considered algorithm allows to signi cantly reduce the memory limits for the substructuring method. The partitioning of the matrix S( 1 ) and its distributed storage makes it possible to realize coarse-grained parallelism at the formation stage of global boundary sti ness matrix, to reduce exchanges and to exclude synchronous access to one memory cell during the local matrices addition. At the same time, an increase in the number of involved computational modules permits achieving an acceleration close to linear at the stage of matrix construction. At the solver stage, the distributed storage of the matrix S( 1 ) allows solving systems with a dimension proportional to the memory size of the involved computational modules. Scalability is limited to assembling the results of scalar products of vectors and depends only on the restrictions on the speed to data exchange over the network.

References

1. Toselli A. , Widlund O.B. : Domain Decomposition Methods | Algorithms and Theory. Springer Series in Computational Mathematics, Vol. 34 ( 2005 )

2. Kopysov S.P. , Kuzmin I.M. , Nedozhogin N.S. , Novikov A.K. , Sagdeeva Y. A. : Hybrid Multi-GPU solver based on Schur complement method . Lecture Notes in Computer Science , vol. 7979 , pp. 65 - 79 ( 2013 )

3. Amestoy P.R. , Buttari A. , L'Excelent J-Y. , Mary T. : On the Complexity of the Block Low-Rank Multifrontal Factorization . SIAM Journal on Scienti c Computing , vol. 39 ( 4 ), pp. A1710 - A1740 ( 2017 )

4. Chen C. , Pouransari H. , Rajamanickam S. , Boman E.G. , Darve E.: A distributedmemory hierarchical solver for general sparse linear systems . Parallel Computing , vol. 74 , pp. 49 - 65 ( 2018 )

5. Novikov A.K. , Piminova N.K. , Kopysov S.P. , Sagdeeva Y.A. : Layer-by-layer ordering in parallel nite composition on shared-memory multiprocessors . IOP Conf. Ser.: Mater. Sci. Eng ., vol. 158 ( 1 ), p. 158 ( 2016 )

6. Sherman J. , Morrison W.J.: Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix . Ann. Math. Statist., vol. 21 ( 1 ), pp. 124 - 127 ( 1950 )

7. Saad Y. : Iterative Methods for Sparse Linear System , Boston: PWS Publiching company, p. 447 ( 1996 )