# Parallel Multiprojection Preconditioned Methods Based on Subspace Compression

Authors: Byron E. Moutafis; Christos K. Filelis-Papadopoulos; George A. Gravvanis
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

During the last decades, the continuous expansion of supercomputing infrastructures necessitates the design of scalable and robust parallel numerical methods for solving large sparse linear systems. A new approach for the additive projection parallel preconditioned iterative method based on semiaggregation and a subspace compression technique, for general sparse linear systems, is presented. The subspace compression technique utilizes a subdomain adjacency matrix and breadth first search to discover and aggregate subdomains to limit the average size of the local linear systems, resulting in reduced memory requirements. The depth of aggregation is controlled by a user defined parameter. The local coefficient matrices use the aggregates computed during the formation of the subdomain adjacency matrix in order to avoid recomputation and improve performance. Moreover, the rows and columns corresponding to the newly formed aggregates are ordered last to further reduce fill-in during the factorization of the local coefficient matrices. Furthermore, the method is based on nonoverlapping domain decomposition in conjunction with algebraic graph partitioning techniques for separating the subdomains. Finally, the applicability and implementation issues are discussed and numerical results along with comparative results are presented.


---

## Body

## 1. Introduction

Let us consider a sparse linear system of the following form:(1)Ax=b,where A is the coefficient matrix, b is the right-hand-side vector, and x is the solution vector. In order to solve the linear system, in (1), a direct or an iterative method can be used. However, a direct method is computationally expensive with excessive memory requirements. Krylov subspace iterative methods, compare [1], have gained the attention of the scientific community during the recent decades, due to their efficiency and memory requirements for solving large sparse linear systems. Preconditioned iterative methods improve the convergence rate and have been used extensively for solving large sparse linear systems. The left preconditioned form of the linear system, in (1), is as follows:(2)M-1Ax=M-1b,where M-1 is a suitable preconditioning matrix. Moreover, there is the right and the symmetric preconditioned form, compare [1]. The preconditioner M-1 has to satisfy the following conditions: (i) M-1A should have a “clustered” spectrum, (ii) M-1 can be efficiently computed in parallel, and (iii) the product of “M-1 by vector” should be fast to compute in parallel, compare [1, 2]. The domain decomposition method has been extensively used for solving (very) large sparse linear systems in parallel systems, compare [3–8]. Furthermore, domain decomposition-type methods have been combined with Krylov subspace methods or have been used in hybrid schemes in conjunction with direct solvers, compare [6, 9–11]. During the last decade, research efforts were directed towards solvers for general sparse linear systems that have improved convergence behavior and are scalable.Various overlapping domain decomposition methods based on Lagrange multipliers have been proposed. The Finite Element Tearing and Interconnect (FETI), compare [12], and its variants, compare [13, 14], are scalable and efficient solvers for second-order and fourth-order partial differential equation problems. Related methods are the Balancing Domain Decomposition (BDD), compare [15], and the Balancing Domain Decomposition by constraints (BDDC), compare [16], which are multiplicative and additive domain decomposition methods, respectively. The aforementioned methods are suitable for solving linear systems derived from the finite element method. Additionally, there are domain decomposition-type solvers using graph partitioning algorithms. These algebraic domain decomposition methods do not require any geometric property of the problem as the partitioning is based on the corresponding graph of the coefficient matrix. Parallel hybrid algebraic domain decomposition-type solvers that combine direct and iterative methods have been proposed, compare [17, 18]. Additionally, algebraic domain decomposition method in conjunction with algebraic multigrid method, compare [19–22], has satisfactory convergence behavior and is suitable for parallel systems. An algebraic domain decomposition method, based on semiaggregation and subspace compression techniques, namely, multiprojection method with subspace compression (MPMSC), for solving large sparse linear systems is proposed. Several efficient software libraries for graph partitioning, such as Metis, compare [23], or Scotch, compare [24], can be used for algebraic partitioning. The aggregation techniques have been previously used in the context of algebraic multigrid methods and have been shown to be efficient for a wide variety of model problems, compare [25–27]. In the MPMSC scheme, a direct method is used for solving the local linear systems derived from the semiaggregated domains. The semiaggregation procedure results in multiple grids composed of a fine part, corresponding to the unknowns of a subdomain, and aggregated part, corresponding to unknowns related to other subdomains. The proposed scheme utilizes approximations for the fine part of each subdomain, obtained with respect to its aggregated components, as opposed to aggregation based algebraic multigrid methods, which obtain approximations for the solution of the entire domain, through smoothing, based on a hierarchy of coarser grids. Moreover, the proposed scheme does not require the computation of a hierarchy of coarser grids. The subspace compression technique is an algorithmic procedure, based on breadth first search with limited depth on the adjacency matrix of the subdomains. By this technique the neighboring subdomains are discovered and aggregated together to an aggregated (coarse) component. The aggregation process reduces the number of aggregated (coarse) unknowns of the local linear systems, resulting in a substantial reduction in memory requirements. The computation of the rows and columns corresponding to the aggregates, during the formation of the local coefficient matrices as well as the subdomain adjacency matrix, is reused to avoid recomputation. Avoiding the recomputation enhances the performance of the proposed scheme. Moreover, the rows and columns corresponding to the newly formed aggregates have a large number of nonzero elements and thus are ordered last to further reduce fill-in during the factorization of the local coefficient matrices. Further, the number of aggregates of the multiprojection method with subspace compression is controlled by a parameter that describes the level sets of aggregates to be compressed in a single aggregate.In Section2, the multiprojection method with subspace compression is presented, implementation details are given, and the aggregation strategy is discussed. In Section 3, numerical results illustrating the applicability and efficiency of the MPMSC scheme are presented. Moreover, comparative results for the performance and convergence behavior are given.


## 2. Multiprojection Method with Subspace Compression

The multiprojection method with subspace compression is a domain decomposition-type method for solving large sparse linear systems (see (1)). The graph, corresponding to the coefficient matrix A, has vertices, the unknowns of the linear system, in (1), and edges, the nonzero elements of the coefficient matrix A. Since A is expressed as a graph, the algebraic properties of the problem can be utilized in order to partition the graph into ndoms nonoverlapping subgraphs (subdomains). The graph partitioning process is performed using the algorithm provided in the Metis software library, compare [23]. Consider the domain Ω partitioned into ndoms nonoverlapping subdomains Ωj with j=0,…,ndoms-1, so that(3)Ω=⋃j=0ndoms-1Ωj.The number of vertices, in each subdomain Ωj, is denoted by mj and Sj is the index set defined by the following expression:(4)Sj=i1,i2,…,imj,where ik are the indices corresponding to the mj vertices of each subdomain.


### 2.1. Subspace Compression

