# Parallel RFSAI-BFGS Preconditioners for Large Symmetric Eigenproblems

Authors: L. Bergamaschi; A. Martínez
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

We propose a parallel preconditioner for the Newton method in the computation of the leftmost eigenpairs of large and sparse symmetric positive definite matrices. A sequence of preconditioners starting from an enhanced approximate inverse RFSAI (Bergamaschi and Martínez, 2012) and enriched by a BFGS-like update formula is proposed to accelerate the preconditioned conjugate gradient solution of
the linearized Newton system to solveAu=q(u)u, q(u) being the Rayleigh quotient. In a previous work (Bergamaschi and Martínez, 2013) the sequence of preconditioned Jacobians is proven to remain close to the identity matrix if the initial preconditioned Jacobian is so. Numerical results onto matrices arising from various realistic problems with size up to 1.5 million unknowns account for the efficiency and the scalability of the proposed low rank update of the RFSAI preconditioner. The overall RFSAI-BFGS preconditioned Newton algorithm has shown comparable efficiencies with a well-established eigenvalue solver on all the test problems.


---

## Body

## 1. Introduction

The computation of thep≪n leftmost eigenpairs of a symmetric positive definite (SPD) matrix A is a common task in many scientific applications. Typical examples are offered by the vibrational analysis of mechanical structures [1]: the spectral superposition approach for the solution of large sets of 1st order linear differential equations [2] and the approximation of the generalized inverse of the graph Laplacian [3], to mention a few.We denote(1)λ1<λ2<⋯<λp<⋯<λn
as the eigenvalues of A and
(2)v1,v2,…,vp,…,vn
as the corresponding (normalized) eigenvectors.In this paper, we propose to use the Newton method in the unit sphere [4, 5] or Newton-Grassman method, which constructs a sequence of vectors {uk} by solving the linear systems
(3)(I-ukuk⊤)(A-θkI)(I-ukuk⊤)s=-(Auk-θkuk),θk=uk⊤Aukuk⊤uk
ensuring that the correction s be orthogonal to uk. Then the next approximation is set as uk+1=t∥t∥-1, where t=uk+s. Linear system (3) is shown to be better conditioned than the one with A-θkI. For SPD matrices, this scheme is shown to converge cubically to the wanted eigenvector, provided that the system (3) is solved exactly. The same linear system represents the correction equation in the well-known Jacobi-Davidson method [6], which in its turn can be viewed as an accelerated inexact Newton method [7]. When A is SPD and the leftmost eigenpairs are being sought, it has been proved in [8] that the preconditioned conjugate gradient (PCG) method can be employed in the solution of the correction equation.Starting from the findings in [9–11] the main contribution of this paper is the development of a sequence of preconditioners {Pk} for the PCG solution of the Newton correction equation (3), based on the BFGS update of a given parallel initial preconditioner for the coefficient matrix A. A similar approach has been used in [12] where a rank-two modification of a given preconditioner is used to accelerated MINRES in the framework of the inexact Rayleigh quotient iteration. Updating (cheaply) a given preconditioner to obtain a new preconditioner for a slightly changed linear system is also common in optimization problems. See for example, [13–15].The initial Newton vector is obtained after a small number of iterations of a conjugate gradient procedure for the minimization of the Rayleigh quotient (DACG, [16]). As the initial preconditioner we choose RFSAI [17], which is built by recursively applying the (factorized sparse approximate inverse) FSAI  preconditioner developed in [18]. We elected a factorized approximate inverse preconditioner (AIP) since it is more naturally parallelizable than preconditioners based on ILU factorizations. Moreover, factorized variants of AIP provide better approximations to A-1 for the same amount of storage than nonfactorized ones, because they can express denser matrices than the total number of nonzeros in their factors [19]. The RFSAI formulation enhances this characteristic, since the final preconditioner is implicitly built by a product of four or more triangular factors.The combined DACG-Newton algorithm is used in the approximation ofp=10 eigenpairs of a number of matrices arising from various realistic applications of size up to 1.5×106 and number of nonzeros up to 6×107. Numerical results show that, in the solution of the correction equation, the PCG method preconditioned by BFGS displays much faster convergence than the same method when the preconditioner is kept fixed during the Newton process, in every test case. Moreover, the proposed approach is shown to outperform the pure DACG method.


## 2. BFGS Sequence of Preconditioners

### 2.1.  Choice of the Initial Preconditioner

Following the developments in [17], we propose an implicit enlargement of the sparsity pattern using a banded target matrix B; the lower factor of the FSAI preconditioner is obtained by minimizing ∥B-GL∥F over the set of matrices G having a fixed sparsity pattern. Denoting with Gout the result of this minimization, we compute explicitly the preconditioned matrix S=GoutAGoutT and then evaluate a second FSAI factor Gin for S. Thus, the final preconditioner can be written as the product GoutTGinTGinGout. This RFSAI—recursive FSAI—procedure can be iterated a number of times to yield a preconditioner written as a product of several triangular factors. This preconditioner has shown its efficiency in the acceleration of the PCG onto realistic geomechanical problems of very large size [17].We denote FSAIout as the procedure which computes theGout factor by minimizing ∥B-GL∥, where B is an arbitrary banded matrix. FSAIout depends on the nband parameter in addition to the usual FSAI parameters: δout prefiltration threshold, dout power of A~ defining the sparsity pattern, and εout postfiltration parameter.The second preconditioner factor,Gin, is the result of the FSAI procedure applied to the whole product matrix GoutAGoutT, with parameters δin,din, and εin. The steps to obtain the final RFSAI preconditioner are summarized in Algorithm 1.Algorithm 1:RFSAI computation.
INPUT:nband, δout,dout,εout,δin,din,εin
Compute the first lower triangular factor:Gout=  FSAIout(A,nband,δout,dout,εout)
Compute the product:A(1)=GoutAGoutT
Compute the second lower triangular factor:Gin=FSAI(A(1),δin,din,εin)Following the idea described in [10], we propose a sequence of preconditioners for the Newton systems using the BFGS rank-two update. To precondition the initial Newton system as follows:
(4)J0s0=-r,
where
(5)J0=(I-u0u0⊤)(A-θ0I)(I-u0u0⊤),r=-(Au0-θ0u0),θ0=u0⊤Au0.
We choose to use a projected RFSAI preconditioner, P0=(I-u0u0⊤)P^0(I-u0u0⊤) with P^0=WTW and W=GinGout.


### 2.2. Update of Initial Preconditioner by BFGS-Like Rank-Two Corrections

A sequence of projected preconditioners for the subsequent linear systemsJk+1sk+1=-rk+1 may be defined by
(6)P^k+1=-ss⊤s⊤r+(I-sr⊤s⊤r)P^k(I-rs⊤s⊤r),
where s≡sk the solution of the previous correction equation and r≡rk. In view of the cubic convergence of the Newton process, we used the residual -rk instead of yk=rk+1-rk.Theorem3 of Section 3 will state that the preconditioner defined in (6) is SPD if P^k is so.


## 2.1.  Choice of the Initial Preconditioner

