# Parallel Rayleigh Quotient Optimization with FSAI-Based Preconditioning

Authors: Luca Bergamaschi; Angeles Martínez; Giorgio Pini
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

The present paper describes a parallel preconditioned algorithm for the solution of partial eigenvalue problems for large sparse symmetric matrices, on parallel computers. Namely, we consider the Deflation-Accelerated Conjugate Gradient (DACG) algorithm accelerated by factorized-sparse-approximate-inverse- (FSAI-) type preconditioners. We present an enhanced parallel implementation of the FSAI preconditioner and make use of the recently developed Block FSAI-IC preconditioner, which combines the FSAI and the Block Jacobi-IC preconditioners. Results onto matrices of large size arising from finite element discretization of geomechanical models reveal that DACG accelerated by these type of preconditioners is competitive with respect to the available public parallelhypre package, especially in the computation of a few of the leftmost eigenpairs. The parallel DACG code accelerated by FSAI is written in MPI-Fortran 90 language and exhibits good scalability up to one thousand processors.


---

## Body

## 1. Introduction

The computation by iterative methods of thes partial eigenspectrum of the generalized eigenproblem:
(1.1)Au=λBu,
where A,B∈ℝn×n are large sparse symmetric positive definite (SPD) matrices, is an important and difficult task in many applications. It has become increasingly widespread owing to the development in the last twenty years of robust and computationally efficient schemes and corresponding software packages. Among the most well-known approaches for the important class of symmetric positive definite (SPD) matrices are the implicitly restarted Arnoldi method (equivalent to the Lanczos technique for this type of matrices) [1, 2], the Jacobi-Davidson (JD) algorithm [3], and schemes based on preconditioned conjugate gradient minimization of the Rayleigh quotient [4, 5].The basic idea of the latter is to minimize the Rayleigh Quotient(1.2)q(x)=xTAxxTBx.
in a subspace which is orthogonal to the previously computed eigenvectors via a preconditioned CG-like procedure. Among the different variants of this technique we chose to use the Deflation-Accelerated Conjugate Gradient (DACG) scheme [4, 6] which has been shown to be competitive with the Jacobi Davidson method and with the PARPACK package [7]. As in any other approach, for our DACG method, the choice of the preconditioning technique is a key factor to accelerate and, in some cases even to allow for, convergence. To accelerate DACG in a parallel environment we selected the Factorized Sparse Approximate inverse (FSAI) preconditioner introduced in [8]. We have developed a parallel implementation of this algorithm which has displayed excellent performances on both the setup phase and the application phase within a Krylov subspace solver [9–11]. The effectiveness of the FSAI preconditioner in the acceleration of DACG is compared to that of the Block FSAI-IC preconditioner, recently developed in [12], which combines the FSAI and the Block Jacobi-IC preconditioners obtaining good results on a small number of processors for the solution of SPD linear systems and for the solution of large eigenproblems [13]. We used the resulting parallel codes to compute a few of the leftmost eigenpairs of a set of test matrices of large size arising from Finite Element discretization of geomechanical models. The reported results show that DACG preconditioned with either FSAI or BFSAI is a scalable and robust algorithm for the partial solution of SPD eigenproblems. The parallel performance of DACG is also compared to that of the publicly available parallel package hypre [14] which implements a number of preconditioners which can be used in combination with the Locally Optimal Block PCG (LOBPCG) iterative eigensolver [15]. The results presented in this paper show that the parallel DACG code accelerated by FSAI exhibits good scalability up to one thousand processors and displays comparable performance with respect to hypre, specially when a low number of eigenpairs is sought.The outline of the paper is as follows: in Section2 we describe the DACG Algorithm; in Sections 3 and 4 we recall the definition and properties of the FSAI and BFSAI preconditioners, respectively. Section 5 contains the numerical results obtained with the proposed algorithm in the eigensolution of very large SPD matrices of size up to almost 7 million unknowns and 3×108 nonzeros. A comparison with the hypre eigensolver code is also included. Section 6 ends the paper with some conclusions.


## 2. The DACG Iterative Eigensolver and Implementation