Let us consider the subspaceK, of the domain Ω, which contains ndoms aggregated (coarse) components, corresponding to the subdomains Ωj, and the prolongation matrix R from K onto Ω, R:Κ→Ω, which is denoted according to the following expression:(5)Ri,j=1mj,ifi∈Sj,i=0,…,n-1,0,ifi∉Sj,j=0,…,ndoms-1.The projection of A to subspace K is the matrix H of dimension (ndoms×ndoms) and is given by the following equation:(6)H=RTAR.The graph G(H)=(v,e), corresponding to the matrix H, concerns the relations between the aggregated (coarse) components. The subspace compression technique is based on the reaggregation of the ndoms (number of subdomains) aggregated (coarse) components into naggs (number of aggregates). The number of aggregates is not known a priori and thus is a result of the subspace compression algorithm. Breadth first search with limited depth (BFS-LD) is performed on G(H) to find the distance d neighborhood. The components in a distance d neighborhood of G(H) are reaggregated together. Then the neighborhood’s vertices are removed from the set v and the BFS-LD is performed again on the new graph. When every vertex is removed from the set v, the procedure is terminated. The order in which the reaggregation process is performed may result in a slight variation of the number of aggregates; however this does not affect the efficiency or the convergence behavior of the proposed scheme significantly. The order of the vertices used in the reaggregation scheme is lexicographical with respect to the enumeration of the subdomains by the graph partitioning scheme. Other approaches can be used for the reaggregation process, compare [26, 27]; however the computational work of these schemes will affect negatively the performance of the preprocessing phase, especially for large numbers of subdomains. The subspace compression algorithm is as follows:(1)
U←K, Agg←Agg0,Agg1,…,Aggnaggs-1=Ø,Ø,…,Ø.(2)
For eachyi∈U oruntilU⊃Ø
Aggk←Nyid∩U, where Nyid is the d-distance neighborhood of component yi and k is the index of the aggregated (coarse) components(7)U←U∖AggkEnd ForThe parameter d of the neighborhood Nyid specifies the maximum distance between the vertex yi and any other vertex in the neighborhood. The value of d is selected with respect to the adjacency matrix G(H), considering that the larger the value of d, the less the number of aggregates. Distance one neighborhood is chosen by default; however for parallel systems with limited memory resources, distance two or distance three neighborhood should be chosen to reduce the memory requirements of the proposed preconditioning scheme. An increase in the value of d might affect negatively the convergence behavior of the proposed scheme. The computational complexity for finding the aggregates depends on the number of subdomains and it is a searching procedure on the graph G(H) with negligible execution time compared to the whole initialization phase of the preconditioning scheme.


### 2.2. Semiaggregation

Let us consider the semiaggregated subdomainsZj, containing the mj fine components of subdomain Ωj and naggs aggregated components. Furthermore, the semiaggregated Zj subdomains are associated with their respective prolongation matrices Vj∈Rn×mj+naggs:Zj→Ω, with j=0,1,…,ndoms-1. The matrices VjT:Ω→Zj are restriction operators mapping an arbitrary full vector of the domain Ω into a vector of the semiaggregated subdomain Zj with size (mj+naggs). The first mj columns of the matrix Vj are the columns of the identity matrix with dimension (n×n), corresponding the indices Sj. The remaining (naggs) column vectors of Vj correspond to the aggregates. The components of the column vector that are included in the corresponding aggregated subdomain have value 1/Mk and the rest of the components are zero, where Mk=∑miSi⊆Aggk, with i≠j and k=0,1,…,naggs-1.Hence, the matrixVj is defined as follows:(8)Vj=ek,k∈Sj∣qk,0≤k<naggs,where qk is the column vectors corresponding to the aggregates as they are described above. Moreover, additional prolongation matrices are required, namely, Wj. The prolongation matrices Wj of dimension (n×mj+naggs) are used to map a vector defined on the subdomain Zj to a vector defined on the domain Ω by prolonging only the components corresponding to the subdomain Ωj and discarding the aggregated (coarse) components. The ith column of the matrix Wj is given as follows:(9)Wj:,i=ei,ifi∈Sj,0,otherwise.The linear systems corresponding to each subdomain Zj are as follows:(10)Ajxj=VjTb=bj,where Aj=VjTAVj is a coefficient matrix of order (mj+naggs).The proposed method is used as a preconditioner to a Krylov subspace iterative method, such as preconditioned GMRES(m), compare [28, 29]. The product of the preconditioning matrix by a vector r can be described by the following compact algorithmic scheme:Forj=0,…,ndoms-1(11)Computeyj=WjWjTVjAj-1VjTrEnd For(12)y=∑j=0ndoms-1yj.Considering that X∗ is the exact solution and that r=b-Ax=Ax∗-x=Ae, the error ek+1 at (k+1)th iteration can be written as follows:(13)ek+1=ek-∑j=0ndoms-1WjWjTVjAj-1VjTAek,or equivalently(14)ek+1=I-∑j=0ndoms-1Qjek,where(15)Qj=WjWjTVjAj-1VjTA.WjWjT is an orthogonal projector with respect to the Euclidean inner product. Considering that if A is symmetric positive definite, VjAj-1VjTA is an orthogonal projector with respect to the A-norm inner product, compare [8, 30]. Qj is a product of the two orthogonal projections, which are orthogonal to different inner products; thus the projection is not orthogonal, compare [31]. Taking into consideration the above, the method is categorized into the class of oblique projection procedures, compare [1].The preconditioned matrix using the MPMSC scheme is given as follows:(16)MMPMSCA=∑j=0ndoms-1Qj,which is a sum of projections.Let us consider the unit square domainΩ, which is discretized with mesh size h=1/7 as depicted in Figure 1. The resulting grid is partitioned into 16 subdomains, containing 4 vertices each. The corresponding matrix A is projected onto the subdomains Zj, in (10), and the grid corresponding to the projected subdomain Z0 is presented schematically in Figure 2.Figure 1
A square domainΩ discretized with mesh size h=1/7 partitioned into 16 subdomains Ωj.Figure 2
The grid corresponding to the projection of the domainΩ onto subdomain Z0.In the MPMSC scheme,ndoms ((mj+naggs)×(mj+naggs)) linear systems of the form in (10) have to be solved in every iteration. Equivalently, the linear system, in (10), can be written in a reordered block form as follows:(17)BFFBFCBCFBCCxFxC=bFbC,where the size of BFF is (mj×mj) and the size of BCC is naggs×naggs. The vector xF is prolonged to the full solution vector x of domain Ω, and the aggregated (coarse) components are used as auxiliary. In the MPMSC preprocessing step each matrix Aj is factorized by sparse LU factorization method, compare [32]:(18)Aj=LjUj.Therefore, the preconditioning operation can be described by the following algorithmic scheme:Forj=0,…,ndoms-1(19)Solve  theLjwj=VjTx(20)Solve  theUjzj=wj(21)Computeyj=WjWjTVjzjEnd For(22)y=∑j=0ndoms-1yj.The solution, in (19)-(20), of the local linear systems is performed inherently in parallel; thus each core computes one backward-forward substitution process. The aggregated (coarse) components are ordered last in the coefficient matrix Aj, because they have many connections with the fine components and consequently many nonzeros in their rows/columns. By this ordering, the fill-in of the L and U matrices is reduced, compare [32]. Direct solvers have memory and computational complexity limitations for large linear systems. In order to keep the size of Aj below these limits, the coefficient matrix A should be partitioned into a suitable number of subdomains. Maintaining the size of Aj almost constant, regardless of the size of A, results in almost constant LU factorization computational work.It should be stated that the MPMSC scheme is inherently parallel and thus suitable for distributed systems. The MPMSC scheme is used as a preconditioner to the Parallel Preconditioned restarted Generalized Minimum RESidual, namely, (PPGMRES(m)), method which is a parallel variant of the preconditioned GMRES(m) method proposed by Douglas et al., compare [28], and is based on Gram–Schmidt orthogonalization, compare [33].