Following the developments in [17], we propose an implicit enlargement of the sparsity pattern using a banded target matrix B; the lower factor of the FSAI preconditioner is obtained by minimizing ∥B-GL∥F over the set of matrices G having a fixed sparsity pattern. Denoting with Gout the result of this minimization, we compute explicitly the preconditioned matrix S=GoutAGoutT and then evaluate a second FSAI factor Gin for S. Thus, the final preconditioner can be written as the product GoutTGinTGinGout. This RFSAI—recursive FSAI—procedure can be iterated a number of times to yield a preconditioner written as a product of several triangular factors. This preconditioner has shown its efficiency in the acceleration of the PCG onto realistic geomechanical problems of very large size [17].We denote FSAIout as the procedure which computes theGout factor by minimizing ∥B-GL∥, where B is an arbitrary banded matrix. FSAIout depends on the nband parameter in addition to the usual FSAI parameters: δout prefiltration threshold, dout power of A~ defining the sparsity pattern, and εout postfiltration parameter.The second preconditioner factor,Gin, is the result of the FSAI procedure applied to the whole product matrix GoutAGoutT, with parameters δin,din, and εin. The steps to obtain the final RFSAI preconditioner are summarized in Algorithm 1.Algorithm 1:RFSAI computation.
INPUT:nband, δout,dout,εout,δin,din,εin
Compute the first lower triangular factor:Gout=  FSAIout(A,nband,δout,dout,εout)
Compute the product:A(1)=GoutAGoutT
Compute the second lower triangular factor:Gin=FSAI(A(1),δin,din,εin)Following the idea described in [10], we propose a sequence of preconditioners for the Newton systems using the BFGS rank-two update. To precondition the initial Newton system as follows:
(4)J0s0=-r,
where
(5)J0=(I-u0u0⊤)(A-θ0I)(I-u0u0⊤),r=-(Au0-θ0u0),θ0=u0⊤Au0.
We choose to use a projected RFSAI preconditioner, P0=(I-u0u0⊤)P^0(I-u0u0⊤) with P^0=WTW and W=GinGout.


## 2.2. Update of Initial Preconditioner by BFGS-Like Rank-Two Corrections

A sequence of projected preconditioners for the subsequent linear systemsJk+1sk+1=-rk+1 may be defined by
(6)P^k+1=-ss⊤s⊤r+(I-sr⊤s⊤r)P^k(I-rs⊤s⊤r),
where s≡sk the solution of the previous correction equation and r≡rk. In view of the cubic convergence of the Newton process, we used the residual -rk instead of yk=rk+1-rk.Theorem3 of Section 3 will state that the preconditioner defined in (6) is SPD if P^k is so.


## 3. Theoretical Analysis of the Preconditioner

### 3.1. Finding the Smallest Eigenpair

We recall in this section the main theoretical results regarding the sequence of preconditioners previously defined. Differently from the classical papers which study BFGS convergence properties, here our Jacobian matrixJ(u)=(I-uu⊤)(A-q(u)I)(I-uu⊤) is singular whatever u, in particular it is singular when u is equal to the exact eigenvector. The theoretical analysis of the “goodness” of the preconditioner will be, therefore, completely different, though obtaining similar results, than that proposed, for example, in [10, 20]. In the following developments we will indicate as v1 the exact eigenvector corresponding to the smallest exact eigenvalue λ1. The error vector at step k is denoted by ek=uk-v1, while the error in the eigenvalue approximation is εk=θk-λ1(>0). It is easily proved that there is a constant M independent of k such that
(7)εk≤M∥ek∥2.Remark 1.
At first sight the Jacobian matrix in the correction equation is singular, but this does not matter since the PCG algorithm is run within the subspace of vectors orthogonal touk (in fact also r⊤uk=0). Thus, notion of positive definiteness, eigenvalue distribution, condition number, norms, and so forth, apply as usual but with respect to matrices restricted to this subspace.The following lemma will bound the extremal eigenvalues ofJk in the subspace orthogonal to uk.Lemma 2.
There is a positive numberδ such that if ∥ek∥<δ then
(8)Jk=(I-ukuk⊤)(A-θkI)(I-ukuk⊤)
is SPD in the subspace orthogonal to uk. Moreover the following bounds hold:
(9)λ2-λ12<z⊤Jkz<λn
for every unit norm vector z orthogonal to uk.Proof.
See [21].The previous lemma allows us to prove that the preconditioner defined in (6) is SPD, as stated in the following theorem.Theorem 3.
If the correction equation is solved exactly, then any matrixP^k defined by (6) is SPD and hence Pk is SPD in the subspace orthogonal to uk.Proof.
See [21].Let us define the difference between the preconditioned Jacobian and the identity matrix as(10)Ek=I-Jk1/2PkJk1/2.
Since by definition we have Jkuk=0 then uk is the eigenvector of Jk corresponding to the zero eigenvalue. Hence, since also Jk1/2uk=0 the error Ek can also be defined as
(11)Ek=I-Jk1/2PkJk1/2=I-Jk1/2(I-ukuk⊤)×P^k(I-ukuk⊤)Jk1/2=I-Jk1/2P^kJk1/2.
The following technical lemma will bound the norm of P^k in terms of that of Ek. Being P^k SPD we can define its norm in the space orthogonal to uk as
(12)∥P^k∥=supw⊥uk,w≠0w⊤P^kww⊤w.Lemma 4.
There is a positive numberδ such that if ∥ek∥<δ then
(13)∥P^k∥≤2λ2-λ1(1+∥Ek∥).Proof.
See [21].The next lemma will relate the norms of the differences and of the error vector ek.Lemma 5.
There exists a positive numberδ s.t. if ∥ek∥<δ then
(14)∥s∥≤3∥ek∥.Proof.
See [21].Before stating Theorem7 we need to state as a last preliminary result that also the difference between the square root of two consecutive Jacobians is bounded in terms of the norm of the error vector.Lemma 6.
LetSk=Jk+11/2-Jk1/2. Then there is a positive number δ s.t. if ∥ek∥<δ then
(15)∥Sk∥≤c3∥ek∥.
for a suitable constant c3.Proof.
See [21].The following theorem will state the main theoretical result of this Section: the so calledbounded deterioration [22] of the preconditioner at step k+1 with respect to that of step k. In other words it can be proved that the distance of the preconditioned matrix from the identity matrix at step k+1 is less or equal than that at step k plus a constant that may be small as desired, depending on the closeness of u0 to the exact eigenvector. We report also the proof of this theorem, which is taken from [21].Theorem 7.
Letδ0 be such that ∥E0∥<δ0, there is a positive number δ s.t. if ∥e0∥<δ then
(16)∥Ek+1∥≤∥Ek∥+K∥ek∥
for a suitable constant K.Proof.
After definingN=SkP^k+1Jk+11/2+Jk+11/2P^k+1Sk+SkP^k+1Sk, we can write
(17)Ek+1=I-Jk+11/2P^k+1Jk+11/2=I-(Jk1/2+Sk)P^k+1(Jk1/2+Sk)=I-Jk1/2P^k+1Jk1/2-N=I-Jk1/2ss⊤s⊤JksJk1/2-Jk1/2(I-ss⊤Jks⊤Jks)P^k(I-Jkss⊤s⊤Jks)×Jk1/2-N=I-Jk1/2ss⊤Jk1/2s⊤Jks-(I-Jk1/2ss⊤Jk1/2s⊤Jks)Jk1/2P^kJk1/2×(I-Jk1/2ss⊤Jk1/2s⊤Jks)-N=W+W(Ek-I)W+N=WEkW-N,
where we set w=Jk1/2s/∥Jk1/2s∥ and W=I-ww⊤; W is an orthogonal projector since ∥w∥=1. To bound ∥N∥, we will use Lemmas 4 and 6 as follows:
(18)∥N∥≤2λ2-λ1(1+∥Ek+1∥)×(2c3∥ek∥λn+c32δ∥ek∥)=c4(1+∥Ek+1∥)∥ek∥.
Now taking norms in (17) yields
(19)∥Ek+1∥≤∥Ek∥+c4(1+∥Ek+1∥)∥ek∥,
which can be rewritten as
(20)∥Ek+1∥(1-c4∥ek∥)≤∥Ek∥+c4∥ek∥.
From (20), we derive a bound for ∥Ek∥. If δ<1/2c4 then
(21)∥Ek∥≤2∥Ek-1∥+1≤⋯≤2k∥E0∥+2k-1≤2k(δ0+1)=c5.
Again from (20) and using the bound (21) we finally have
(22)∥Ek-1∥≤∥Ek∥+c4∥ek∥1-c4∥ek∥≤(∥Ek∥+c4∥ek∥)·(1+2c4∥ek∥)≤∥Ek∥+2c4∥ek∥c5+c4(1+2c4δ)∥ek∥=∥Ek∥+c4(2c5+1+2c4δ)∥ek∥.
Setting K=c4(2c5+1+2c4δ) completes the proof.