The DACG algorithm sequentially computes the eigenpairs, starting from the leftmost one(λ1,u1). To evaluate the jth eigenpair, j>1, DACG minimizes the Rayleigh Quotient (RQ) in a subspace orthogonal to the j-1 eigenvectors previously computed. More precisely, DACG minimizes the Rayleigh Quotient:
(2.1)q(z)=zTAzzTz,
where
(2.2)z=x-Uj(UjTx),Uj=[u1,…,uj-1],x∈Rn.
The first eigenpair (λ1,u1) is obtained by minimization of (2.1) with z=x(U1=∅). Indicating with M the preconditioning matrix, that is, M≈A-1, the s leftmost eigenpairs are computed by the conjugate gradient procedure [6] described in Algorithm 1.Algorithm 1:DACG Algorithm.
Choose toleranceε, set U=0.DOj=1,s(1) Choose x0 such that UTx0=0; set k=0, β0=0;(2)  x0A=Ax0,γ=x0Tx0A,η=x0Tx0,q0≡q(x0)=γ/η,r0=x0A-q0x0;(3) REPEAT(3.1)gk≡∇q(xk)=(2/η)rk;(3.2)gkM=Mgk;(3.3)IFk>0THENβk=gkT(gkM-gk-1M)/gk-1Tgk-1M;(3.4)p̃k=gkM+βkpk-1;(3.5)pk=p̃k-U(UTp̃k)(3.6)pkA=Apk,(3.7)αk=argmint{q(xk+tpk)}=(ηd-γb+Δ)/2(bc-ad), witha=pkTxkA,b=pkTpkA,c=pkTxk,d=pkTpk,Δ=(ηd-γb)2-4(bc-ad)(γa-ηc);(3.8)xk+1=xk+αkpk,xk+1A=xkA+αkpkA;(3.9)γ=γ+2aαk+bαk2,η=η+2cαk+dαk2;(3.10)qk+1≡q(xk+1)=γ/η;(3.11)k=k+1;(3.12)rk=xkA-qkxkUNTIL(qk+1-qk)/qk+1<tol;(4)  λj=qk,uj=xk/η,U=[U,uj].END DOThe schemes relying on the Rayleigh quotient optimization are quite attractive for parallel computations; however preconditioning is an essential feature to ensure practical convergence. When seeking for an eigenpair(λj,uj) it can be proved that the number of iterations is proportional to the square root of the condition number ξj=κ(Hj) of the Hessian of the Rayleigh quotient in the stationary point uj [4]. It turns out that Hj is similar to (A-λjI)M which is not SPD. However, Hj operates on the orthogonal space spanned by the previous eigenvectors, so that the only important eigenvalues are the positive ones. In the non-preconditioned case (i.e., M=I) we would have
(2.3)κ(Hj)≈λNλj+1-λj.
where in the ideal case M≡A-1, we have
(2.4)κ(Hj)≈λjλj+1-λj=ξj≪λNλj+1-λj.
Therefore, even though A-1 is not the optimal preconditioner for A-λjI, however, if M is a good preconditioner of A then the condition number κ(Hj) will approach ξj.


## 3. The FSAI Preconditioner

The FSAI preconditioner, initially proposed in [8, 16], has been later developed and implemented in parallel by Bergamaschi and Martínez in [9]. Here, we only shortly recall the main features of this preconditioner. Given an SPD matrix A the FSAI preconditioner approximately factorizes its inverse as a product of two sparse triangular matrices as
(3.1)A-1≈WTW.
The choice of nonzeros in W is based on a sparsity pattern which in our work may be the same as Ãd where Ã is the result of prefiltration [10] of A, that is, dropping of all elements below of a threshold parameter δ. The entries of W are computed by minimizing the Frobenius norm of I-WL, where L is the exact Cholesky factor of A, without forming explicitly the matrix L. The computed W is then sparsified by dropping all the elements which are below a second tolerance parameter (ε). The final FSAI preconditioner is therefore related to the following three parameters: δ, prefiltration threshold; d, power of A generating the sparsity pattern (we allow d∈{1,2,4} in our experiments); ε, postfiltration threshold.


### 3.1. Parallel Implementation of FSAI-DACG

We have developed a parallel code written in FORTRAN 90 and which exploits the MPI library for exchanging data among the processors. We used a block row distribution of all matrices (A, W, and WT), that is, with complete rows assigned to different processors. All these matrices are stored in static data structures in CSR format.Regarding the preconditioner computation, we stress that any rowi of matrix W of FSAI preconditioner is computed independently of each other, by solving a small SPD dense linear system of size ni equal to the number of nonzeros allowed in row i of W. Some of the rows which contribute to form this linear system may be nonlocal to processor i and should be received from other processors. To this aim we implemented a routine called get_extra_rows which carries out all the row exchanges among the processors, before starting the computation of W, which proceed afterwards entirely in parallel. Since the number of nonlocal rows needed by each processor is relatively small we chose to temporarily replicate these rows on auxiliary data structures. Once W is obtained a parallel transposition routine provides every processor with its part of WT.The DACG iterative solver is essentially based on scalar and matrix-vector products. We made use of an optimized parallel matrix-vector product which has been developed in [17] showing its effectiveness up to 1024 processors.


## 3.1. Parallel Implementation of FSAI-DACG

We have developed a parallel code written in FORTRAN 90 and which exploits the MPI library for exchanging data among the processors. We used a block row distribution of all matrices (A, W, and WT), that is, with complete rows assigned to different processors. All these matrices are stored in static data structures in CSR format.Regarding the preconditioner computation, we stress that any rowi of matrix W of FSAI preconditioner is computed independently of each other, by solving a small SPD dense linear system of size ni equal to the number of nonzeros allowed in row i of W. Some of the rows which contribute to form this linear system may be nonlocal to processor i and should be received from other processors. To this aim we implemented a routine called get_extra_rows which carries out all the row exchanges among the processors, before starting the computation of W, which proceed afterwards entirely in parallel. Since the number of nonlocal rows needed by each processor is relatively small we chose to temporarily replicate these rows on auxiliary data structures. Once W is obtained a parallel transposition routine provides every processor with its part of WT.The DACG iterative solver is essentially based on scalar and matrix-vector products. We made use of an optimized parallel matrix-vector product which has been developed in [17] showing its effectiveness up to 1024 processors.