### 2.3. Memory Requirements

The memory requirements of a sparse matrix of sizen with nnz nonzero elements stored in CSR (Compressed Sparse Row) format are denoted by(23)csrn,nnz=n+1∗sizeofint+nnz∗sizeofint+nnz∗sizeofdouble.The size of the coefficient matrix A is denoted by “n”, E is the expected value operator, and nnz() is an operator returning the number of nonzero elements of a matrix. Also, let us denote by “nnode” the average number of components that correspond to a node, which usually is equal to cores  per  node∗Emj. The restart parameter of PPGMRES(m) is denoted by “m.” MPMSC requires the following matrices and vectors in the preprocessing phase.The coefficient matrixA in CSR with memory requirements csrn,nnzA, the right-hand-side vector bn∗sizeofdouble, the partitioning vector (n∗sizeof(int)) and (number of subdomains/nodes) vectors that map the local index of the components of each subdomain to the global index ((number  of  subdomains/nodes)∗Emj∗sizeof(int)), and also the (number of subdomains/nodes) sparse (CSR) matrices Vj with memory requirements csrn,n are required. The matrices Wj are not stored, since they can be applied explicitly. Moreover, the LU factorization process requires a large fraction of the total memory requirements. Each node requires (number of subdomains/nodes) Lj and Uj sparse (CSR) matrices with memory requirements csrmj+naggs+1,EnnzLj and csrmj+naggs+1,EnnzUj. The MPMSC requires also the sparse (CSR) matrices Hcsrndoms,nnzH and Rcsrn,n for the subspace compression process; however they are temporary as they are deleted before the LU factorization step. The PPGMRES(m) algorithm has the following extra memory requirements: six auxiliary vectors (nnode∗sizeof(double)), the dense matrix retaining the basis vectors of the Krylov subspace n∗m+1∗sizeofdouble, the upper Hessenberg matrix Tm+1∗m∗sizeofdouble, and three auxiliary vectors for the Givens rotation m+1∗sizeofdouble.


## 2.1. Subspace Compression

Let us consider the subspaceK, of the domain Ω, which contains ndoms aggregated (coarse) components, corresponding to the subdomains Ωj, and the prolongation matrix R from K onto Ω, R:Κ→Ω, which is denoted according to the following expression:(5)Ri,j=1mj,ifi∈Sj,i=0,…,n-1,0,ifi∉Sj,j=0,…,ndoms-1.The projection of A to subspace K is the matrix H of dimension (ndoms×ndoms) and is given by the following equation:(6)H=RTAR.The graph G(H)=(v,e), corresponding to the matrix H, concerns the relations between the aggregated (coarse) components. The subspace compression technique is based on the reaggregation of the ndoms (number of subdomains) aggregated (coarse) components into naggs (number of aggregates). The number of aggregates is not known a priori and thus is a result of the subspace compression algorithm. Breadth first search with limited depth (BFS-LD) is performed on G(H) to find the distance d neighborhood. The components in a distance d neighborhood of G(H) are reaggregated together. Then the neighborhood’s vertices are removed from the set v and the BFS-LD is performed again on the new graph. When every vertex is removed from the set v, the procedure is terminated. The order in which the reaggregation process is performed may result in a slight variation of the number of aggregates; however this does not affect the efficiency or the convergence behavior of the proposed scheme significantly. The order of the vertices used in the reaggregation scheme is lexicographical with respect to the enumeration of the subdomains by the graph partitioning scheme. Other approaches can be used for the reaggregation process, compare [26, 27]; however the computational work of these schemes will affect negatively the performance of the preprocessing phase, especially for large numbers of subdomains. The subspace compression algorithm is as follows:(1)
U←K, Agg←Agg0,Agg1,…,Aggnaggs-1=Ø,Ø,…,Ø.(2)
For eachyi∈U oruntilU⊃Ø
Aggk←Nyid∩U, where Nyid is the d-distance neighborhood of component yi and k is the index of the aggregated (coarse) components(7)U←U∖AggkEnd ForThe parameter d of the neighborhood Nyid specifies the maximum distance between the vertex yi and any other vertex in the neighborhood. The value of d is selected with respect to the adjacency matrix G(H), considering that the larger the value of d, the less the number of aggregates. Distance one neighborhood is chosen by default; however for parallel systems with limited memory resources, distance two or distance three neighborhood should be chosen to reduce the memory requirements of the proposed preconditioning scheme. An increase in the value of d might affect negatively the convergence behavior of the proposed scheme. The computational complexity for finding the aggregates depends on the number of subdomains and it is a searching procedure on the graph G(H) with negligible execution time compared to the whole initialization phase of the preconditioning scheme.


## 2.2. Semiaggregation