### 3.2.  Computing Several Eigenpairs

When seeking an eigenvalue different fromλ1, say λj, the Jacobian matrix changes as
(23)Jk=(I-QQ⊤)(A-θkI)(I-QQ⊤),
where Q=[v1v2⋯vjuk] is the matrix whose first j columns are the previously computed eigenvectors. Analogously, also the preconditioner must be chosen orthogonal to Q as
(24)Pk+1=(I-QQ⊤)P^k+1(I-QQ⊤).
The theoretical analysis developed in Section 3.1 applies with small technical variants also in this case since it is readily proved that Jk+11/2Pk+1Jk+11/2=Jk+11/2P^k+1Jk+11/2. The most significant changes regard the definition of εk=θk-λj,  ek=uk-vj, and the statement of Lemma 2 (and the proof of Lemma 4 that uses its results); namely, the bound for the smallest eigenvalue of Jk which in the general case becomes
(25)z⊤Jkz>λj+1-λj2
for every unit norm vector zsuch that Q⊤z=0.


## 3.1. Finding the Smallest Eigenpair

We recall in this section the main theoretical results regarding the sequence of preconditioners previously defined. Differently from the classical papers which study BFGS convergence properties, here our Jacobian matrixJ(u)=(I-uu⊤)(A-q(u)I)(I-uu⊤) is singular whatever u, in particular it is singular when u is equal to the exact eigenvector. The theoretical analysis of the “goodness” of the preconditioner will be, therefore, completely different, though obtaining similar results, than that proposed, for example, in [10, 20]. In the following developments we will indicate as v1 the exact eigenvector corresponding to the smallest exact eigenvalue λ1. The error vector at step k is denoted by ek=uk-v1, while the error in the eigenvalue approximation is εk=θk-λ1(>0). It is easily proved that there is a constant M independent of k such that
(7)εk≤M∥ek∥2.Remark 1.
At first sight the Jacobian matrix in the correction equation is singular, but this does not matter since the PCG algorithm is run within the subspace of vectors orthogonal touk (in fact also r⊤uk=0). Thus, notion of positive definiteness, eigenvalue distribution, condition number, norms, and so forth, apply as usual but with respect to matrices restricted to this subspace.The following lemma will bound the extremal eigenvalues ofJk in the subspace orthogonal to uk.Lemma 2.
There is a positive numberδ such that if ∥ek∥<δ then
(8)Jk=(I-ukuk⊤)(A-θkI)(I-ukuk⊤)
is SPD in the subspace orthogonal to uk. Moreover the following bounds hold:
(9)λ2-λ12<z⊤Jkz<λn
for every unit norm vector z orthogonal to uk.Proof.
See [21].The previous lemma allows us to prove that the preconditioner defined in (6) is SPD, as stated in the following theorem.Theorem 3.
If the correction equation is solved exactly, then any matrixP^k defined by (6) is SPD and hence Pk is SPD in the subspace orthogonal to uk.Proof.
See [21].Let us define the difference between the preconditioned Jacobian and the identity matrix as(10)Ek=I-Jk1/2PkJk1/2.
Since by definition we have Jkuk=0 then uk is the eigenvector of Jk corresponding to the zero eigenvalue. Hence, since also Jk1/2uk=0 the error Ek can also be defined as
(11)Ek=I-Jk1/2PkJk1/2=I-Jk1/2(I-ukuk⊤)×P^k(I-ukuk⊤)Jk1/2=I-Jk1/2P^kJk1/2.
The following technical lemma will bound the norm of P^k in terms of that of Ek. Being P^k SPD we can define its norm in the space orthogonal to uk as
(12)∥P^k∥=supw⊥uk,w≠0w⊤P^kww⊤w.Lemma 4.
There is a positive numberδ such that if ∥ek∥<δ then
(13)∥P^k∥≤2λ2-λ1(1+∥Ek∥).Proof.
See [21].The next lemma will relate the norms of the differences and of the error vector ek.Lemma 5.
There exists a positive numberδ s.t. if ∥ek∥<δ then
(14)∥s∥≤3∥ek∥.Proof.
See [21].Before stating Theorem7 we need to state as a last preliminary result that also the difference between the square root of two consecutive Jacobians is bounded in terms of the norm of the error vector.Lemma 6.
LetSk=Jk+11/2-Jk1/2. Then there is a positive number δ s.t. if ∥ek∥<δ then
(15)∥Sk∥≤c3∥ek∥.
for a suitable constant c3.Proof.
See [21].The following theorem will state the main theoretical result of this Section: the so calledbounded deterioration [22] of the preconditioner at step k+1 with respect to that of step k. In other words it can be proved that the distance of the preconditioned matrix from the identity matrix at step k+1 is less or equal than that at step k plus a constant that may be small as desired, depending on the closeness of u0 to the exact eigenvector. We report also the proof of this theorem, which is taken from [21].Theorem 7.
Letδ0 be such that ∥E0∥<δ0, there is a positive number δ s.t. if ∥e0∥<δ then
(16)∥Ek+1∥≤∥Ek∥+K∥ek∥
for a suitable constant K.Proof.
After definingN=SkP^k+1Jk+11/2+Jk+11/2P^k+1Sk+SkP^k+1Sk, we can write
(17)Ek+1=I-Jk+11/2P^k+1Jk+11/2=I-(Jk1/2+Sk)P^k+1(Jk1/2+Sk)=I-Jk1/2P^k+1Jk1/2-N=I-Jk1/2ss⊤s⊤JksJk1/2-Jk1/2(I-ss⊤Jks⊤Jks)P^k(I-Jkss⊤s⊤Jks)×Jk1/2-N=I-Jk1/2ss⊤Jk1/2s⊤Jks-(I-Jk1/2ss⊤Jk1/2s⊤Jks)Jk1/2P^kJk1/2×(I-Jk1/2ss⊤Jk1/2s⊤Jks)-N=W+W(Ek-I)W+N=WEkW-N,
where we set w=Jk1/2s/∥Jk1/2s∥ and W=I-ww⊤; W is an orthogonal projector since ∥w∥=1. To bound ∥N∥, we will use Lemmas 4 and 6 as follows:
(18)∥N∥≤2λ2-λ1(1+∥Ek+1∥)×(2c3∥ek∥λn+c32δ∥ek∥)=c4(1+∥Ek+1∥)∥ek∥.
Now taking norms in (17) yields
(19)∥Ek+1∥≤∥Ek∥+c4(1+∥Ek+1∥)∥ek∥,
which can be rewritten as
(20)∥Ek+1∥(1-c4∥ek∥)≤∥Ek∥+c4∥ek∥.
From (20), we derive a bound for ∥Ek∥. If δ<1/2c4 then
(21)∥Ek∥≤2∥Ek-1∥+1≤⋯≤2k∥E0∥+2k-1≤2k(δ0+1)=c5.
Again from (20) and using the bound (21) we finally have
(22)∥Ek-1∥≤∥Ek∥+c4∥ek∥1-c4∥ek∥≤(∥Ek∥+c4∥ek∥)·(1+2c4∥ek∥)≤∥Ek∥+2c4∥ek∥c5+c4(1+2c4δ)∥ek∥=∥Ek∥+c4(2c5+1+2c4δ)∥ek∥.
Setting K=c4(2c5+1+2c4δ) completes the proof.