## 4. Block FSAI-IC Preconditioning

The Block FSAI-IC preconditioner, BFSAI-IC in the following, is a recent development for the parallel solution to Symmetric Positive Definite (SPD) linear systems. Assume thatD is an arbitrary nonsingular block diagonal matrix consisting of nb equal size blocks.Let𝒮L and 𝒮BD be a sparse lower triangular and a dense block diagonal nonzero pattern, respectively, for an n×n matrix. Even though not strictly necessary, for the sake of simplicity assume that 𝒮BD consists of nb diagonal blocks with equal size m=n/nb and let D∈ℝn×n be an arbitrary full-rank matrix with nonzero pattern 𝒮BD.Consider the set of lower block triangular matricesF with a prescribed nonzero pattern SBL and minimize over F the Frobenius norm:
(4.1)‖D-FL‖F,
where L is the exact lower Cholesky factor of an SPD matrix A. A matrix F satisfying the minimality condition (4.1) for a given D is the lower block triangular factor of BFSAI-IC. Recalling the definition of the classical FSAI preconditioner, it can be noticed that BFSAI-IC is a block generalization of the FSAI concept.The differentiation of (4.1) with respect to the unknown entries [F]ij, (i,j)∈SBL, yields the solution to n independent dense subsystems which, in the standard FSAI case, do not require the explicit knowledge of L. The effect of applying F to A is to concentrate the largest entries of the preconditioned matrix FAFT into nb diagonal blocks. However, as D is arbitrary, it is still not ensured that FAFT is better than A in an iterative method, so it is necessary to precondition FAFT again. As FAFT resembles a block diagonal matrix, an efficient technique relies on using a block diagonal matrix which collects an approximation of the inverse of each diagonal block Bib of FAFT.It is easy to show thatF is guaranteed to exist with SPD matrices and Bib is SPD, too [12]. Using an IC decomposition with partial fill-in for each block Bib and collecting in J the lower IC factors, the resulting preconditioned matrix reads
(4.2)J-1FAFTJ-T=WAWT
with the final preconditioner
(4.3)M=WTW=FTJ-TJ-1F.M in (4.3) is the BFSAI-IC preconditioner of A.For its computation BFSAI-IC needs the selection ofnb and 𝒮L¯. The basic requirement for the number of blocks nb is to be larger than or equal to the number of computing cores p. From a practical viewpoint, however, the most efficient choice in terms of both wall clock time and iteration count is to keep the blocks as large as possible, thus implying nb=p. Hence, nb is by default set equal to p. By distinction, the choice of 𝒮L¯ is theoretically more challenging and still not completely clear. A widely accepted option for other approximate inverses, such as FSAI or SPAI, is to select the nonzero pattern of Ad for small values of d on the basis of the Neumann series expansion of A-1. Using a similar approach, in the BFSAI construction we select 𝒮L¯ as the lower block triangular pattern of Ad. As the nonzeros located in the diagonal blocks are not used for the computation of F a larger value of d, say 3 or 4, can still be used.Though theoretically not necessary, three additional user-specified parameters are worth introducing in order to better control the memory occupation and the BFSAI-IC density:(1)
ε is a postfiltration parameter that allows for dropping the smallest entries of F. In particular, [F]ij is neglected if [F]ij<ε∥fi∥2, where fi is the ith row of F;(2)
ρB is a parameter that controls the fill-in of Bib and determines the maximum allowable number of nonzeros for each row of Bib in addition to the corresponding entries of A. Quite obviously, the largest ρB entries only are retained;(3)
ρL is a parameter that controls the fill-in of each IC factor L̃ib denoting the maximum allowable number of nonzeros for each row of L̃ib in addition to the corresponding entries of Bib.An OpenMP implementation of the algorithms above is available in [18].


## 5. Numerical Results