Let us consider the semiaggregated subdomainsZj, containing the mj fine components of subdomain Ωj and naggs aggregated components. Furthermore, the semiaggregated Zj subdomains are associated with their respective prolongation matrices Vj∈Rn×mj+naggs:Zj→Ω, with j=0,1,…,ndoms-1. The matrices VjT:Ω→Zj are restriction operators mapping an arbitrary full vector of the domain Ω into a vector of the semiaggregated subdomain Zj with size (mj+naggs). The first mj columns of the matrix Vj are the columns of the identity matrix with dimension (n×n), corresponding the indices Sj. The remaining (naggs) column vectors of Vj correspond to the aggregates. The components of the column vector that are included in the corresponding aggregated subdomain have value 1/Mk and the rest of the components are zero, where Mk=∑miSi⊆Aggk, with i≠j and k=0,1,…,naggs-1.Hence, the matrixVj is defined as follows:(8)Vj=ek,k∈Sj∣qk,0≤k<naggs,where qk is the column vectors corresponding to the aggregates as they are described above. Moreover, additional prolongation matrices are required, namely, Wj. The prolongation matrices Wj of dimension (n×mj+naggs) are used to map a vector defined on the subdomain Zj to a vector defined on the domain Ω by prolonging only the components corresponding to the subdomain Ωj and discarding the aggregated (coarse) components. The ith column of the matrix Wj is given as follows:(9)Wj:,i=ei,ifi∈Sj,0,otherwise.The linear systems corresponding to each subdomain Zj are as follows:(10)Ajxj=VjTb=bj,where Aj=VjTAVj is a coefficient matrix of order (mj+naggs).The proposed method is used as a preconditioner to a Krylov subspace iterative method, such as preconditioned GMRES(m), compare [28, 29]. The product of the preconditioning matrix by a vector r can be described by the following compact algorithmic scheme:Forj=0,…,ndoms-1(11)Computeyj=WjWjTVjAj-1VjTrEnd For(12)y=∑j=0ndoms-1yj.Considering that X∗ is the exact solution and that r=b-Ax=Ax∗-x=Ae, the error ek+1 at (k+1)th iteration can be written as follows:(13)ek+1=ek-∑j=0ndoms-1WjWjTVjAj-1VjTAek,or equivalently(14)ek+1=I-∑j=0ndoms-1Qjek,where(15)Qj=WjWjTVjAj-1VjTA.WjWjT is an orthogonal projector with respect to the Euclidean inner product. Considering that if A is symmetric positive definite, VjAj-1VjTA is an orthogonal projector with respect to the A-norm inner product, compare [8, 30]. Qj is a product of the two orthogonal projections, which are orthogonal to different inner products; thus the projection is not orthogonal, compare [31]. Taking into consideration the above, the method is categorized into the class of oblique projection procedures, compare [1].The preconditioned matrix using the MPMSC scheme is given as follows:(16)MMPMSCA=∑j=0ndoms-1Qj,which is a sum of projections.Let us consider the unit square domainΩ, which is discretized with mesh size h=1/7 as depicted in Figure 1. The resulting grid is partitioned into 16 subdomains, containing 4 vertices each. The corresponding matrix A is projected onto the subdomains Zj, in (10), and the grid corresponding to the projected subdomain Z0 is presented schematically in Figure 2.Figure 1
A square domainΩ discretized with mesh size h=1/7 partitioned into 16 subdomains Ωj.Figure 2
The grid corresponding to the projection of the domainΩ onto subdomain Z0.In the MPMSC scheme,ndoms ((mj+naggs)×(mj+naggs)) linear systems of the form in (10) have to be solved in every iteration. Equivalently, the linear system, in (10), can be written in a reordered block form as follows:(17)BFFBFCBCFBCCxFxC=bFbC,where the size of BFF is (mj×mj) and the size of BCC is naggs×naggs. The vector xF is prolonged to the full solution vector x of domain Ω, and the aggregated (coarse) components are used as auxiliary. In the MPMSC preprocessing step each matrix Aj is factorized by sparse LU factorization method, compare [32]:(18)Aj=LjUj.Therefore, the preconditioning operation can be described by the following algorithmic scheme:Forj=0,…,ndoms-1(19)Solve  theLjwj=VjTx(20)Solve  theUjzj=wj(21)Computeyj=WjWjTVjzjEnd For(22)y=∑j=0ndoms-1yj.The solution, in (19)-(20), of the local linear systems is performed inherently in parallel; thus each core computes one backward-forward substitution process. The aggregated (coarse) components are ordered last in the coefficient matrix Aj, because they have many connections with the fine components and consequently many nonzeros in their rows/columns. By this ordering, the fill-in of the L and U matrices is reduced, compare [32]. Direct solvers have memory and computational complexity limitations for large linear systems. In order to keep the size of Aj below these limits, the coefficient matrix A should be partitioned into a suitable number of subdomains. Maintaining the size of Aj almost constant, regardless of the size of A, results in almost constant LU factorization computational work.It should be stated that the MPMSC scheme is inherently parallel and thus suitable for distributed systems. The MPMSC scheme is used as a preconditioner to the Parallel Preconditioned restarted Generalized Minimum RESidual, namely, (PPGMRES(m)), method which is a parallel variant of the preconditioned GMRES(m) method proposed by Douglas et al., compare [28], and is based on Gram–Schmidt orthogonalization, compare [33].


## 2.3. Memory Requirements

The memory requirements of a sparse matrix of sizen with nnz nonzero elements stored in CSR (Compressed Sparse Row) format are denoted by(23)csrn,nnz=n+1∗sizeofint+nnz∗sizeofint+nnz∗sizeofdouble.The size of the coefficient matrix A is denoted by “n”, E is the expected value operator, and nnz() is an operator returning the number of nonzero elements of a matrix. Also, let us denote by “nnode” the average number of components that correspond to a node, which usually is equal to cores  per  node∗Emj. The restart parameter of PPGMRES(m) is denoted by “m.” MPMSC requires the following matrices and vectors in the preprocessing phase.The coefficient matrixA in CSR with memory requirements csrn,nnzA, the right-hand-side vector bn∗sizeofdouble, the partitioning vector (n∗sizeof(int)) and (number of subdomains/nodes) vectors that map the local index of the components of each subdomain to the global index ((number  of  subdomains/nodes)∗Emj∗sizeof(int)), and also the (number of subdomains/nodes) sparse (CSR) matrices Vj with memory requirements csrn,n are required. The matrices Wj are not stored, since they can be applied explicitly. Moreover, the LU factorization process requires a large fraction of the total memory requirements. Each node requires (number of subdomains/nodes) Lj and Uj sparse (CSR) matrices with memory requirements csrmj+naggs+1,EnnzLj and csrmj+naggs+1,EnnzUj. The MPMSC requires also the sparse (CSR) matrices Hcsrndoms,nnzH and Rcsrn,n for the subspace compression process; however they are temporary as they are deleted before the LU factorization step. The PPGMRES(m) algorithm has the following extra memory requirements: six auxiliary vectors (nnode∗sizeof(double)), the dense matrix retaining the basis vectors of the Krylov subspace n∗m+1∗sizeofdouble, the upper Hessenberg matrix Tm+1∗m∗sizeofdouble, and three auxiliary vectors for the Givens rotation m+1∗sizeofdouble.


## 3. Numerical Results