## 3.2.  Computing Several Eigenpairs

When seeking an eigenvalue different fromλ1, say λj, the Jacobian matrix changes as
(23)Jk=(I-QQ⊤)(A-θkI)(I-QQ⊤),
where Q=[v1v2⋯vjuk] is the matrix whose first j columns are the previously computed eigenvectors. Analogously, also the preconditioner must be chosen orthogonal to Q as
(24)Pk+1=(I-QQ⊤)P^k+1(I-QQ⊤).
The theoretical analysis developed in Section 3.1 applies with small technical variants also in this case since it is readily proved that Jk+11/2Pk+1Jk+11/2=Jk+11/2P^k+1Jk+11/2. The most significant changes regard the definition of εk=θk-λj,  ek=uk-vj, and the statement of Lemma 2 (and the proof of Lemma 4 that uses its results); namely, the bound for the smallest eigenvalue of Jk which in the general case becomes
(25)z⊤Jkz>λj+1-λj2
for every unit norm vector zsuch that Q⊤z=0.


## 4. Implementation

### 4.1. Choosing an Initial Eigenvector Guess

As mentioned in Section1, another important issue in the efficiency of the Newton approach for eigenvalue computation is represented by the appropriate choice of the initial guess. We propose here to perform some preliminary iterations using the DACG eigensolver [16, 23, 24], which is based on the preconditioned conjugate gradient (nonlinear) minimization of the Rayleigh quotient. This method has proven very robust, and not particularly sensitive to the initial vector, in the computation of a few eigenpairs of large SPD matrices. Recently in [25], the DACG method has been compared with the Jacobi-Davidson method in parallel environments, showing similar performances when seeking a small number of eigenpairs (1÷5).


### 4.2. Implementation of the BFGS Preconditioner Update

At a certain nonlinear iteration level,k, we need to compute c=Pkgl, where gl is the residual of the linear system at iteration l. Let us suppose we compute an initial preconditioner P0. Then, at the initial nonlinear iteration k=0, we simply have c=P0zl. At step k+1 the preconditioner P^k+1 is defined recursively by (6) while Pk+1 using (24) can be written as
(26)Pk+1=(I-QQ⊤)P^k+1(I-QQ⊤)=(I-QQ⊤){(I-sr⊤s⊤r)P^k=(I-QQ⊤)×(I-rs⊤s⊤r)-ss⊤s⊤r}×(I-QQ⊤).
To compute vector c first we observe that gl is orthogonal to Q so there is no need to apply matrix I-QQ⊤ on the right of (26). Application of preconditioner P^k+1 to the vector gl can be performed at the price of 2k dot products and 2k daxpys as described in Algorithm 2. The scalar products αk=sk⊤rk, which appear at the denominator of P^k+1, can be computed once and for all before starting the solution of the (k+1)th linear system. Last, the obtained vector c must be orthogonalized against the columns of Q by a classical Gram-Schimdt procedure.Algorithm 2:Computation of c=Pkgl for the BFGS preconditioner.
INPUT: Vectorgl, scalar products αs=ss⊤rs,s=0,…,k-1.
w
=
g
l
fors:=k-1   to  0
(1) as:=ss⊤w/αs
(2) w:=w-asrs
end for
c
=
P
^
0
w
fors:=0   to  k-1
(1) b:=rs⊤c/αs
(2) c:=c-(as+b)ss
end for
z
:
=
Q
⊤
c
c
:
=
c
-
Q
z


### 4.3. PCG Solution of the Correction Equation

As a Krylov subspace solver for the correction equation, we chose the preconditioned conjugate gradient (PCG) method since the JacobianJk has been shown to be SPD in the subspace orthogonal to uk. Regarding the implementation of PCG, we mainly refer to the work [8], where the author shows that it is possible to solve the linear system in the subspace orthogonal to uk and hence, the projection step needed in the application of Jk can be skipped. Moreover, we adopted the exit strategy for the linear system solution described in the previous paper, which allows for stopping the PCG iteration, in addition to the classical exit test based on a tolerance on the relative residual and on the maximum number of iterations, whenever the current solution xl satisfies
(27)∥rk,l∥=∥Axl-xl⊤Axlxl⊤xlxl∥<τ(xl⊤Axl),
or when the decrease of ∥rk,l∥ is slower than the decrease of ∥gl∥, because in this case further iterating does not improve the accuracy of the eigenvector. Note that this dynamic exit strategy implicitly defines an inexact Newton method since the correction equation is not solved “exactly”, that is, up to machine precision.We have implemented the PCG method as described in Algorithm 5.1 of [8] with the obvious difference in the application of the preconditioner which in this case is accomplished as described in Algorithm 2.


### 4.4. Implementation of the DACG-Newton Method