In this section we examine the performance of the parallel DACG preconditioned by both FSAI and BFSAI in the partial solution of four large-size sparse eigenproblems. The test cases, which we briefly describe below, are taken from different real engineering mechanical applications. In detail, they are as follows.(i)
FAULT-639 is obtained from a structural problem discretizing a faulted gas reservoir with tetrahedral finite elements and triangular interface elements [19]. The interface elements are used with a penalty formulation to simulate the faults behavior. The problem arises from a 3D discretization with three displacement unknowns associated to each node of the grid.(ii)
PO-878 arises in the simulation of the consolidation of a real gas reservoir of the Po Valley, Italy, used for underground gas storage purposes (for details, see [20]).(iii)
GEO-1438 is obtained from a geomechanical problem discretizing a region of the earth crust subject to underground deformation. The computational domain is a box with an areal extent of 50 × 50 km and 10 km deep consisting of regularly shaped tetrahedral finite elements. The problem arises from a 3D discretization with three displacement unknowns associated to each node of the grid [21].(iv)
CUBE-6091 arises from the equilibrium of a concrete cube discretized by a regular unstructured tetrahedral grid.Matrices FAULT-639 and GEO-1438 are publicly available in the University of Florida Sparse Matrix Collection athttp://www.cise.ufl.edu/research/sparse/matrices/.In Table1 we report sizes and nonzeros of the four matrices together with three of the most significant eigenvalues for each problem.Table 1
Size, number of nonzeros, and three representative eigenvalues of the test matrices.
SizeNonzerosλ1λ10λNFAULT-639638,80228,614,5646.99·1061.73·1072.52·1016PO-878878,35538,847,9151.46·1064.45·1065.42·1015GEO-14381,437,96063,156,6907.81·1051.32·1061.11·1013CUBE-60916,091,008270,800,5861.82·1013.84·1021.05·1007The computational performance of FSAI is compared to the one obtained by using BFSAI as implemented in [12]. The comparison is done evaluating the number of iterations niter to converge at the same tolerance, the wall clock time in seconds Tprec and Titer for the preconditioner computation, and the eigensolver to converge, respectively, with the total time Ttot=Tprec+Titer. All tests are performed on the IBM SP6/5376 cluster at the CINECA Centre for High Performance Computing, equipped with IBM Power6 processors at 4.7 GHz with 168 nodes, 5376 computing cores, and 21 Tbyte of internal network RAM. The FSAI-DACG code is written in Fortran 90 and compiled with -O4 -q64 -qarch=pwr6 -qtune=pwr6 -qnoipa -qstrict -bmaxdata:0x70000000 options. For the BFSAI-IC code only an OpenMP implementation is presently available.To study parallel performance we will use a strong scaling measure to see how the CPU times vary with the number of processors for a fixed total problem size. Denote withTp the total CPU elapsed times expressed in seconds on p processors. We introduce a relative measure of the parallel efficiency achieved by the code, Sp(p̅), which is the pseudo speedup computed with respect to the smallest number of processors (p̅) used to solve a given problem. Accordingly, we will denote by Ep(p̅) the corresponding efficiency:
(5.1)Sp(p̅)=Tp̅p̅Tp,Ep(p̅)=Sp(p̅)p=Tp̅p̅Tpp.


### 5.1. FSAI-DACG Results

In this section we report the results of our FSAI-DACG implementation in the computation of the 10 leftmost eigenpairs of the 4 test problems. We used the exit test described in the DACG algorithm (see Algorithm1) with tol=10-10. The results are summarized in Table 2. As the FSAI parameters, we choose δ=0.1, d=4, and ε=0.1 for all the test matrices. This combination of parameters produces, on the average, the best (or close to the best) performance of the iterative procedure. Note that the number of iterations does not change with the number of processors, for a fixed problem. The scalability of the code is very satisfactory in both the setup stage (preconditioner computation) and the iterative phase.Table 2
Number of iterations, timings, and scalability indices for FSAI-DACG in the computation of the 10 leftmost eigenpairs of the four test problems.
piterTprecTiterTtotSp(4)Ep(4)FAULT-6394444825.9261.4287.38444813.2139.0152.27.60.941644486.669.476.015.10.953244484.028.232.235.71.116444481.915.517.466.11.0312844481.19.410.5109.00.85PO-8784587648.1722.5770.68587625.2399.8425.07.30.9116587611.4130.2141.621.81.363258766.865.872.542.51.336458764.130.134.190.31.4112858761.919.121.0146.81.15GEO-14374621690.3901.5991.78621647.5478.9526.47.50.9416621624.7239.4264.115.00.9432621613.6121.0134.629.50.926462168.260.969.157.40.9012862164.229.533.8117.50.9225662162.319.121.4185.40.72piterTprecTiterTtotSp(16)Ep(16)CUBE-60911615796121.52624.82746.2321579662.21343.81406.031.30.98641579632.5737.0769.557.10.891281579617.3388.4405.7108.30.85256157969.1183.9192.9227.80.89512157965.7106.0111.7393.50.771024157963.876.680.4546.60.53


### 5.2. BFSAI-IC-DACG Results