In this section, the effectiveness and applicability of the MPMSC scheme are examined. The numerical experiments were performed on a BlueGene/P (BG/P) supercomputer with the following specifications: CPU: 1024x Quad-Core PowerPC-450 850 Mhz; RAM: 4 GB/node; interconnect: 3D-Torus network with bandwidth 5.1 GBps. The methods are designed for distributed systems with multicore nodes, using MPI and OpenMP environments, compare [34]. Each multicore node executes an MPI process, which undertakes the computations for as many subdomains as the number of its cores. The parallelization in each multicore node is achieved with OpenMP, so that every subdomain is mapped to one core. Having equal number of cores and subdomains allows the parallel computation of every LU factorization. Moreover, the software libraries BLAS (Basic Linear Algebra Subprograms), compare [35], and Metis, compare [23], have been used. The restart parameter of PPGMRES(m) was set to 20 and the maximum number of iterations was set to 500 iterations. The termination criterion was r2<10-8b2. The execution time is given in “seconds.” The number of nonzero elements in the matrix A is denoted by nnzA. The convergence behavior of the PPGMRES(m) is given by outer (inner) iterations. The inner iterations are the number of iterations after the last restart of the PPGMRES(m) and outer iterations are the number of times PPGMRES(m) was restarted. The neighborhood distance parameter d has been set to one, unless stated otherwise explicitly.The formation of the prolongation operators and theLU factorization of the local matrices are included in the initialization phase of the preconditioning scheme, which is part of the preprocessing time (Pre-Time). The time required for the graph partitioning, using Metis, compare [23], is not included in the Pre-Time. It should be noted that parallel graph partitioning schemes, such as ParMETIS, compare [36], can be used to increase parallel performance. The communication time of the PPGMRES(m) is given as Comm-Time and the GMRES-Time is the execution time of the PPGMRES(m). It should be mentioned that the Comm-Time is included in the GMRES-Time. Total-Time is the sum of the Pre-Time and the GMRES-Time. The speedup is computed using the execution on 64 cores as a baseline unless stated otherwise explicitly.


### 3.1. 3D Poisson Problem

The PPGMRES(m) in conjunction with MPMSC scheme has been used for solving the Poisson equation in three space variables subject to Dirichlet boundary conditions. The Poisson equation was discretized with the 7-point stencil. The right-hand side of the linear systems was computed as the product of coefficient matrix A and the solution vector set to 0,1,…,n-1T, where n denotes the order of the linear system. In Table 1, the convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000) are presented for various numbers of cores-subdomains. The communication time is relative to the required iterations for convergence to the prescribed tolerance and the number of cores. It can be easily seen that as the number of cores and subdomains increases, the number of required iterations for convergence to the prescribed tolerance decreases for the 3D Poisson problem. In Table 2, the convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the Block Jacobi preconditioner, for the 3D Poisson problem (n=1,000,000) with various numbers of cores-subdomains are given. It should be stated that the MPMSC scheme presents better convergence behavior than the Block Jacobi preconditioning scheme and therefore requires less overall communication time, especially for many cores where the convergence behavior of the Block Jacobi method is downgraded. Therefore, it should be noted that the MPMSC scheme has better scalability than the Block Jacobi preconditioner. In Figure 3, the speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000) is depicted. In Figure 4, the weak scalability diagram of time versus cores for MPMSC scheme is presented. The order of the linear system increases as the number of cores increases since the local linear system should be of constant order. This results in augmenting, by 10,000 more unknowns, the linear system for every added core. The Total-Time should ideally remain constant. It should be noted that the Pre-Time is almost constant, as the initialization phase of the preconditioner is an inherently parallel process. From Figures 3 and 4 it is evident that the proposed method is scalable up to large number of cores, for the efficient solution of large linear systems. It should be mentioned that partitioning the domain Ω into more subdomains improves the convergence behavior of the method, since the grid corresponding to the aggregated (coarse) components becomes finer, which improves the accuracy of the local approximations to the solution. The improved convergence behavior reduces the communications and hence improves the performance of the proposed scheme. Due to the reduced communications, while the local computational work remains almost constant, it results in superlinear speedups, especially in the case of moderate number of subdomains. For small number of cores the Pre-Time is dominant compared to the GMRES-Time and as the preprocessing phase is an inherently parallel process, both methods have satisfactory speedup. When the number of cores is 1024, the Pre-Time is a small percentage of the Total-Time; therefore the speedup reduces, since it relies on the parallelization of the PPGMRES(m). It should be noted that the GMRES-Time increases proportionally to the problem size, as global communications, related to the synchronization of distributed vectors, and sequential redundant operations of PPGMRES(m) depend on the size of global vectors, compare [28], as it can be seen in Figure 4. However, the inherently parallel preprocessing phase scales, with satisfactory speedup, up to large number cores-subdomains. Thus, the Total-Time required by the proposed scheme is improved.Table 1
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000).
Poisson 3D
Cores/subdomains
64
128
256
512
1024
MPMSC
Pre-Time
367.38
70.04
21.45
6.31
3.64
Comm-Time
23.54
13.63
10.92
9.81
8.39
GMRES-Time
122.45
73.77
56.45
50.85
41.51
Total-Time
489.84
143.81
77.90
57.16
45.14
Iterations
4 (8)
4 (13)
4 (12)
4 (7)
3 (15)
Aggregates
13
20
40
84
162
Preconditioned Relative Residual
9.24
E
-
09
9.73
E
-
09
8.90
E
-
09
9.99
E
-
09
8.91
E
-
09
Relative Residual
1.38
E
-
08
1.48
E
-
08
1.28
E
-
08
1.54
E
-
08
1.47
E
-
08Table 2
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the Block Jacobi preconditioner, for the 3D Poisson problem (n=1,000,000).
Poisson 3D
Cores/subdomains
64
128
256
512
1024
Block Jacobi
Pre-Time
368.53
69.43
21.42
6.26
3.58
Comm-Time
24.14
16.63
12.21
11.72
12.62
GMRES-Time
133.28
90.17
66.13
64.25
61.81
Total-Time
501.81
159.60
87.55
70.51
65.39
Iterations
4 (18)
5 (15)
5 (9)
5 (15)
5 (12)
Preconditioned Relative Residual
9.52
E
-
09
9.33
E
-
09
9.81
E
-
09
9.85
E
-
09
9.90
E
-
09
Relative Residual
1.21
E
-
08
1.57
E
-
08
1.60
E
-
08
1.41
E
-
08
1.42
E
-
08Figure 3
The speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000).Figure 4
The weak scalability (time versus cores) of the MPMSC scheme for 3D Poisson problem.Furthermore, the convergence behavior of MPMSC scheme for 3D Poisson problem with 10,000 unknowns per subdomain is shown in Figure5. The total number of iterations is computed as outer∗m+inner. It should be stated that the outer iterations of MPMSC scheme are almost constant for large number of subdomains. Moreover, the convergence behavior of the MPMSC scheme is related to the quality of the graph partitioning, provided by Metis, compare [23]. In Figure 6, the number of aggregates (naggs) of the MPMSC scheme against the number of subdomains for the 3D Poisson problem with 10,000 unknowns per subdomain is depicted.Figure 5
The convergence behavior of MPMSC scheme for 3D Poisson problem with various numbers of subdomains.Figure 6
The number of aggregated (coarse) components (naggs) of the MPMSC scheme for 3D Poisson problem with various numbers of subdomains.In Table3, the performance and convergence behavior of MPMSC scheme with distance two neighborhood for the 3D Poisson problem (n=1,000,000) with various numbers of cores-subdomains are presented. It should be mentioned that using distance two neighborhood in the subspace compression technique results in less aggregates compared to distance one neighborhood, which decreases the memory requirements. However, since the required number of iterations for convergence increases, the performance of the iterative scheme is downgraded.Table 3
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme with distance two neighborhood for the 3D Poisson problem (n=1,000,000).
Poisson 3D
Cores/subdomains
64
128
256
512
1024
MPMSCd=2
Pre-Time
368.53
69.43
21.42
6.26
3.58
Comm-Time
24.14
16.63
12.21
11.72
12.62
GMRES-Time
133.28
90.17
66.13
64.25
61.81
Total-Time
501.81
159.60
87.55
70.51
65.39
Iterations
4 (18)
5 (15)
5 (9)
5 (15)
5 (12)
Aggregates
6
9
19
32
51
Preconditioned Relative Residual
9.18
E
-
09
8.62
E
-
09
9.57
E
-
09
8.80
E
-
09
8.56
E
-
09
Relative Residual
1.57
E
-
08
1.18
E
-
08
1.59
E
-
08
1.38
E
-
08
1.46
E
-
08Furthermore, the parallel performance of the proposed method is compared to sequential preconditioned GMRES(m) in conjunction with dual threshold incomplete LU factorization (ILUT(p,τ)) preconditioner, compare [1]. The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=50 and τ=10-3, and various numbers of cores, for the 3D Poisson problem (n=1,000,000), is presented in Figure 7. The preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner required 2 (14) iterations for convergence to the prescribed tolerance. The PPGMRES(m) in conjunction with the MPMSC scheme performed better, in terms of total execution time, compared to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, especially for large numbers of cores-subdomains.Figure 7
The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=50 and τ=10-3, and various numbers of cores for the 3D Poisson problem (n=1,000,000).