The BFGS preconditioner defined in Algorithm2 suffers from two main drawbacks, namely, increasing costs of memory for storing sand r and the increasing cost of preconditioner application with the iteration index k. Note that these drawbacks are common to many iterative schemes, such as sparse (Limited Memory) Broyden implementations [26], GMRES [27], and Arnoldi method for eigenvalue problems [28]. If the number of nonlinear iterations is high the application of BFGS update may be too costly as compared with the expected reduction in the iteration number. To this aim we define kmax the maximum number of rank-two corrections we allow. When the nonlinear iteration counter k is larger than kmax, the vectors si, ri, i=k-kmax are substituted with the last computed sk,  rk. Vectors {si,ri} are stored in a matrix V with n rows and 2×kmax columns.The implementation of our DACG-Newton method for computing the leftmost eigenpairs of large SPD matrices is described in Algorithm3.Algorithm 3:DACG-Newton algorithm.
INPUT: MatrixA;
number of sought eigenpairsneig;
tolerance and maximum number of its for the outer iteration:τ,  ITMAX;
tolerance for the initial eigenvector guessτDACG;
tolerance and maximum number of its for the inner iteration:τPCG.
ITMAXPCG;
parameters for the RFSAI preconditioner,
δout,dout,εout,  nband, for the 1st FSAI factor and,
δin,din,εin for the 2nd factor;
maximum allowed rank-two update in the BFGS preconditioner:kmax.
Q
~
:
=
[
].
ComputeP^0, an RFSAI preconditioner for A,
forj:=1   to  neig
(1) Choose x0 such that Q~⊤x0=0;
(2) Compute u0, an approximation to vj by the DACG procedure with initial
vectorx0, preconditioner P^0 and tolerance τDACG;
(3) k:=0, θk:=uk⊤Auk;
(4) while  ∥Auk-θkuk∥>τθk   and  k<IMAX   do
(1)Q:=[Q~uk].
(2) SolveJksk=-rk for sk⊥Q by the PCG method with preconditioner Pk and tolerance τPCG.
(3)uk+1:=(uk+sk)/∥uk+sk∥,  θk+1=uk+1⊤Auk+1.
(4)k1=k  MOD kmax; V(*,2k1+1):=sk,V(*,2k1+2):=rk.
(5)k:=k+1
(6) end while
(7) Assume vj=uk and λj=θk. Set Q~:=[Q~vj]
end for


### 4.5. Parallel Implementation

We have developed a parallel code which implements the construction and application inside parallel PCG, of the RFSAI algorithm along with the BFGS updates. The resulting program is written in Fortran 90 and exploits the MPI library for exchanging data among the processors. We use a block row distribution of all matrices, that is, with complete rows assigned to different processors. All the matrices are stored in static data structures in CSR format.Parallelization of the FSAI preconditioner, which is the basis of the parallel RFSAI construction, has been performed and tested for example, in [29–31] where prefiltration and postfiltration have been implemented together with a priori sparsity pattern based on nonzeros of Ad with d≤4.The code makes also use of an optimized parallel matrix-vector product which has been developed in [32] showing its effectiveness up to 1024 processors.


## 4.1. Choosing an Initial Eigenvector Guess

As mentioned in Section1, another important issue in the efficiency of the Newton approach for eigenvalue computation is represented by the appropriate choice of the initial guess. We propose here to perform some preliminary iterations using the DACG eigensolver [16, 23, 24], which is based on the preconditioned conjugate gradient (nonlinear) minimization of the Rayleigh quotient. This method has proven very robust, and not particularly sensitive to the initial vector, in the computation of a few eigenpairs of large SPD matrices. Recently in [25], the DACG method has been compared with the Jacobi-Davidson method in parallel environments, showing similar performances when seeking a small number of eigenpairs (1÷5).


## 4.2. Implementation of the BFGS Preconditioner Update

At a certain nonlinear iteration level,k, we need to compute c=Pkgl, where gl is the residual of the linear system at iteration l. Let us suppose we compute an initial preconditioner P0. Then, at the initial nonlinear iteration k=0, we simply have c=P0zl. At step k+1 the preconditioner P^k+1 is defined recursively by (6) while Pk+1 using (24) can be written as
(26)Pk+1=(I-QQ⊤)P^k+1(I-QQ⊤)=(I-QQ⊤){(I-sr⊤s⊤r)P^k=(I-QQ⊤)×(I-rs⊤s⊤r)-ss⊤s⊤r}×(I-QQ⊤).
To compute vector c first we observe that gl is orthogonal to Q so there is no need to apply matrix I-QQ⊤ on the right of (26). Application of preconditioner P^k+1 to the vector gl can be performed at the price of 2k dot products and 2k daxpys as described in Algorithm 2. The scalar products αk=sk⊤rk, which appear at the denominator of P^k+1, can be computed once and for all before starting the solution of the (k+1)th linear system. Last, the obtained vector c must be orthogonalized against the columns of Q by a classical Gram-Schimdt procedure.Algorithm 2:Computation of c=Pkgl for the BFGS preconditioner.
INPUT: Vectorgl, scalar products αs=ss⊤rs,s=0,…,k-1.
w
=
g
l
fors:=k-1   to  0
(1) as:=ss⊤w/αs
(2) w:=w-asrs
end for
c
=
P
^
0
w
fors:=0   to  k-1
(1) b:=rs⊤c/αs
(2) c:=c-(as+b)ss
end for
z
:
=
Q
⊤
c
c
:
=
c
-
Q
z


## 4.3. PCG Solution of the Correction Equation

As a Krylov subspace solver for the correction equation, we chose the preconditioned conjugate gradient (PCG) method since the JacobianJk has been shown to be SPD in the subspace orthogonal to uk. Regarding the implementation of PCG, we mainly refer to the work [8], where the author shows that it is possible to solve the linear system in the subspace orthogonal to uk and hence, the projection step needed in the application of Jk can be skipped. Moreover, we adopted the exit strategy for the linear system solution described in the previous paper, which allows for stopping the PCG iteration, in addition to the classical exit test based on a tolerance on the relative residual and on the maximum number of iterations, whenever the current solution xl satisfies
(27)∥rk,l∥=∥Axl-xl⊤Axlxl⊤xlxl∥<τ(xl⊤Axl),
or when the decrease of ∥rk,l∥ is slower than the decrease of ∥gl∥, because in this case further iterating does not improve the accuracy of the eigenvector. Note that this dynamic exit strategy implicitly defines an inexact Newton method since the correction equation is not solved “exactly”, that is, up to machine precision.We have implemented the PCG method as described in Algorithm 5.1 of [8] with the obvious difference in the application of the preconditioner which in this case is accomplished as described in Algorithm 2.


## 4.4. Implementation of the DACG-Newton Method

