Parallel Rayleigh Quotient optimization with FSAI-based preconditioning

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 parallel hypre 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.


Introduction
The computation by iterative methods of the s leftmost partial eigenspectrum of the generalized eigenproblem: where A, B ∈ R n×n are large sparse symmetric positive definite (SPD) matrices, is an important and difficult task in many applications.It has become increasingly wide spread 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) [2,19], the Jacobi-Davidson (JD) algorithm [21] and schemes based on preconditioned Conjugate Gradient minimization of the Rayleigh quotient [3,16].
The basic idea of the later is to minimize the Rayleigh quotient: 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 [3,7] which has been shown to be competitive with the Jacobi Davidson method and with the PARPACK package [8].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 [18].
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 [4,6,5].The effectiveness of the FSAI preconditioner in the acceleration of DACG is compared to that of the Block FSAI-IC preconditioner, recently developed in [13], 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 [11].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 [1] which implements a number of preconditioners which can be used in combination with the Locally Optimal Block PCG (LOBPCG) iterative eigensolver [14].
The results presented in this paper shows 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 Section 2 we describe the DACG Algorithm; in Section 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×10 8 nonzeros.A comparison with the hypre eigensolver code is also included.Section 6 ends the paper with some conclusions.
Indicating with M the preconditioning matrix, i. e. M ≈ A −1 , the s leftmost eigenpairs are computed by the conjugate gradient procedure [7] described in Algorithm 1.
, with

end do
The 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 , u j ) it can be proved that the number of iterations is proportional to the square root of the condition number ξ j = κ(H j ) of the Hessian of the Rayleigh Quotient in the stationary point u j [3].It turns out that H j is similar to (A − λ j I)M which is not SPD.However, H j operates on the orthognal 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 where in the ideal case M ≡ A −1 , we have Therefore, even though A −1 is not the optimal preconditioner for A − λ j I, however, if M is a good preconditioner of A then the condition number κ(H j ) will approach ξ j .

The FSAI Preconditioner
The FSAI preconditioner, initially proposed in [17] and [18], has been later developed and implemented in parallel by Bergamaschi et al. in [4].Here, we only shortly recall the main features of this preconditioner.Given and SPD matrix A the FSAI preconditioner approximately factorize its inverse as a product of two sparse triangular matrices as The choice of nonzeros in W are based on a sparsity pattern which in our work may be the same as Ãd where Ã is the result of prefiltration [6] of A i.e. dropping of all elements below of a threshold parameter δ.The entries of W are computed by minimizing the Frobenius norm of I − W L 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.

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 W T ), 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 row i of matrix W of FSAI preconditioner is computed independently of each other, by solving a small SPD dense linear system of size n i 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 non local 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 non local 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 W T .
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 [20] showing its effectiveness up to 1024 processors.

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 that D is an arbitrary non-singular block diagonal matrix consisting of n b equal size blocks.
Let S L and S BD be a sparse lower triangular and a dense block diagonal non-zero pattern, respectively, for a n×n matrix.Even though not strictly necessary, for the sake of simplicity assume that S BD consists of n b diagonal blocks with equal size m = n/n b and let D ∈ R n×n be an arbitrary full-rank matrix with non-zero pattern S BD .
Consider the set of lower block triangular matrices F with a prescribed non-zero pattern S BL and minimize over F the Frobenius norm: where L is the exact lower Cholesky factor of an SPD matrix A. A matrix F satisfying the minimality condition (4) 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) with respect to the unknown entries [F ] ij , (i, j) ∈ S BL , yields the solution to n independent dense subsystems which 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 F AF T into n b diagonal blocks.However, as D is arbitrary, it is still not ensured that F AF T is better than A in an iterative method, so it is necessary to precondition F AF T again.As F AF T 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 It is easy to show that F is guaranteed to exist with SPD matrices and B i b is SPD, too [13].Using an IC decomposition with partial fill-in for each block B i b and collecting in J the lower IC factors, the resulting preconditioned matrix reads: with the final preconditioner: M in equation ( 6) is the BFSAI-IC preconditioner of A.
For its computation BFSAI-IC needs the selection of n b and S L .The basic requirement for the number of blocks n b 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 n b = p.Hence, n b is by default set equal to p.By distinction, the choice of S 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 non-zero pattern of A d 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 S L as the lower block triangular pattern of A d .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 post-filtration parameter that allows for dropping the smallest entries of F .
In particular, An OpenMP implementation of the algorithms above is available in [12].

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, • FAULT-639: is obtained from a structural problem discretizing a faulted gas reservoir with tetrahedral Finite Elements and triangular Interface Elements [10].
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.
• 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 [9]).
• 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 x 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 [22].
• 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 at http://www.cise.ufl.edu/research/sparse/matrices.The computational performance of FSAI is compared to the one obtained by using BFSAI as implemented in [13].The comparison is done evaluating the number of iterations n iter to converge at the same tolerance, the wall clock time in seconds T prec and T iter for the preconditioner computation and the eigensolver to converge, respectively, with the total time T tot = T prec + T iter .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 with T p the total CPU elapsed times expressed in seconds on p processors.We introduce a relative measure of the parallel efficiency achieved by the code, S (p) 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 E (p) p the corresponding efficiency:

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 Fig. 1) 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.

BFSAI-IC-DACG results
We present in this Section the results of DACG accelerated by the BFSAI-IC preconditioner for the approximation of the s = 10 leftmost eigenpairs of the matrices described above.
Table 3 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: The user-specified parameters for BFSAI-IC given above provide evidence that it is important to build a dense preconditioner based on the lower non-zero pattern of A 3 (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 for p ≤ 8.For this reason the number of iterations reported in Table 4 are obtained with increasing number of blocks n b and with p = 8 processors.This iteration numbers 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.

Comparison with the LOBPCG eigensolver provided by hypre
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 Tables 2 and 4 are compared with those obtained by the schemes implemented in the hypre software package [1].The Locally Optimal Block Preconditioned Conjugate Gradient method (LOBPCG) [14] is experimented with, using the different preconditioners developed in the hypre project, i.e. 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 [15].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 size bl parameter, i.e. 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 obtain analogous results.6 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, i.e.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.
All 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 than 10 −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 Figure 1 the total CPU time vs 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.
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 simulta-neous 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 Table 7.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.
The parameters used to construct the FSAI preconditioner for these experiments are as follows: 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 Table 7, 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 .

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.
a parameter that controls the fill-in of B i b and determines the maximum allowable number of nonzeros for each row of B i b 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 Li b denoting the maximum allowable number of nonzeros for each row of Li b in addition to the corresponding entries of B i b .

Figure 1 :
Figure 1: Comparison between FSAI-DACG and LOBPCG-hypre in terms of total CPU time for different number of processors.FAULT-639 PO-878

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.

Table 3 :
Performance of BFSAI-DACG for matrix PO-878 with 2 to 8 processors and different parameter values.

Table 4 :
Number of iterations for BFSAI-DACG in the computations of the 10 leftmost eigenpairs.

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.

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.