### 3.2. UFL Matrices

In order to assess the applicability, scalability, and performance of MPMSC scheme for solving general large sparse linear systems, seven matrices from the University of Florida sparse matrix collection, compare [37], have been considered. Information concerning the aforementioned matrices is given in Table 4. The right-hand-side vectors for atmosmodd, atmosmodl, and thermal2 problems were given by the University of Florida sparse matrix collection, compare [37], whereas, for the rest of the problems, they were computed as the product of coefficient matrix A and the solution vector set to 0,1,…,n-1T.Table 4
Problems from the University of Florida sparse matrix collection.
Problem
Application
Size(A)
nnz(A)
atmosmodd
Atmospheric modeling
1270432
8814880
atmosmodl
1489752
10319760
tmt_sym
Electromagnetics problem
726713
5080961
apache2
Structural problem
715176
4817870
Cage14
DNA electrophoresis
1505785
27130349
ecology2
Ecology problem
999999
4995991
thermal2
Thermal problem
1228045
8580313The convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for the atmosmodd, atmosmodl, tmt_sym, apache2, and cage14 problems, are presented in Table 5. In Tables 6 and 7, the convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for the ecology2 and the thermal2 problems, are given. In Figures 8 and 9, the speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme for the ecology2 and the thermal2 problem, is shown.Table 5
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for various problems.
MPMSC
Cores/subdomains
64
128
256
512
1024
atmosmodd
Iterations
5 (12)
6 (12)
6 (10)
6 (17)
7 (10)
Total-Time
759.43
255.30
137.97
108.40
109.37
atmosmodl
Iterations
1 (10)
1 (16)
1 (15)
1 (15)
1 (16)
Total-Time
439.56
271.09
77.05
43.85
36.38
tmt_sym
Iterations
22 (17)
19 (7)
13 (18)
12 (16)
11 (18)
Total-Time
272.72
175.99
117.74
105.86
96.67
apache2
Iterations
8 (19)
8 (15)
7 (20)
7 (16)
5 (20)
Total-Time
189.46
105.15
74.75
66.76
50.49
cage14
Iterations
—
—
0 (10)
0 (10)
0 (11)
Total-Time
—
—
1,296.23
196.70
62.69
“—”: memory limits were reached in the initialization phase.Table 6
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for ecology2 problem.
ecology2
Cores/subdomains
64
128
256
512
1024
MPMSC
Pre-Time
20.09
6.60
3.50
2.69
2.51
Comm-Time
84.46
80.90
47.66
32.48
32.33
GMRES-Time
547.96
434.61
260.92
164.42
153.88
Total-Time
568.06
441.21
264.42
167.11
156.40
Iterations
34 (18)
35 (4)
23 (11)
14 (14)
14 (7)
Aggregates
19
39
62
139
263
Preconditioned Relative Residual
9.70
E
-
09
9.26
E
-
09
9.96
E
-
09
9.64
E
-
09
9.91
E
-
09
Relative Residual
9.10
E
-
08
6.99
E
-
08
8.67
E
-
08
7.19
E
-
08
1.32
E
-
07Table 7
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for thermal2 problem.
thermal2
Cores/subdomains
64
128
256
512
1024
MPMSC
Pre-Time
1491.13
198.63
51.66
10.12
5.51
Comm-Time
176.09
42.41
45.64
44.51
26.97
GMRES-Time
738.37
350.66
297.53
251.73
149.12
Total-Time
2229.50
549.29
349.19
261.84
154.64
Iterations
17 (12)
16 (18)
18 (6)
17 (16)
10 (16)
Aggregates
21
38
70
135
268
Preconditioned Relative Residual
9.98
E
-
09
9.61
E
-
09
1.00
E
-
08
9.97
E
-
09
9.87
E
-
09
Relative Residual
5.36
E
-
09
9.94
E
-
09
5.17
E
-
09
5.46
E
-
09
7.80
E
-
09Figure 8
The speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the ecology2 problem.Figure 9
The speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the thermal2 problem.For the ecology2 problem, the local coefficient matrices have reduced number of nonzero elements and hence Pre-Time is a small fraction of the Total-Time, not affecting the Total-Speedup curve, since theLU factorization requires reduced number of floating point operations. In contrast, for the thermal2 problem compared to ecology2 problem, the local coefficient matrices have more nonzero elements resulting in increased Pre-Time, since the LU factorization requires increased number of floating point operations, especially for small number of cores/subdomains. As the number of cores/subdomains increases, the Pre-Time decreases rapidly, because there are less nonzero elements in the local coefficient matrices resulting in improved performance of the LU factorization. Therefore, superlinear speedup is achieved for the thermal2 problem and not for the ecology2 problem.The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=100 and τ=10-6, and various numbers of cores, for the thermal2 problem, is presented in Figure 10. The preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner required 11 (2) iterations for convergence to the prescribed tolerance. The PPGMRES(m) based on MPMSC scheme performed better, in terms of total execution time, compared to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, especially for large numbers of cores-subdomains.Figure 10
The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=100 and τ=10-6, and various numbers of cores, for the thermal2 problem.Ιt should be noted that the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner for the ecology2 problem exceeds the available memory resources and does not converge to the prescribed tolerance for any compatible choice of p or τ not exceeding available memory resources.