The BFGS preconditioner defined in Algorithm2 suffers from two main drawbacks, namely, increasing costs of memory for storing sand r and the increasing cost of preconditioner application with the iteration index k. Note that these drawbacks are common to many iterative schemes, such as sparse (Limited Memory) Broyden implementations [26], GMRES [27], and Arnoldi method for eigenvalue problems [28]. If the number of nonlinear iterations is high the application of BFGS update may be too costly as compared with the expected reduction in the iteration number. To this aim we define kmax the maximum number of rank-two corrections we allow. When the nonlinear iteration counter k is larger than kmax, the vectors si, ri, i=k-kmax are substituted with the last computed sk,  rk. Vectors {si,ri} are stored in a matrix V with n rows and 2×kmax columns.The implementation of our DACG-Newton method for computing the leftmost eigenpairs of large SPD matrices is described in Algorithm3.Algorithm 3:DACG-Newton algorithm.
INPUT: MatrixA;
number of sought eigenpairsneig;
tolerance and maximum number of its for the outer iteration:τ,  ITMAX;
tolerance for the initial eigenvector guessτDACG;
tolerance and maximum number of its for the inner iteration:τPCG.
ITMAXPCG;
parameters for the RFSAI preconditioner,
δout,dout,εout,  nband, for the 1st FSAI factor and,
δin,din,εin for the 2nd factor;
maximum allowed rank-two update in the BFGS preconditioner:kmax.
Q
~
:
=
[
].
ComputeP^0, an RFSAI preconditioner for A,
forj:=1   to  neig
(1) Choose x0 such that Q~⊤x0=0;
(2) Compute u0, an approximation to vj by the DACG procedure with initial
vectorx0, preconditioner P^0 and tolerance τDACG;
(3) k:=0, θk:=uk⊤Auk;
(4) while  ∥Auk-θkuk∥>τθk   and  k<IMAX   do
(1)Q:=[Q~uk].
(2) SolveJksk=-rk for sk⊥Q by the PCG method with preconditioner Pk and tolerance τPCG.
(3)uk+1:=(uk+sk)/∥uk+sk∥,  θk+1=uk+1⊤Auk+1.
(4)k1=k  MOD kmax; V(*,2k1+1):=sk,V(*,2k1+2):=rk.
(5)k:=k+1
(6) end while
(7) Assume vj=uk and λj=θk. Set Q~:=[Q~vj]
end for


## 4.5. Parallel Implementation

We have developed a parallel code which implements the construction and application inside parallel PCG, of the RFSAI algorithm along with the BFGS updates. The resulting program is written in Fortran 90 and exploits the MPI library for exchanging data among the processors. We use a block row distribution of all matrices, that is, with complete rows assigned to different processors. All the matrices are stored in static data structures in CSR format.Parallelization of the FSAI preconditioner, which is the basis of the parallel RFSAI construction, has been performed and tested for example, in [29–31] where prefiltration and postfiltration have been implemented together with a priori sparsity pattern based on nonzeros of Ad with d≤4.The code makes also use of an optimized parallel matrix-vector product which has been developed in [32] showing its effectiveness up to 1024 processors.


## 5. Numerical Results

### 5.1. Test Problems

We report the results of our experiments with the RFSAI preconditioner in the solution of a set of problems of large size.The test matrices are realistic examples arising from 3D FE discretization of fluid flow and geomechanical models. Described in detail as follows:(1)
FLOW3D-663 arises from tetrahedral FE discretization of an elliptic PDE describing fluid flow in porous media.(2)
EMILIA-923 arises from the regional geomechanical model of a deep hydrocarbon reservoir. It is obtained discretizing the structural problem with tetrahedral finite elements. Due to the complex geometry of the geological formation, it was not possible to obtain a computational grid characterized by regularly shaped elements.(3)
LONG-1103 is the structural SPD block obtained from a 3D coupled consolidation problem of a geological formation, discretized with tetrahedral finite elements on a higly irregular computational grid.(4)
GEO-1438 arises from a regional geomechanical model of the sedimentary basin underlying the Venice lagoon discretized by a linear FE with randomly heterogeneous properties. 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.Matrices 2 to 4 are publicly available in the University of Florida Sparse Matrix Collection athttp://www.cise.ufl.edu/research/sparse/matrices/. The size and number of nonzero elements for each matrix is provided in Table 1.Table 1
Sizen and number of nonzeros (nnz)  of the test matrices.
Name
n
nnz
F
LOW3D-663
663 390
9 762 648
EMILIA-923
923 136
41 005 206
LONG-1103
1 102 614
48 987 558
GEO-1438
1 437 960
63 156 690All tests are performed on the IBM BlueGene∖Q cluster at the CINECA Centre for HPC, equipped with IBM PowerA2 processors at 1.6 GHz with 10240 nodes, 163 840 computing cores, and 164 Tbytes of internal network RAM. The Fortran 90 code is compiled with the native IBM xlf compiler using the -O3 -qarch=qp -qtune=qp options.


### 5.2. Results on a Fixed Number of Processors

We now compare the performance of the DACG-Newton algorithm for variouskmax values. We tested the proposed algorithm in the computation of the 10 smallest eigenpairs of the afore mentioned large matrices arising from various realistic applications. In all the runs, unless differently specified, we selected the values of the parameters as reported in Table 2.Table 2
Default values of parameters.
Number of eigenpairs to compute
n
eig
=
10
Parameters for the outer iteration
τ
=
1
0
-
8,   ITMAX =50
Tolerance for the initial eigenvector guess
τ
DACG
=
1
0
-
2
Parameters for the PCG iteration
τ
PCG
=
1
0
-
2,   ITMAXPCG=50We will also compute the fill-inρ of the initial preconditioner defined as
(28)ρ=nnz(Gout)+nnz(Gin)nnz(lower  triangular  part  ofA).
In Table 3, we report the parameters for the RFSAI preconditioner which we selected for the four-test problems.Table 3
RFSAI parameters used for the eigensolution of the four test problems.
Matrix name
G
out
n
band
G
in
ρ
δ
d
ε
δ
d
ε
F
LOW3D-663
0.01
4
0.05
1
0.01
2
0.05
1.16
EMILIA-923
0.1
4
0.05
1
0.05
2
0.01
0.74
LONG-1103
0.15
4
0.1
1
0.1
4
0.1
0.29
GEO-1438
0.15
4
0.1
1
0.1
2
0.05
0.34We first analyze the influence of parameterkmax in the convergence of the Newton-DACG method. To this aim we provide, in Table 4 the values of outer iteration number, total matrix-vector products (MVP) and CPU time in evaluating 10 eigenpairs of matrix EMILIA-923 for different values of kmax. From the table, we notice how iteration number and CPU time monotonically decreases with increasing kmax. Moreover, for kmax≤2, the iterative process reached the maximum number of iterations before fulfilling the exit test on the relative residual. Another evidence of the efficiency of the BFGS update is given in Figure 1, where we plot the relative residual norm ∥rk,l∥/(xlTAxl) versus number of linear iterations in converging to eigenpair number 3 of matrix  FLOW3D-663. Note that, with kmax=0, that is, using the constant RFSAI preconditioner, the Newton phase could not reach the 10-8 accuracy within the prescribed ITMAX=50 maximum number of outer iterations.Table 4
Influence of parameterkmax in the number of iterations of the Newton phase. Problem   EMILIA-923  with p=64 processors.
k
max
Outerits
MVP
CPU
Comments
0
431
15035
421.25
ITMAXPCG=100. No convergence within  ITMAX =50. Final (average) residual norm of 10-6
1
438
10899
336.74
No convergence within  ITMAX=50.  Final (average) residual norm of 10-6
2
343
8329
258.22
No convergence within  ITMAX=50. Final (average) residual norm of 10-7
5
267
6839
198.19
10
204
5629
178.04
20
185
5616
164.18Figure 1
Convergence profile of the relative residual norm versus cumulative inner iterations for eigenvalue number 3, matrixFLOW3D-663.In Table5, we report the number of matrix vector products and CPU times for Newton-DACG as compared with “pure” DACG, that is, with DACG run with a final tolerance of 10-8. In every test problem DACG-Newton reveals superior to DACG.Table 5
Comparison between DACG and Newton-DACG for the four test problems onp=64 processors.  Results on problem LONG-1103 are obtained using τDACG=0.5.
Matrix name
DACG
Newton-DACG
MVP
CPU time
k
max
MVP
CPU time
F
LOW3D-663
51427
658.55
10
15833
229.33
EMILIA-923
13404
410.93
20
10649
346.83
LONG-1103
15222
390.05
20
8857
256.59
GEO-1438
12667
354.63
20
10441
326.99In particular, for problemFLOW3D-663, the improvement provided by DACG-Newton is impressive. The DACG method, if run until the final tolerance τ=10-8, is very slow due to the very small relative separation between consecutive eigenvalues ξj=(λj+1-λj)/λj, which drives convergence of this PCG-like solver [16, 25]. In fact the ratio ξj also influences the convergence of Newton’s method being related to the condition number of the Jacobian. However, the (RFSAI-BFGS preconditioned) DACG-Newton algorithm seems less sensitive to this ill-conditioning.


