SLIDE 1 Preconditioned Locally Harmonic Residual Methods for Interior Eigenvalue Computations
Eugene Vecharynski Lawrence Berkeley National Laboratory Computational Research Division
−0.4 −0.3 −0.2 −0.1 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 Re Im
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.34 0.34 0.79 0.79 0.00 0.00 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.23 0.12 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.23 0.12 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.29 0.15 0.00 0.39 0.24 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.29 0.15 0.00 0.39 0.24 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.09 0.06 0.00 0,44 0.81 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.09 0.06 0.00 0,44 0.81 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.01 0.00 0.00 0.21 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.01 0.00 0.00 0.21 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.00 0.01 0.00 0.03 0.02 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 0.00 0.01 0.00 0.03 0.02 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.34 0.34 0.79 0.79 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.29 0.00 0.29 0.00 0.92 0.92 0.00 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.00 0.09 0.09 0.00 0.00 0.00 0.00 0.00 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.00 0.00 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.00 0.01 0.00 0.00 0.00 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 0.00 0.01 0.00 0.00 0.00 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
SIAM Conference on Applied Linear Algebra, October 29, 2015 Atlanta, GA
SLIDE 2
The Preconditioned Locally Harmonic Residual methods
A new class of robust block preconditioned iterative eigensolvers for computing interior eigenpairs
I Hermitian eigenproblems
The PLHR algorithm (EV, A. Knyazev)
I Non-Hermitian eigenproblems
The Generalized PLHR (GPLHR) algorithm (EV, C. Yang, F. Xue)
SLIDE 3 Problem
Compute a subset of k right eigenpairs (λ, x) of a non-Hermitian matrix pair (A, B) that are closest to a given shift σ ∈ C Ax = λBx, A, B ∈ Cn×n
−0.4 −0.3 −0.2 −0.1 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 Re Im
A, B can be extremely large A, B may not be explicitly given (A, B) is regular σ can point to the spectrum’s interior or exterior
SLIDE 4
Many origins of the problem
Study of excited states of molecules (complex scaling, EOM-CC) Chemical reactions, stability analysis in fluid dynamics, crystal growth simulations, power systems, Markov chains, plasmas, ETC.
SLIDE 5
Available eigensolvers
I (Inexact) Inverse subspace iteration
Sufficiently accurate (A − σB)−1 is needed
I Arnoldi (ARPACK)
Sufficiently accurate (A − σB)−1 is needed
I Block Generalized Davidson (Morgan-92)
Robustness issues, may need a larger search subspace
I JDQR/JDQZ (Fokkema et. al. SISC-98)
“One-by-one” eigenpair computation, lack of BLAS3
SLIDE 6
Need an eigensolver that ...
I Avoids “shift-and-invert” I Takes advantage of preconditioning, T ≈ (A − σB)−1 I Performs block iteration (parallel performance, BLAS3) I Is robust and reliable if memory is limited I Handles standard and generalized, interior and exterior,
eigenproblems in a uniform manner
I Is simple (to the extent a non-Hermitian solver can be ;-))
SLIDE 7
An example from the “Hermitian world”: LOBPCG
The LOBPCG algorithm (Knyazev SISC-01): X (i+1) ← col n X (i), T(AX (i) − BX (i)Λ(i)), X (i−1)o
I Block iterations I No inversions I Short-term recurrence I The Rayleigh–Ritz procedure over low-dimensional subspaces
SLIDE 8
An example from the “Hermitian world”: LOBPCG
The LOBPCG algorithm (Knyazev SISC-01): X (i+1) ← col n X (i), T(AX (i) − BX (i)Λ(i)), X (i−1)o
I Block iterations I No inversions I Short-term recurrence I The Rayleigh–Ritz procedure over low-dimensional subspaces
Does NOT extend to non-Hermitian eigenproblems
SLIDE 9 Difficulties
- 1. The eigendecomposition AX = BXΛ may not exist or X can
be ill-conditioned. Partial generalized Schur form should be computed
⇢ AV = QRA BV = QRB , λj = RA(j, j)/RB(j, j)
- 2. A correct definition of a block residual is not clear
- 3. The LOBPCG-like trial subspace is “too small”
- 4. The standard Rayleigh–Ritz may not be appropriate for
interior eigenvalue computations
SLIDE 10
The “Q-free” partial generalized Schur form
I Partial generalized Schur form
⇢ AV = QRA BV = QRB , λj = RA(j, j)/RB(j, j)
I The “Q-free” partial generalized Schur form
AVMB = BVMA, λj = MA(j, j)/MB(j, j), where MA, MB are upper triangular. – Gives a numerically stable analogue of AXΛB = BXΛA – Allows defining a block residual AVMB − BVMA
SLIDE 11
The “Q-free” partial generalized Schur form
I Partial generalized Schur form
⇢ AV = QRA BV = QRB , λj = RA(j, j)/RB(j, j)
I The “Q-free” partial generalized Schur form
AVMB = BVMA, λj = MA(j, j)/MB(j, j), where MA, MB are upper triangular. – Gives a numerically stable analogue of AXΛB = BXΛA – Allows defining a block residual AVMB − BVMA Does the “Q-free” form exist for any regular pair (A, B)?
SLIDE 12 The “Q-free” partial generalized Schur form
For any regular pair (A, B), there exist AV MB = BV MA, λj = MA(j, j)/MB(j, j) where MA = G2G −1RA, MB = I − G1G −1RA, with an upper triangular G = RAG1 + RBG2, G(j, j) = 1, and G1, G2 diagonal
G1(j, j) = 8 < : 0, |RA(j, j)| < |RB(j, j)| 1 − RB(j, j) RA(j, j) ,
G2(j, j) = ⇢ 1/RB(j, j), |RA(j, j)| < |RB(j, j)| 1,
SLIDE 13
The simplest GPLHR trial subspace
Given approximations V (i), R(i)
A , and R(i) B , the “Q-free” partial
generalized Schur form allows
I Defining a preconditioned residual
W (i) = T(AV (i)M(i)
B − BV (i)M(i) A ),
where T is a preconditioning operator.
I Constructing a LOBPCG-like trial subspace
Z = col{V (i), W (i), P(i)} where P(i) is a block of additional search directions.
SLIDE 14
The simplest GPLHR trial subspace
Given approximations V (i), R(i)
A , and R(i) B , the “Q-free” partial
generalized Schur form allows
I Defining a preconditioned residual
W (i) = T(AV (i)M(i)
B − BV (i)M(i) A ),
where T is a preconditioning operator.
I Constructing a LOBPCG-like trial subspace
Z = col{V (i), W (i), P(i)} where P(i) is a block of additional search directions. How to construct a larger trial subspace?
SLIDE 15 The GPLHR trial subspace
I Find the corrections, such that
A(V (i)+C (i))(M(i)
B +∆(i) B ) = B(V (i)+C (i))(M(i) A +∆(i) A ),
V (i)∗C (i) = 0
SLIDE 16 The GPLHR trial subspace
I Find the corrections, such that
A(V (i)+C (i))(M(i)
B +∆(i) B ) = B(V (i)+C (i))(M(i) A +∆(i) A ),
V (i)∗C (i) = 0
I The correction C (i) can be approximated by C that satisfies
the generalized Sylvester equation
L(C) = F, V (i)∗C = 0,
where
L(C) ≡ (PQ⊥APV ⊥)CM(i)
B − (PQ⊥BPV ⊥)CM(i) A
and
F ≡ −PQ⊥(AV (i)M(i)
B − BV (i)M(i) A ).
Here PV ⊥ = I − V (i)V (i)∗, PQ⊥ = I − Q(i)Q(i)∗, and Q(i) is an
- rthonormal basis of col{(A − σB)V (i)}.
SLIDE 17 The GPLHR trial subspace
Define a preconditioner
TL = (I − VV ∗)T(I − QQ∗), T ≈ (A − σB)−1
for the generalized Sylvester’s operator
L(C) ≡ (PQ⊥APV ⊥)CM(i)
B − (PQ⊥BPV ⊥)CM(i) A
Apply m steps of a preconditioned block Arnoldi to L(C) = F
1: Set V ← V (i), MA ← M(i)
A , and MB ← M(i) B ;
2: Set W ← (I − VV ∗)T(I − QQ∗)(AVMB − BVMA);(≡ TL(F)) 3: W ← orth(W ); S0 ← W ; S ← [ ]; 4: for l = 1 → m do 5: Sl ← (I − VV ∗)T(I − QQ∗)(ASl−1MB − BSl−1MA); (≡ TL(L(Sl−1))) 6: Sl ← Sl − W (W ∗Sl); Sl ← Sl − S(S∗Sl); 7: Sl ← orth(Sl); S ← [S Sl]; 8: end for
⇒ The GPLHR trial subspace
Z = col{V (i), W (i), S(i)
1 , S(i) 2 , . . . , S(i) m
| {z }
captures C
, P(i)}.
SLIDE 18
The harmonic Schur–Rayleigh–Ritz (SRR)
I Find V , MB, and MA, such that (Petrov–Galerkin)
AVMB − BVMA ⊥ (A − σB)Z, λj ≈ MA(j, j)/MB(j, j)
where columns of V are in Z .
SLIDE 19
The harmonic Schur–Rayleigh–Ritz (SRR)
I Find V , MB, and MA, such that (Petrov–Galerkin)
AVMB − BVMA ⊥ (A − σB)Z, λj ≈ MA(j, j)/MB(j, j)
where columns of V are in Z .
I Solve the projected problem (the ordered QZ algorithm)
(U∗AZ)YMB = (U∗BZ)YMA,
where Z = orth{Z} and U = orth{(A − σI)Z}.
SLIDE 20
The harmonic Schur–Rayleigh–Ritz (SRR)
I Find V , MB, and MA, such that (Petrov–Galerkin)
AVMB − BVMA ⊥ (A − σB)Z, λj ≈ MA(j, j)/MB(j, j)
where columns of V are in Z .
I Solve the projected problem (the ordered QZ algorithm)
(U∗AZ)YMB = (U∗BZ)YMA,
where Z = orth{Z} and U = orth{(A − σI)Z}.
I Set M(i+1) A
← MA, M(i+1)
B
← MB, and V (i+1) ← ZY .
SLIDE 21
The GPLHR algorithm
Input: A, B, σ, V (0), T ≈ (A − σI)−1, m. Output: k eigenpairs of (A, B) closest to σ
1: Use V (0) to construct initial V , MA, and MB; P ← [ ]; 2: while convergence not reached do 3:
W ← (I − VV ∗)T(I − QQ∗)(AVMB − BVMA); W ← orth(W );
4:
Apply m block Arnoldi steps to generate blocks S1, . . . , Sm;
5:
Set Z ← orth([V , W , S1, . . . , Sm, P]);
6:
Set U ← orth((A − σB)Z);
7:
[¯ RA, ¯ RB, ¯ YL, ¯ YR] ← ordqz(U∗AZ, U∗BZ, σ);
8:
Y ← ¯ YR(:,1:k), YL ← ¯ YL(:,1:k);
9:
RA ← ¯ RA(1:k,1:k), RB ← ¯ RB(1:k,1:k);
10:
V ← ZY ; Q ← UYL; P ← Z ¯ YR(:,k+1:2k);
11:
G ← RAG1 + RBG2; MA ← G2G −1RA; MB ← I − G1G −1RA;
12: end while 13: Extract wanted eigenpairs (X, Λ) using RR for (A, B) over col{V }.
SLIDE 22
The GPLHR algorithm
Input: A, B, σ, V (0), T ≈ (A − σI)−1, m. Output: k eigenpairs of (A, B) closest to σ
1: Use V (0) to construct initial V , MA, and MB; P ← [ ]; 2: while convergence not reached do 3:
W ← (I − VV ∗)T(I − QQ∗)(AVMB − BVMA); W ← orth(W );
4:
Apply m block Arnoldi steps to generate blocks S1, . . . , Sm;
5:
Set Z ← orth([V , W , S1, . . . , Sm, P]);
6:
Set U ← orth((A − σB)Z);
7:
[¯ RA, ¯ RB, ¯ YL, ¯ YR] ← ordqz(U∗AZ, U∗BZ, σ);
8:
Y ← ¯ YR(:,1:k), YL ← ¯ YL(:,1:k);
9:
RA ← ¯ RA(1:k,1:k), RB ← ¯ RB(1:k,1:k);
10:
V ← ZY ; Q ← UYL; P ← Z ¯ YR(:,k+1:2k);
11:
G ← RAG1 + RBG2; MA ← G2G −1RA; MB ← I − G1G −1RA;
12: end while 13: Extract wanted eigenpairs (X, Λ) using RR for (A, B) over col{V }.
SLIDE 23 Choice of the additional search directions P
Problem P(i)
thick
P(i)
LOBPCG
W/o P(i) A15428 12 14 16 AF23560 10 DNC 153 CRY10000 27 DNC DNC DW8192 118 55 329 LSTAB NS 38 DNC DNC MHD4800 137 DNC DNC PDE2961 13 DNC DNC QH882 15 DNC 30 RDB3200L 10 12 12 UTM1700 16 23 21 Table: Numbers of iterations performed by GPLHR to compute k = 5 eigenpairs with different choices of the search directions P(i); m = 1.
SLIDE 24
The GPLHR algorithm
Input: A, B, σ, V (0), T ≈ (A − σI)−1, m. Output: k eigenpairs of (A, B) closest to σ
1: Use V (0) to construct initial V , MA, and MB; P ← [ ]; 2: while convergence not reached do 3:
W ← (I − VV ∗)T(I − QQ∗)(AVMB − BVMA); W ← orth(W );
4:
Apply m block Arnoldi steps to generate blocks S1, . . . , Sm;
5:
Set Z ← orth([V , W , S1, . . . , Sm, P]);
6:
Set U ← orth((A − σB)Z);
7:
[¯ RA, ¯ RB, ¯ YL, ¯ YR] ← ordqz(U∗AZ, U∗BZ, σ);
8:
Y ← ¯ YR(:,1:k), YL ← ¯ YL(:,1:k);
9:
RA ← ¯ RA(1:k,1:k), RB ← ¯ RB(1:k,1:k);
10:
V ← ZY ; Q ← UYL; P ← Z ¯ YR(:,k+1:2k);
11:
G ← RAG1 + RBG2; MA ← G2G −1RA; MB ← I − G1G −1RA;
12: end while 13: Extract wanted eigenpairs (X, Λ) using RR for (A, B) over col{V }.
SLIDE 25
Influence of the projectors
Problem With proj. W/o proj. A15428 14 DNC AF23560 9 9 CRY10000 13 17 DW8192 44 45 LSTAB NS 49 84 MHD4800 12 12 PDE2961 10 11 RDB3200L 12 21 UTM1700 25 26
Table: Numbers of iterations performed by GPLHR with and without the projection at the preconditioning step to compute k = 10 eigenpairs; m = 1.
SLIDE 26
Few algorithmic details
I The preconditioner T ≈ (A − σB)−1 can be given by a
general procedure, e.g., by several steps of a preconditioned iterative linear solver applied to (A − σB)w = r.
I The dimension of the search subspace is (m + 3)k, where k is
the block size. In practice, m is between 1 and 5.
I The algorithms can be implemented with m + 1 matrix-block
multiplications and m + 1 block preconditioning operations per iteration.
I Converged Schur vectors should be locked. I The parameter m can be adaptively chosen, e.g., increased
after a number of eigenpairs vectors has converged
SLIDE 27
Numerical examples: benchmark systems
Experiments in Q-Chem 4.2: compare GPLHR and Davidson with the EOM-IP method for calculation of ionized potentials Benchmark systems:
I Hydrated photoactive yellow protein chromophore PYPa-Wp
(left). The 6-31+G(d,p) basis set (292 basis functions).
I Dihydrated 1,3-dimethyluracil (mU)2-(H2O)2 (right) The
6-311+G(d,p)basis set (336 basis functions).
SLIDE 28
Davidson vs GPLHR: iteration count
Figure: Left: PYPa-Wp/6-31+G(d,p), convergence for the eigenvalues 4.11 and 4.20 eV; Right: (mU)2-(H2O)2/6-311+G(d,p), convergence for eigenvalues 8.89 and 10.04 eV
I The GPLHR iteration count decreases as m grows I GPLHR requires less iterations than Davidson I Each GPLHR iterations requires more (nev ×(m + 1))
matvecs per iteration
SLIDE 29 Davidson vs GPLHR: low-lying eigenvalues (σ = 0a.u.)
(mU)2-(H2O)2/6-311+G(d,p) Davidson nev niters
# matvec 1 6 14 7 3 8 48 24 5 8 76 38 10 15 120 99 12 13 120 118 GPLHR nev niters m
# matvec 1 4 1 8 9 3 4 1 24 27 5 4 1 40 43 10 7 1 80 118 12 5 2 120 151
SLIDE 30 Davidson vs GPLHR: the interior eigenvalues (σ = 11a.u.)
PYPa-Wp/6-31+G(d,p) Davidson nev niters
# matvec 1 DNC – – 2 DNC – – 3 DNC – – 5 DNC – – GPLHR nev niters m
# matvec 1 4 1 8 9 2 4 1 16 18 3 4 1 24 27 5 8 1 40 63
SLIDE 31 Block Generalized Davidson vs GPLHR: other examples
200 400 600 800 1000 10
−8
10
−6
10
−4
10
−2
CRY10000, k=5 Largest relative eigenresidual norm Number of preconditioned matvecs GPLHR BGD(4k) BGD(6k) BGD(8k) 0.5 1 1.5 2 x 10
4
10
−8
10
−6
10
−4
10
−2
10 LSTAB_NS, k=7 Largest relative eigenresidual norm Number of preconditioned matvecs GPLHR BGD(4k) BGD(6k) BGD(8k)
Figure: σ = 8 for CRY10000 (left) and σ = 0 for LSTAB NS (right)
SLIDE 32
JDQZ vs GPLHR
GPLHR JDQZ Problem k itG #it #mv #prec #it #mv #prec MHD4800 1 20 41 40 DNC DNC DNC 3 15 79 76 DNC DNC DNC 5 12 101 96 DNC DNC DNC 10 9 165 155 DNC DNC DNC LSTAB NS 1 25 5 261 260 DNC DNC DNC 3 25 49 5307 5304 DNC DNC DNC 5 25 28 4893 4888 DNC DNC DNC 10 25 36 12984 12974 DNC DNC DNC UTM1700 1 15 5 161 160 MCV MCV MCV 3 15 13 963 960 16 263 280 5 15 13 1669 1664 31 511 545 10 15 21 5034 5024 166 2691 2865
SLIDE 33 JDQZ vs GPLHR
−0.4 −0.3 −0.2 −0.1 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 Twenty eigenvalues of MHD4800 computed by GPLHR Re Im Exact GPLHR Shift −0.4 −0.3 −0.2 −0.1 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 Twenty eigenvalues of MHD4800 computed by JDQZ Re Im Exact JDQZ Shift
SLIDE 34
ARPACK vs GPLHR
Problem k GPLHR ARPACK A15876 10 194 2895 A15428 10 335 3006 AF23560 5 2 8 Table: Time required by GPLHR with the preconditioner T given by ILU(10−3) and ARPACK with the shift-and-invert operator computed through the LU decomposition of A − σI to find k eigenvalues closest to σ and their associated eigenvectors.
SLIDE 35 The PLHR algorithm for Hermitian problems
I A “shorter” trial subspace I Requires SPD preconditioning I Can incorporate preconditioning directly into the harmonic RR
50 100 150 200 250 300 10
−6
10
−4
10
−2
10
Si2H4, σ = 0.5 maxj ||A vj
(i) − λj (i) vj (i)||
Iteration number, i
BPLHR BPLHR−HARM BGD(80) 50 100 150 200 250 300 350 10
−6
10
−4
10
−2
10
Si2H4, σ = 0.7 maxj ||A vj
(i) − λj (i) vj (i)||
Iteration number, i
BPLHR BPLHR−HARM BGD(70)
EV and A.Knyazev, to appear in SISC Copper Mountain Special Issue
SLIDE 36
Conclusions
I GPLHR: a new block preconditioned eigensolver for
non-Hermitian problems;
I Robust and efficient under limited memory; I Solves standard and generalized eigenproblems in a uniform
manner;
I Well suited for high-performance parallel computations; I A pilot version implemented in Q-Chem
SLIDE 37
Thank you!
I EV, C. Yang, and F. Xue: “Generalized preconditioned locally harmonic
residual method for non-Hermitian eigenproblems”, submitted, available at http://arxiv.org/abs/1506.06829
I D. Zuev, EV, C. Yang, N. Orms, and A. Krylov: “New algorithms for
iterative matrix-free eigensolvers in quantum chemistry”, J. Comp. Chem., 2015
I EV and A. Knyazev: “Preconditioned Locally Harmonic Residual Method
for Computing Interior Eigenpairs of Certain Classes of Hermitian Matrices”, to appear in SISC, available at http://arxiv.org/abs/1408.0042
SLIDE 38 Test problems
Problem Type n σ Preconditioner
A15428 stand. 15428 −40.871 ILU(10−3) interior AF23560 stand. 23560 −40 + 300i ILU(10−3) largest mod. CRY10000 stand. 10000 8.0 ILU(10−3) largest Re DW8192 stand. 8192 1.0 ILU(10−3) rightmost LSTAB NS gen. 12619 GMRES(25)+LSC rightmost MHD4800 gen. 4800 −0.1 + 0.5i (A − σB)−1 interior PDE2961 stand. 2961 10.0 ILU(10−3) largest Re QH882 stand. 882 −150 + 180i GMRES(5)+ILU(10−5) interior RDB3200L stand. 3200 2i ILU(10−3) rightmost UTM1700 gen. 1700 GMRES(15)+ILU(10−4) leftmost
Table: Test problems.