## 3.1. 3D Poisson Problem

The PPGMRES(m) in conjunction with MPMSC scheme has been used for solving the Poisson equation in three space variables subject to Dirichlet boundary conditions. The Poisson equation was discretized with the 7-point stencil. The right-hand side of the linear systems was computed as the product of coefficient matrix A and the solution vector set to 0,1,…,n-1T, where n denotes the order of the linear system. In Table 1, the convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000) are presented for various numbers of cores-subdomains. The communication time is relative to the required iterations for convergence to the prescribed tolerance and the number of cores. It can be easily seen that as the number of cores and subdomains increases, the number of required iterations for convergence to the prescribed tolerance decreases for the 3D Poisson problem. In Table 2, the convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the Block Jacobi preconditioner, for the 3D Poisson problem (n=1,000,000) with various numbers of cores-subdomains are given. It should be stated that the MPMSC scheme presents better convergence behavior than the Block Jacobi preconditioning scheme and therefore requires less overall communication time, especially for many cores where the convergence behavior of the Block Jacobi method is downgraded. Therefore, it should be noted that the MPMSC scheme has better scalability than the Block Jacobi preconditioner. In Figure 3, the speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000) is depicted. In Figure 4, the weak scalability diagram of time versus cores for MPMSC scheme is presented. The order of the linear system increases as the number of cores increases since the local linear system should be of constant order. This results in augmenting, by 10,000 more unknowns, the linear system for every added core. The Total-Time should ideally remain constant. It should be noted that the Pre-Time is almost constant, as the initialization phase of the preconditioner is an inherently parallel process. From Figures 3 and 4 it is evident that the proposed method is scalable up to large number of cores, for the efficient solution of large linear systems. It should be mentioned that partitioning the domain Ω into more subdomains improves the convergence behavior of the method, since the grid corresponding to the aggregated (coarse) components becomes finer, which improves the accuracy of the local approximations to the solution. The improved convergence behavior reduces the communications and hence improves the performance of the proposed scheme. Due to the reduced communications, while the local computational work remains almost constant, it results in superlinear speedups, especially in the case of moderate number of subdomains. For small number of cores the Pre-Time is dominant compared to the GMRES-Time and as the preprocessing phase is an inherently parallel process, both methods have satisfactory speedup. When the number of cores is 1024, the Pre-Time is a small percentage of the Total-Time; therefore the speedup reduces, since it relies on the parallelization of the PPGMRES(m). It should be noted that the GMRES-Time increases proportionally to the problem size, as global communications, related to the synchronization of distributed vectors, and sequential redundant operations of PPGMRES(m) depend on the size of global vectors, compare [28], as it can be seen in Figure 4. However, the inherently parallel preprocessing phase scales, with satisfactory speedup, up to large number cores-subdomains. Thus, the Total-Time required by the proposed scheme is improved.Table 1
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000).
Poisson 3D
Cores/subdomains
64
128
256
512
1024
MPMSC
Pre-Time
367.38
70.04
21.45
6.31
3.64
Comm-Time
23.54
13.63
10.92
9.81
8.39
GMRES-Time
122.45
73.77
56.45
50.85
41.51
Total-Time
489.84
143.81
77.90
57.16
45.14
Iterations
4 (8)
4 (13)
4 (12)
4 (7)
3 (15)
Aggregates
13
20
40
84
162
Preconditioned Relative Residual
9.24
E
-
09
9.73
E
-
09
8.90
E
-
09
9.99
E
-
09
8.91
E
-
09
Relative Residual
1.38
E
-
08
1.48
E
-
08
1.28
E
-
08
1.54
E
-
08
1.47
E
-
08Table 2
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the Block Jacobi preconditioner, for the 3D Poisson problem (n=1,000,000).
Poisson 3D
Cores/subdomains
64
128
256
512
1024
Block Jacobi
Pre-Time
368.53
69.43
21.42
6.26
3.58
Comm-Time
24.14
16.63
12.21
11.72
12.62
GMRES-Time
133.28
90.17
66.13
64.25
61.81
Total-Time
501.81
159.60
87.55
70.51
65.39
Iterations
4 (18)
5 (15)
5 (9)
5 (15)
5 (12)
Preconditioned Relative Residual
9.52
E
-
09
9.33
E
-
09
9.81
E
-
09
9.85
E
-
09
9.90
E
-
09
Relative Residual
1.21
E
-
08
1.57
E
-
08
1.60
E
-
08
1.41
E
-
08
1.42
E
-
08Figure 3
The speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the 3D Poisson problem (n=1,000,000).Figure 4
The weak scalability (time versus cores) of the MPMSC scheme for 3D Poisson problem.Furthermore, the convergence behavior of MPMSC scheme for 3D Poisson problem with 10,000 unknowns per subdomain is shown in Figure5. The total number of iterations is computed as outer∗m+inner. It should be stated that the outer iterations of MPMSC scheme are almost constant for large number of subdomains. Moreover, the convergence behavior of the MPMSC scheme is related to the quality of the graph partitioning, provided by Metis, compare [23]. In Figure 6, the number of aggregates (naggs) of the MPMSC scheme against the number of subdomains for the 3D Poisson problem with 10,000 unknowns per subdomain is depicted.Figure 5
The convergence behavior of MPMSC scheme for 3D Poisson problem with various numbers of subdomains.Figure 6
The number of aggregated (coarse) components (naggs) of the MPMSC scheme for 3D Poisson problem with various numbers of subdomains.In Table3, the performance and convergence behavior of MPMSC scheme with distance two neighborhood for the 3D Poisson problem (n=1,000,000) with various numbers of cores-subdomains are presented. It should be mentioned that using distance two neighborhood in the subspace compression technique results in less aggregates compared to distance one neighborhood, which decreases the memory requirements. However, since the required number of iterations for convergence increases, the performance of the iterative scheme is downgraded.Table 3
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme with distance two neighborhood for the 3D Poisson problem (n=1,000,000).
Poisson 3D
Cores/subdomains
64
128
256
512
1024
MPMSCd=2
Pre-Time
368.53
69.43
21.42
6.26
3.58
Comm-Time
24.14
16.63
12.21
11.72
12.62
GMRES-Time
133.28
90.17
66.13
64.25
61.81
Total-Time
501.81
159.60
87.55
70.51
65.39
Iterations
4 (18)
5 (15)
5 (9)
5 (15)
5 (12)
Aggregates
6
9
19
32
51
Preconditioned Relative Residual
9.18
E
-
09
8.62
E
-
09
9.57
E
-
09
8.80
E
-
09
8.56
E
-
09
Relative Residual
1.57
E
-
08
1.18
E
-
08
1.59
E
-
08
1.38
E
-
08
1.46
E
-
08Furthermore, the parallel performance of the proposed method is compared to sequential preconditioned GMRES(m) in conjunction with dual threshold incomplete LU factorization (ILUT(p,τ)) preconditioner, compare [1]. The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=50 and τ=10-3, and various numbers of cores, for the 3D Poisson problem (n=1,000,000), is presented in Figure 7. The preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner required 2 (14) iterations for convergence to the prescribed tolerance. The PPGMRES(m) in conjunction with the MPMSC scheme performed better, in terms of total execution time, compared to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, especially for large numbers of cores-subdomains.Figure 7
The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=50 and τ=10-3, and various numbers of cores for the 3D Poisson problem (n=1,000,000).