### 5.3. Parallel Results and Scalability

We now analyze the parallel efficiency of the previously described eigensolver preconditioned by RFSAI-BFGS. We will use a strong scaling measure to see how CPU times vary with the number of processors for a fixed total problem size. Denoting withTp the total CPU elapsed times expressed in seconds on p processors, we define relative measures of the parallel efficiency and speedup of our code. We define Sp(p-) as the pseudo speedup computed with respect to the smallest number of processors (p-) used to solve a given problem and Ep(p-) the corresponding efficiency:
(29)Sp(p-)=Tp-p-Tp,Ep(p-)=Sp(p-)p=Tp-p-Tpp.The results obtained with the four test matrices are presented in Tables6, 7, 8, and 9. As a general comment, we may observe that the overall code scales well for the three largest matrices, up to 256 processors (512 for problem GEO-1438). The parallel pseudo efficiency, as expected, decreases with growing number of processors but it is roughly 50% for the largest number of processors employed.Table 6
Iterations, timings, and pseudo efficiencies for problemFLOW3D-663.
p
precond
DACG
Newton
Total
E
p
(
8
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
8
42.43
238.34
4291
704.78
11112
985.55
15403
16
27.00
131.29
4291
408.05
11681
566.33
15972
0.87
32
18.17
79.85
4292
241.79
11275
339.81
15567
0.73
64
89.32
54.57
4291
166.77
11542
230.67
15833
0.53
128
84.89
40.68
4292
119.82
11177
165.39
15469
0.37Table 7
Iterations, timings, and pseudo efficiencies for problemEMILIA-923.
p
precond
DACG
Newton
Total
E
p
(
16
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
16
85.91
478.29
5379
527.05
5616
1091.24
10995
32
55.44
254.57
5377
286.06
5620
596.07
10997
0.92
64
35.47
149.71
5378
164.18
5616
349.36
10994
0.78
128
22.30
88.96
5378
98.61
5618
209.87
10996
0.65
256
13.53
58.98
5378
66.68
5616
139.19
10994
0.49Table 8
Iterations, timings, and pseudo efficiencies for problemLONG-1103.
p
precond
DACG
Newton
Total
E
p
(
32
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
32
24.87
127.47
2734
346.87
5972
499.20
8706
64
15.16
68.12
2734
172.24
5823
255.52
8557
0.97
128
8.57
42.90
2734
109.31
5928
160.78
8662
0.78
256
5.12
31.10
2734
77.41
5806
113.63
8540
0.55Table 9
Iterations, timings, and pseudo efficiencies for problemGEO-1438.
p
precond
DACG
Newton
Total
E
p
(
32
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
32
23.87
289.92
5002
358.20
5432
671.99
10434
64
13.16
139.61
5002
172.14
5439
324.91
10441
1.03
128
7.66
80.15
5002
97.98
5427
185.79
10429
0.90
256
5.38
49.86
5002
61.66
5431
116.90
10433
0.72
512
6.82
37.78
5002
47.31
5429
91.91
10431
0.46Regarding theFLOW3D-663 problem, which is characterized by a relatively small dimension and high sparsity (≈15 nonzeros per row), the parallel efficiency starts to worsen with p=128 processors.


## 5.1. Test Problems

We report the results of our experiments with the RFSAI preconditioner in the solution of a set of problems of large size.The test matrices are realistic examples arising from 3D FE discretization of fluid flow and geomechanical models. Described in detail as follows:(1)
FLOW3D-663 arises from tetrahedral FE discretization of an elliptic PDE describing fluid flow in porous media.(2)
EMILIA-923 arises from the regional geomechanical model of a deep hydrocarbon reservoir. It is obtained discretizing the structural problem with tetrahedral finite elements. Due to the complex geometry of the geological formation, it was not possible to obtain a computational grid characterized by regularly shaped elements.(3)
LONG-1103 is the structural SPD block obtained from a 3D coupled consolidation problem of a geological formation, discretized with tetrahedral finite elements on a higly irregular computational grid.(4)
GEO-1438 arises from a regional geomechanical model of the sedimentary basin underlying the Venice lagoon discretized by a linear FE with randomly heterogeneous properties. 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.Matrices 2 to 4 are publicly available in the University of Florida Sparse Matrix Collection athttp://www.cise.ufl.edu/research/sparse/matrices/. The size and number of nonzero elements for each matrix is provided in Table 1.Table 1
Sizen and number of nonzeros (nnz)  of the test matrices.
Name
n
nnz
F
LOW3D-663
663 390
9 762 648
EMILIA-923
923 136
41 005 206
LONG-1103
1 102 614
48 987 558
GEO-1438
1 437 960
63 156 690All tests are performed on the IBM BlueGene∖Q cluster at the CINECA Centre for HPC, equipped with IBM PowerA2 processors at 1.6 GHz with 10240 nodes, 163 840 computing cores, and 164 Tbytes of internal network RAM. The Fortran 90 code is compiled with the native IBM xlf compiler using the -O3 -qarch=qp -qtune=qp options.


## 5.2. Results on a Fixed Number of Processors

We now compare the performance of the DACG-Newton algorithm for variouskmax values. We tested the proposed algorithm in the computation of the 10 smallest eigenpairs of the afore mentioned large matrices arising from various realistic applications. In all the runs, unless differently specified, we selected the values of the parameters as reported in Table 2.Table 2
Default values of parameters.
Number of eigenpairs to compute
n
eig
=
10
Parameters for the outer iteration
τ
=
1
0
-
8,   ITMAX =50
Tolerance for the initial eigenvector guess
τ
DACG
=
1
0
-
2
Parameters for the PCG iteration
τ
PCG
=
1
0
-
2,   ITMAXPCG=50We will also compute the fill-inρ of the initial preconditioner defined as
(28)ρ=nnz(Gout)+nnz(Gin)nnz(lower  triangular  part  ofA).
In Table 3, we report the parameters for the RFSAI preconditioner which we selected for the four-test problems.Table 3
RFSAI parameters used for the eigensolution of the four test problems.
Matrix name
G
out
n
band
G
in
ρ
δ
d
ε
δ
d
ε
F
LOW3D-663
0.01
4
0.05
1
0.01
2
0.05
1.16
EMILIA-923
0.1
4
0.05
1
0.05
2
0.01
0.74
LONG-1103
0.15
4
0.1
1
0.1
4
0.1
0.29
GEO-1438
0.15
4
0.1
1
0.1
2
0.05
0.34We first analyze the influence of parameterkmax in the convergence of the Newton-DACG method. To this aim we provide, in Table 4 the values of outer iteration number, total matrix-vector products (MVP) and CPU time in evaluating 10 eigenpairs of matrix EMILIA-923 for different values of kmax. From the table, we notice how iteration number and CPU time monotonically decreases with increasing kmax. Moreover, for kmax≤2, the iterative process reached the maximum number of iterations before fulfilling the exit test on the relative residual. Another evidence of the efficiency of the BFGS update is given in Figure 1, where we plot the relative residual norm ∥rk,l∥/(xlTAxl) versus number of linear iterations in converging to eigenpair number 3 of matrix  FLOW3D-663. Note that, with kmax=0, that is, using the constant RFSAI preconditioner, the Newton phase could not reach the 10-8 accuracy within the prescribed ITMAX=50 maximum number of outer iterations.Table 4
Influence of parameterkmax in the number of iterations of the Newton phase. Problem   EMILIA-923  with p=64 processors.
k
max
Outerits
MVP
CPU
Comments
0
431
15035
421.25
ITMAXPCG=100. No convergence within  ITMAX =50. Final (average) residual norm of 10-6
1
438
10899
336.74
No convergence within  ITMAX=50.  Final (average) residual norm of 10-6
2
343
8329
258.22
No convergence within  ITMAX=50. Final (average) residual norm of 10-7
5
267
6839
198.19
10
204
5629
178.04
20
185
5616
164.18Figure 1
Convergence profile of the relative residual norm versus cumulative inner iterations for eigenvalue number 3, matrixFLOW3D-663.In Table5, we report the number of matrix vector products and CPU times for Newton-DACG as compared with “pure” DACG, that is, with DACG run with a final tolerance of 10-8. In every test problem DACG-Newton reveals superior to DACG.Table 5
Comparison between DACG and Newton-DACG for the four test problems onp=64 processors.  Results on problem LONG-1103 are obtained using τDACG=0.5.
Matrix name
DACG
Newton-DACG
MVP
CPU time
k
max
MVP
CPU time
F
LOW3D-663
51427
658.55
10
15833
229.33
EMILIA-923
13404
410.93
20
10649
346.83
LONG-1103
15222
390.05
20
8857
256.59
GEO-1438
12667
354.63
20
10441
326.99In particular, for problemFLOW3D-663, the improvement provided by DACG-Newton is impressive. The DACG method, if run until the final tolerance τ=10-8, is very slow due to the very small relative separation between consecutive eigenvalues ξj=(λj+1-λj)/λj, which drives convergence of this PCG-like solver [16, 25]. In fact the ratio ξj also influences the convergence of Newton’s method being related to the condition number of the Jacobian. However, the (RFSAI-BFGS preconditioned) DACG-Newton algorithm seems less sensitive to this ill-conditioning.


## 5.3. Parallel Results and Scalability

We now analyze the parallel efficiency of the previously described eigensolver preconditioned by RFSAI-BFGS. We will use a strong scaling measure to see how CPU times vary with the number of processors for a fixed total problem size. Denoting withTp the total CPU elapsed times expressed in seconds on p processors, we define relative measures of the parallel efficiency and speedup of our code. We define Sp(p-) as the pseudo speedup computed with respect to the smallest number of processors (p-) used to solve a given problem and Ep(p-) the corresponding efficiency:
(29)Sp(p-)=Tp-p-Tp,Ep(p-)=Sp(p-)p=Tp-p-Tpp.The results obtained with the four test matrices are presented in Tables6, 7, 8, and 9. As a general comment, we may observe that the overall code scales well for the three largest matrices, up to 256 processors (512 for problem GEO-1438). The parallel pseudo efficiency, as expected, decreases with growing number of processors but it is roughly 50% for the largest number of processors employed.Table 6
Iterations, timings, and pseudo efficiencies for problemFLOW3D-663.
p
precond
DACG
Newton
Total
E
p
(
8
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
8
42.43
238.34
4291
704.78
11112
985.55
15403
16
27.00
131.29
4291
408.05
11681
566.33
15972
0.87
32
18.17
79.85
4292
241.79
11275
339.81
15567
0.73
64
89.32
54.57
4291
166.77
11542
230.67
15833
0.53
128
84.89
40.68
4292
119.82
11177
165.39
15469
0.37Table 7
Iterations, timings, and pseudo efficiencies for problemEMILIA-923.
p
precond
DACG
Newton
Total
E
p
(
16
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
16
85.91
478.29
5379
527.05
5616
1091.24
10995
32
55.44
254.57
5377
286.06
5620
596.07
10997
0.92
64
35.47
149.71
5378
164.18
5616
349.36
10994
0.78
128
22.30
88.96
5378
98.61
5618
209.87
10996
0.65
256
13.53
58.98
5378
66.68
5616
139.19
10994
0.49Table 8
Iterations, timings, and pseudo efficiencies for problemLONG-1103.
p
precond
DACG
Newton
Total
E
p
(
32
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
32
24.87
127.47
2734
346.87
5972
499.20
8706
64
15.16
68.12
2734
172.24
5823
255.52
8557
0.97
128
8.57
42.90
2734
109.31
5928
160.78
8662
0.78
256
5.12
31.10
2734
77.41
5806
113.63
8540
0.55Table 9
Iterations, timings, and pseudo efficiencies for problemGEO-1438.
p
precond
DACG
Newton
Total
E
p
(
32
)
CPU
CPU
MVP
CPU
MVP
CPU
MVP
32
23.87
289.92
5002
358.20
5432
671.99
10434
64
13.16
139.61
5002
172.14
5439
324.91
10441
1.03
128
7.66
80.15
5002
97.98
5427
185.79
10429
0.90
256
5.38
49.86
5002
61.66
5431
116.90
10433
0.72
512
6.82
37.78
5002
47.31
5429
91.91
10431
0.46Regarding theFLOW3D-663 problem, which is characterized by a relatively small dimension and high sparsity (≈15 nonzeros per row), the parallel efficiency starts to worsen with p=128 processors.


## 6. Conclusion

We have proposed a parallel RFSAI-BFGS preconditioner for the acceleration of the PCG method in the approximate solution of the linearized Newton systems in the evaluation of a number of the leftmost eigenpairs of large SPD matrices. We have shown that updating an initial preconditioner (here RFSAI) by a low-rank correction using the BFGS formula, produces significant savings in the total number of iterations of the inner solver (the PCG method). The Newton’s algorithm, preconditioner by RFSAI-BFGS, with the aid of a number of initial iterations of DACG to obtain a good initial eigenvector guess, has been completely parallelized and run on the new IBM BlueGene∖Q, located at CINECA, Bologna, Italy. The scalability results are very satisfactory, as compared to the size and nonzeros of the problems selected. In particular, for the largest problem, an efficiency of 72% is obtained with 256 processors. Future research is aimed at investigating the relations between the proposed accelerated Newton method and the well-established Jacobi-Davidson method.



---

*Source: 767042-2013-09-26.xml*