We present in this section the results of DACG accelerated by the BFSAI-IC preconditioner for the approximation of thes=10 leftmost eigenpairs of the matrices described above.Table3 provides iteration count and total CPU time for BFSAI-DACG with different combinations of the parameters needed to construct the BFSAI-IC preconditioner for matrix PO-878 and using from 2 to 8 processors. It can be seen from Table 3 that the assessment of the optimal parameters, ε, ρB, and ρL, is not an easy task, since the number of iterations may highly vary depending on the number of processors. We chose in this case the combination of parameters producing the second smallest total time with p=2,4,8 processors. After intensive testing for all the test problems, we selected similarly the “optimal” values which are used in the numerical experiments reported in Table 4:(i)
FAULT-639:d=3, ε=0.05, ρB=10, ρL=60,(ii)
PO-878d=3, ε=0.05, ρB=10, ρL=10,(iii)
GEO-1438:d=3, ε=0.05, ρB=10, ρL=50,(iv)
CUBE-6091:d=2, ε=0.01, ρB=0, ρL=10.Table 3
Performance of BFSAI-DACG for matrix PO-878 with 2 to 8 processors and different parameter values.
dρBερLp=2p=4p=8iterTtotiterTtotiterTtot2100.01102333385.762877286.243753273.772100.05102345415.812803245.423815142.932100.05202186370.862921276.413445257.182100.00102328445.162880241.233392269.412200.05102340418.202918224.323720253.983100.05102122375.172638228.393366149.593100.05201946433.042560304.433254263.513100.05301822411.002481321.303176179.673100.05401729439.472528346.823019188.134100.05102035499.452469350.033057280.31Table 4
Number of iterations for BFSAI-DACG in the computations of the 10 leftmost eigenpairs.
Matrixnb2481632641282565121024FAULT-6391357143415942002305333363553PO-8782122263833664157482851545373GEO-143814581797211327783947464748504996CUBE-6091585765577746860894439996101899965The user-specified parameters for BFSAI-IC given above provide evidence that it is important to build a dense preconditioner based on the lower nonzero pattern ofA3 (except for CUBE-6091, which is built on a regular discretization) with the aim at decreasing the number of DACG iterations. Anyway, the cost for computing such a dense preconditioner appears to be almost negligible with respect to the wall clock time needed to iterate to convergence.We recall that, presently, the code BFSAI-IC is implemented in OpenMP, and the results in terms of CPU time are significant only forp≤8. For this reason the number of iterations reported in Table 4 is obtained with increasing number of blocks nb and with p=8 processors. This iteration number accounts for a potential implementation of BFSAI-DACG under the MPI (or hybrid OpenMP-MPI) environment as the number of iterations depends only on the number of blocks, irrespective of the number of processors.The only meaningful comparison between FSAI-DACG and BFSAI-DACG can be carried out in terms of iteration numbers which are smaller for BFSAI-DACG for a small number of processors. The gap between FSAI and BFSAI iterations reduces when the number of processors increases.


### 5.3. Comparison with the LOBPCG Eigensolver Provided byhypre