## 3.2. UFL Matrices

In order to assess the applicability, scalability, and performance of MPMSC scheme for solving general large sparse linear systems, seven matrices from the University of Florida sparse matrix collection, compare [37], have been considered. Information concerning the aforementioned matrices is given in Table 4. The right-hand-side vectors for atmosmodd, atmosmodl, and thermal2 problems were given by the University of Florida sparse matrix collection, compare [37], whereas, for the rest of the problems, they were computed as the product of coefficient matrix A and the solution vector set to 0,1,…,n-1T.Table 4
Problems from the University of Florida sparse matrix collection.
Problem
Application
Size(A)
nnz(A)
atmosmodd
Atmospheric modeling
1270432
8814880
atmosmodl
1489752
10319760
tmt_sym
Electromagnetics problem
726713
5080961
apache2
Structural problem
715176
4817870
Cage14
DNA electrophoresis
1505785
27130349
ecology2
Ecology problem
999999
4995991
thermal2
Thermal problem
1228045
8580313The convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for the atmosmodd, atmosmodl, tmt_sym, apache2, and cage14 problems, are presented in Table 5. In Tables 6 and 7, the convergence behavior and the parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for the ecology2 and the thermal2 problems, are given. In Figures 8 and 9, the speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme for the ecology2 and the thermal2 problem, is shown.Table 5
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for various problems.
MPMSC
Cores/subdomains
64
128
256
512
1024
atmosmodd
Iterations
5 (12)
6 (12)
6 (10)
6 (17)
7 (10)
Total-Time
759.43
255.30
137.97
108.40
109.37
atmosmodl
Iterations
1 (10)
1 (16)
1 (15)
1 (15)
1 (16)
Total-Time
439.56
271.09
77.05
43.85
36.38
tmt_sym
Iterations
22 (17)
19 (7)
13 (18)
12 (16)
11 (18)
Total-Time
272.72
175.99
117.74
105.86
96.67
apache2
Iterations
8 (19)
8 (15)
7 (20)
7 (16)
5 (20)
Total-Time
189.46
105.15
74.75
66.76
50.49
cage14
Iterations
—
—
0 (10)
0 (10)
0 (11)
Total-Time
—
—
1,296.23
196.70
62.69
“—”: memory limits were reached in the initialization phase.Table 6
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for ecology2 problem.
ecology2
Cores/subdomains
64
128
256
512
1024
MPMSC
Pre-Time
20.09
6.60
3.50
2.69
2.51
Comm-Time
84.46
80.90
47.66
32.48
32.33
GMRES-Time
547.96
434.61
260.92
164.42
153.88
Total-Time
568.06
441.21
264.42
167.11
156.40
Iterations
34 (18)
35 (4)
23 (11)
14 (14)
14 (7)
Aggregates
19
39
62
139
263
Preconditioned Relative Residual
9.70
E
-
09
9.26
E
-
09
9.96
E
-
09
9.64
E
-
09
9.91
E
-
09
Relative Residual
9.10
E
-
08
6.99
E
-
08
8.67
E
-
08
7.19
E
-
08
1.32
E
-
07Table 7
Convergence behavior and parallel performance of the PPGMRES(m) method, in conjunction with the MPMSC scheme for thermal2 problem.
thermal2
Cores/subdomains
64
128
256
512
1024
MPMSC
Pre-Time
1491.13
198.63
51.66
10.12
5.51
Comm-Time
176.09
42.41
45.64
44.51
26.97
GMRES-Time
738.37
350.66
297.53
251.73
149.12
Total-Time
2229.50
549.29
349.19
261.84
154.64
Iterations
17 (12)
16 (18)
18 (6)
17 (16)
10 (16)
Aggregates
21
38
70
135
268
Preconditioned Relative Residual
9.98
E
-
09
9.61
E
-
09
1.00
E
-
08
9.97
E
-
09
9.87
E
-
09
Relative Residual
5.36
E
-
09
9.94
E
-
09
5.17
E
-
09
5.46
E
-
09
7.80
E
-
09Figure 8
The speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the ecology2 problem.Figure 9
The speedup of the PPGMRES(m) method, in conjunction with the MPMSC scheme, for the thermal2 problem.For the ecology2 problem, the local coefficient matrices have reduced number of nonzero elements and hence Pre-Time is a small fraction of the Total-Time, not affecting the Total-Speedup curve, since theLU factorization requires reduced number of floating point operations. In contrast, for the thermal2 problem compared to ecology2 problem, the local coefficient matrices have more nonzero elements resulting in increased Pre-Time, since the LU factorization requires increased number of floating point operations, especially for small number of cores/subdomains. As the number of cores/subdomains increases, the Pre-Time decreases rapidly, because there are less nonzero elements in the local coefficient matrices resulting in improved performance of the LU factorization. Therefore, superlinear speedup is achieved for the thermal2 problem and not for the ecology2 problem.The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=100 and τ=10-6, and various numbers of cores, for the thermal2 problem, is presented in Figure 10. The preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner required 11 (2) iterations for convergence to the prescribed tolerance. The PPGMRES(m) based on MPMSC scheme performed better, in terms of total execution time, compared to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, especially for large numbers of cores-subdomains.Figure 10
The speedup of PPGMRES(m) in conjunction with MPMSC scheme with respect to the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner, with p=100 and τ=10-6, and various numbers of cores, for the thermal2 problem.Ιt should be noted that the preconditioned GMRES(m) in conjunction with the ILUT(p,τ) preconditioner for the ecology2 problem exceeds the available memory resources and does not converge to the prescribed tolerance for any compatible choice of p or τ not exceeding available memory resources.


## 4. Conclusion

The proposed scheme, namely, multiprojection method with subspace compression, has been shown to be effective in solving various problems in distributed memory parallel systems. Further, the PPGMRES(m) in conjunction with the MPMSC scheme scales up to a large number of cores exhibiting superlinear speedups, especially for problems that require large computational work in the initialization phase of the preconditioning scheme. The proposed scheme has improved convergence behavior as the number of subdomains increases compared to the extant domain decomposition based preconditioners. Future work will be concentrated towards using asynchronous and neighboring communications between heterogeneous compute nodes, in order to assess the performance of the method in modern cloud environments. Moreover, the choice of the most suitable parallel graph partitioning algorithm, according to the performance and the partitioning quality, will be studied.



---

*Source: 2580820-2017-07-30.xml*