In order to validate the effectiveness of our preconditioning in the proposed DACG algorithm with respect to already available public parallel eigensolvers, the results given in Tables2 and 4 are compared with those obtained by the schemes implemented in the hypre software package [14]. The Locally Optimal Block Preconditioned Conjugate Gradient method (LOBPCG) [15] is experimented with, using the different preconditioners developed in the hypre project, that is, algebraic multigrid (AMG), diagonal scaling (DS), approximate inverse (ParaSails), additive Schwarz (Schwarz), and incomplete LU (Euclid). The hypre preconditioned CG is used for the inner iterations within LOBPCG. For details on the implementation of the LOBPCG algorithm, see, for instance, [22]. The selected preconditioner, ParaSails, is on its turn based on the FSAI preconditioner, so that the different FSAI-DACG and ParaSails-LOBPCG performances should be ascribed mainly to the different eigensolvers rather than to the preconditioners.We first carried out a preliminary set of runs with the aim of assessing the optimal value of the block sizebl parameter, that is, the size of the subspace where to seek for the eigenvectors. Obviously it must be bl ≥ s = 10. We fixed to 16 the number of processors and obtained the results summarized in Table 5 with different values of bl  ∈  [10,15]. We found that, only in problem CUBE-6091, a value of bl larger than 10, namely, bl = 12, yields an improvement in the CPU time. Note that we also made this comparison with different number of processors, and we obtained analogous results.Table 5
Iterations and CPU time for the iterative solver of LOBPCG-hypre preconditioned by Parasails with different values of bl and p=16 processors.
Matrixbl = 10bl = 11bl = 12bl = 15iterTiteriterTiteriterTiteriterTiterFAULT-63915679.515785.315796.1160128.1PO-87845117.041131.638151.335192.6GEO-143823123.772173.730152.5121291.1CUBE-60911011670.51432414.0381536.7351680.9Table6 presents the number of iterations and timings using the LOBPCG algorithm in the hypre package. The LOBPCG wall clock time is obtained with the preconditioner allowing for the best performance in the specific problem at hand, that is, ParaSails for all the problems. Using AMG as the preconditioner did not allow for convergence in three cases out of four, with the only exception of the FAULT-639 problem, in which the CPU timings were however very much larger than using ParaSails.Table 6
Number of iterations, timings, and scalability of LOBPCG-hypre preconditioned by Parasails.
piterTprecTiterTtotSp(4)Ep(4)FAULT-639 pcgitr =541552.5331.2333.781561.3167.6168.97.90.99161560.879.580.316.61.04321500.538.839.334.01.06641450.322.222.559.40.931281570.114.814.989.70.70PO-878 pcgitr = 304453.3438.4441.78501.3232.3234.07.60.9416451.0117.0118.015.00.9432450.763.263.927.60.8664470.434.434.850.80.79128410.319.4419.7489.50.70GEO-1438pcgitr = 304267.7478.0485.78224.0256.8260.87.50.9316232.1123.7125.815.40.9632281.273.174.326.20.8264230.835.536.353.50.84128250.520.320.893.20.73256260.312.913.2147.20.57piterTprecTiterTtotSp(16)Ep(16)CUBE-609116389.21536.71545.932364.7807.5812.230.50.9564383.2408.2411.460.10.94128411.6251.4253.097.80.76256350.9105.9106.8231.60.90512390.665.365.9375.30.731024370.337.738.0650.90.64All matrices have to be preliminarily scaled by their maximum coefficient in order to allow for convergence. To make the comparison meaningful, the outer iterations of the different methods are stopped when the average relative error measure of the computed leftmost eigenpairs gets smaller than10-10, in order to obtain a comparable accuracy as in the other codes. We also report in Table 6 the number of inner preconditioned CG iterations (pcgitr).To better compare our FSAI DACG with the LOBPCG method, we depict in Figure1 the total CPU time versus the number of processor for the two codes. FSAI-DACG and LOBPCG provide very similar scalability, being the latter code a little bit more performing on the average. On the FAULT-639 problem, DACG reveals faster than LOBPCG, irrespective of the number of processors employed.Comparison between FSAI-DACG and LOBPCG-hypre in terms of total CPU time for different number of processors.
(a)(b)(c)(d)Finally, we have carried out a comparison of the two eigensolvers in the computation of only the leftmost eigenpair. Differently from LOBPCG, which performs a simultaneous approximation of all the selected eigenpairs, DACG proceeds in the computation of the selected eigenpairs in a sequential way. For this reason, DACG should be the better choice, at least in principle, when just one eigenpair is sought. We investigate this feature, and the results are summarized in Table7. We include the total CPU time and iteration count needed by LOBPCG and FSAI-DACG to compute the leftmost eigenpair with 16 processors. For the LOBPCG code we report only the number of outer iterations.Table 7
Performance of LOBPCG-hypre with Parasails and 16 processors in the computation of the smallest eigenvalue using bl = 1 and bl = 2 and FSAI-DACG.
LOBPCG,bl = 1LOBPCG,bl = 2FSAI-DACGiterTtotiterTtotiterTtotFAULT-63914410.113217.5103015.4PO-8789943.23429.199320.4GEO-14935540.22637.375427.0CUBE-60911445218.858522.43257561.1The parameters used to construct the FSAI preconditioner for these experiments are as follows:(1)
FAULT-639.δ=0.1, d=2, ε=0.05,(2)
PO-878.δ=0.2, d=4, ε=0.1,(3)
GEO-1438.δ=0.1, d=2, ε=0.1,(4)
CUBE-6091.δ=0.0, d=1, ε=0.05.These parameters differ from those employed to compute the FSAI preconditioner in the assessment of the 10 leftmost eigenpairs and have been selected in order to produce a preconditioner relatively cheap to compute. This is so because otherwise the setup time would prevail over the iteration time. Similarly, to compute just one eigenpair with LOBPCG we need to setup a different value for pcgitr, the number of inner iterations. As it can be seen from Table7, in the majority of the test cases, LOBPCG takes less time to compute 2 eigenpairs than just only 1. FSAI-DACG reveals more efficient than the best LOBPCG on problems PO-878 and GEO-1438. On the remaining two problems the slow convergence exhibited by DACG is probably due to the small relative separation ξ1 between λ1 and λ2.


## 5.1. FSAI-DACG Results

In this section we report the results of our FSAI-DACG implementation in the computation of the 10 leftmost eigenpairs of the 4 test problems. We used the exit test described in the DACG algorithm (see Algorithm1) with tol=10-10. The results are summarized in Table 2. As the FSAI parameters, we choose δ=0.1, d=4, and ε=0.1 for all the test matrices. This combination of parameters produces, on the average, the best (or close to the best) performance of the iterative procedure. Note that the number of iterations does not change with the number of processors, for a fixed problem. The scalability of the code is very satisfactory in both the setup stage (preconditioner computation) and the iterative phase.Table 2
Number of iterations, timings, and scalability indices for FSAI-DACG in the computation of the 10 leftmost eigenpairs of the four test problems.
piterTprecTiterTtotSp(4)Ep(4)FAULT-6394444825.9261.4287.38444813.2139.0152.27.60.941644486.669.476.015.10.953244484.028.232.235.71.116444481.915.517.466.11.0312844481.19.410.5109.00.85PO-8784587648.1722.5770.68587625.2399.8425.07.30.9116587611.4130.2141.621.81.363258766.865.872.542.51.336458764.130.134.190.31.4112858761.919.121.0146.81.15GEO-14374621690.3901.5991.78621647.5478.9526.47.50.9416621624.7239.4264.115.00.9432621613.6121.0134.629.50.926462168.260.969.157.40.9012862164.229.533.8117.50.9225662162.319.121.4185.40.72piterTprecTiterTtotSp(16)Ep(16)CUBE-60911615796121.52624.82746.2321579662.21343.81406.031.30.98641579632.5737.0769.557.10.891281579617.3388.4405.7108.30.85256157969.1183.9192.9227.80.89512157965.7106.0111.7393.50.771024157963.876.680.4546.60.53


## 5.2. BFSAI-IC-DACG Results

We present in this section the results of DACG accelerated by the BFSAI-IC preconditioner for the approximation of thes=10 leftmost eigenpairs of the matrices described above.Table3 provides iteration count and total CPU time for BFSAI-DACG with different combinations of the parameters needed to construct the BFSAI-IC preconditioner for matrix PO-878 and using from 2 to 8 processors. It can be seen from Table 3 that the assessment of the optimal parameters, ε, ρB, and ρL, is not an easy task, since the number of iterations may highly vary depending on the number of processors. We chose in this case the combination of parameters producing the second smallest total time with p=2,4,8 processors. After intensive testing for all the test problems, we selected similarly the “optimal” values which are used in the numerical experiments reported in Table 4:(i)
FAULT-639:d=3, ε=0.05, ρB=10, ρL=60,(ii)
PO-878d=3, ε=0.05, ρB=10, ρL=10,(iii)
GEO-1438:d=3, ε=0.05, ρB=10, ρL=50,(iv)
CUBE-6091:d=2, ε=0.01, ρB=0, ρL=10.Table 3
Performance of BFSAI-DACG for matrix PO-878 with 2 to 8 processors and different parameter values.
dρBερLp=2p=4p=8iterTtotiterTtotiterTtot2100.01102333385.762877286.243753273.772100.05102345415.812803245.423815142.932100.05202186370.862921276.413445257.182100.00102328445.162880241.233392269.412200.05102340418.202918224.323720253.983100.05102122375.172638228.393366149.593100.05201946433.042560304.433254263.513100.05301822411.002481321.303176179.673100.05401729439.472528346.823019188.134100.05102035499.452469350.033057280.31Table 4
Number of iterations for BFSAI-DACG in the computations of the 10 leftmost eigenpairs.
Matrixnb2481632641282565121024FAULT-6391357143415942002305333363553PO-8782122263833664157482851545373GEO-143814581797211327783947464748504996CUBE-6091585765577746860894439996101899965The user-specified parameters for BFSAI-IC given above provide evidence that it is important to build a dense preconditioner based on the lower nonzero pattern ofA3 (except for CUBE-6091, which is built on a regular discretization) with the aim at decreasing the number of DACG iterations. Anyway, the cost for computing such a dense preconditioner appears to be almost negligible with respect to the wall clock time needed to iterate to convergence.We recall that, presently, the code BFSAI-IC is implemented in OpenMP, and the results in terms of CPU time are significant only forp≤8. For this reason the number of iterations reported in Table 4 is obtained with increasing number of blocks nb and with p=8 processors. This iteration number accounts for a potential implementation of BFSAI-DACG under the MPI (or hybrid OpenMP-MPI) environment as the number of iterations depends only on the number of blocks, irrespective of the number of processors.The only meaningful comparison between FSAI-DACG and BFSAI-DACG can be carried out in terms of iteration numbers which are smaller for BFSAI-DACG for a small number of processors. The gap between FSAI and BFSAI iterations reduces when the number of processors increases.


## 5.3. Comparison with the LOBPCG Eigensolver Provided byhypre

In order to validate the effectiveness of our preconditioning in the proposed DACG algorithm with respect to already available public parallel eigensolvers, the results given in Tables2 and 4 are compared with those obtained by the schemes implemented in the hypre software package [14]. The Locally Optimal Block Preconditioned Conjugate Gradient method (LOBPCG) [15] is experimented with, using the different preconditioners developed in the hypre project, that is, algebraic multigrid (AMG), diagonal scaling (DS), approximate inverse (ParaSails), additive Schwarz (Schwarz), and incomplete LU (Euclid). The hypre preconditioned CG is used for the inner iterations within LOBPCG. For details on the implementation of the LOBPCG algorithm, see, for instance, [22]. The selected preconditioner, ParaSails, is on its turn based on the FSAI preconditioner, so that the different FSAI-DACG and ParaSails-LOBPCG performances should be ascribed mainly to the different eigensolvers rather than to the preconditioners.We first carried out a preliminary set of runs with the aim of assessing the optimal value of the block sizebl parameter, that is, the size of the subspace where to seek for the eigenvectors. Obviously it must be bl ≥ s = 10. We fixed to 16 the number of processors and obtained the results summarized in Table 5 with different values of bl  ∈  [10,15]. We found that, only in problem CUBE-6091, a value of bl larger than 10, namely, bl = 12, yields an improvement in the CPU time. Note that we also made this comparison with different number of processors, and we obtained analogous results.Table 5
Iterations and CPU time for the iterative solver of LOBPCG-hypre preconditioned by Parasails with different values of bl and p=16 processors.
Matrixbl = 10bl = 11bl = 12bl = 15iterTiteriterTiteriterTiteriterTiterFAULT-63915679.515785.315796.1160128.1PO-87845117.041131.638151.335192.6GEO-143823123.772173.730152.5121291.1CUBE-60911011670.51432414.0381536.7351680.9Table6 presents the number of iterations and timings using the LOBPCG algorithm in the hypre package. The LOBPCG wall clock time is obtained with the preconditioner allowing for the best performance in the specific problem at hand, that is, ParaSails for all the problems. Using AMG as the preconditioner did not allow for convergence in three cases out of four, with the only exception of the FAULT-639 problem, in which the CPU timings were however very much larger than using ParaSails.Table 6
Number of iterations, timings, and scalability of LOBPCG-hypre preconditioned by Parasails.
piterTprecTiterTtotSp(4)Ep(4)FAULT-639 pcgitr =541552.5331.2333.781561.3167.6168.97.90.99161560.879.580.316.61.04321500.538.839.334.01.06641450.322.222.559.40.931281570.114.814.989.70.70PO-878 pcgitr = 304453.3438.4441.78501.3232.3234.07.60.9416451.0117.0118.015.00.9432450.763.263.927.60.8664470.434.434.850.80.79128410.319.4419.7489.50.70GEO-1438pcgitr = 304267.7478.0485.78224.0256.8260.87.50.9316232.1123.7125.815.40.9632281.273.174.326.20.8264230.835.536.353.50.84128250.520.320.893.20.73256260.312.913.2147.20.57piterTprecTiterTtotSp(16)Ep(16)CUBE-609116389.21536.71545.932364.7807.5812.230.50.9564383.2408.2411.460.10.94128411.6251.4253.097.80.76256350.9105.9106.8231.60.90512390.665.365.9375.30.731024370.337.738.0650.90.64All matrices have to be preliminarily scaled by their maximum coefficient in order to allow for convergence. To make the comparison meaningful, the outer iterations of the different methods are stopped when the average relative error measure of the computed leftmost eigenpairs gets smaller than10-10, in order to obtain a comparable accuracy as in the other codes. We also report in Table 6 the number of inner preconditioned CG iterations (pcgitr).To better compare our FSAI DACG with the LOBPCG method, we depict in Figure1 the total CPU time versus the number of processor for the two codes. FSAI-DACG and LOBPCG provide very similar scalability, being the latter code a little bit more performing on the average. On the FAULT-639 problem, DACG reveals faster than LOBPCG, irrespective of the number of processors employed.Comparison between FSAI-DACG and LOBPCG-hypre in terms of total CPU time for different number of processors.
(a)(b)(c)(d)Finally, we have carried out a comparison of the two eigensolvers in the computation of only the leftmost eigenpair. Differently from LOBPCG, which performs a simultaneous approximation of all the selected eigenpairs, DACG proceeds in the computation of the selected eigenpairs in a sequential way. For this reason, DACG should be the better choice, at least in principle, when just one eigenpair is sought. We investigate this feature, and the results are summarized in Table7. We include the total CPU time and iteration count needed by LOBPCG and FSAI-DACG to compute the leftmost eigenpair with 16 processors. For the LOBPCG code we report only the number of outer iterations.Table 7
Performance of LOBPCG-hypre with Parasails and 16 processors in the computation of the smallest eigenvalue using bl = 1 and bl = 2 and FSAI-DACG.
LOBPCG,bl = 1LOBPCG,bl = 2FSAI-DACGiterTtotiterTtotiterTtotFAULT-63914410.113217.5103015.4PO-8789943.23429.199320.4GEO-14935540.22637.375427.0CUBE-60911445218.858522.43257561.1The parameters used to construct the FSAI preconditioner for these experiments are as follows:(1)
FAULT-639.δ=0.1, d=2, ε=0.05,(2)
PO-878.δ=0.2, d=4, ε=0.1,(3)
GEO-1438.δ=0.1, d=2, ε=0.1,(4)
CUBE-6091.δ=0.0, d=1, ε=0.05.These parameters differ from those employed to compute the FSAI preconditioner in the assessment of the 10 leftmost eigenpairs and have been selected in order to produce a preconditioner relatively cheap to compute. This is so because otherwise the setup time would prevail over the iteration time. Similarly, to compute just one eigenpair with LOBPCG we need to setup a different value for pcgitr, the number of inner iterations. As it can be seen from Table7, in the majority of the test cases, LOBPCG takes less time to compute 2 eigenpairs than just only 1. FSAI-DACG reveals more efficient than the best LOBPCG on problems PO-878 and GEO-1438. On the remaining two problems the slow convergence exhibited by DACG is probably due to the small relative separation ξ1 between λ1 and λ2.


## 6. Conclusions

We have presented the parallel DACG algorithm for the partial eigensolution of large and sparse SPD matrices. The scalability of DACG, accelerated with FSAI-type preconditioners, has been studied on a set of test matrices of very large size arising from real engineering mechanical applications. Our FSAI-DACG code has shown comparable performances with the LOBPCG eigensolver within the well-known public domain package,hypre. Numerical results reveal that not only the scalability achieved by our code is roughly identical to that of hypre but also, in some instances, FSAI-DACG proves more efficient in terms of absolute CPU time. In particular, for the computation of the leftmost eigenpair, FSAI-DACG is more convenient in 2 problems out of 4.



---

*Source: 872901-2012-04-29.xml*