Mitglied der Helmholtz-Gemeinschaft Block Iterative Eigensolvers for Sequences of Dense Correlated Eigenvalue Problems Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli
Motivation and Goals Reverse Simulation total energy ⇐ = Mathematical model band energy gap ← − Simulations conductivity = ⇒ Algorithmic structure forces, etc. Goal Increasing the performance of large legacy codes by exploiting physical information extracted from the simulations that can be used to speed-up the algorithms used in such codes Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 2
Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution A LGORITHM ⇐ S IM – Exploiting eigenvector collinearity: block iterative eigensolvers Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 3
Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution A LGORITHM ⇐ S IM – Exploiting eigenvector collinearity: block iterative eigensolvers Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 4
The Foundations Investigative framework Quantum Mechanics and its ingredients n n Z α 1 h 2 H = − ¯ ∇ 2 ∑ i = 1 ∑ ∑ | x i − a α | + ∑ i − Hamiltonian 2 m | x i − x j | α i = 1 i < j Φ ( x 1 ; s 1 , x 2 ; s 2 ,..., x n ; s n ) Wavefunction � n � R 3 ×{± 1 Φ : 2 } − → R high-dimensional anti-symmetric function – describes the orbitals of atoms and molecules. In the Born-Oppenheimer approximation, it is the solution of the Electronic Schrödinger Equation H Φ ( x 1 ; s 1 , x 2 ; s 2 ,..., x n ; s n ) = E Φ ( x 1 ; s 1 , x 2 ; s 2 ,..., x n ; s n ) Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 5
Density Functional Theory (DFT) 1 Φ ( x 1 ; s 1 , x 2 ; s 2 ,..., x n ; s n ) = ⇒ Λ i , a φ a ( x i ; s i ) 2 density of states n ( r ) = ∑ a | φ a ( r ) | 2 3 In the Schrödinger equation the exact Coulomb interaction is substituted with an effective potential V 0 ( r ) = V I ( r )+ V H ( r )+ V xc ( r ) Hohenberg-Kohn theorem ∃ one-to-one correspondence n ( r ) ↔ V 0 ( r ) = ⇒ V 0 ( r ) = V 0 ( r )[ n ] ∃ ! a functional E [ n ] : E 0 = min n E [ n ] The high-dimensional Schrödinger equation translates into a set of coupled non-linear low-dimensional self-consistent Kohn-Sham (KS) equation � � h 2 − ¯ 2 m ∇ 2 + V 0 ( r ) ˆ ∀ a H KS φ a ( r ) = φ a ( r ) = ε a φ a ( r ) solve Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 6
Kohn-Sham scheme and DFT Self-consistent cycle Typically this set of equations is solved using an iterative self-consistent cycle Solve KS equations Initial guess Compute KS potential ˆ n init ( r ) = ⇒ V 0 ( r )[ n ] − → H KS φ a ( r ) = ε a φ a ( r ) ↑ No ↓ Converged? Compute new density OUTPUT Yes | n ( ℓ + 1 ) − n ( ℓ ) | < η n ( r ) = ∑ a | φ a ( r ) | 2 ⇐ = ← − Energy, forces, ... In practice this iterative cycle is much more computationally challenging and requires some form of broadly defined discretization Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 7
Generalized eigenvalue problems Ax = λ Bx A common way of discretizing the KS equations is to expand the wave functions φ a ( r ) on a basis set ∑ c G φ a ( r ) − → φ k , ν ( r ) = k , ν ψ G ( k , r ) | G + k |≤ K max This expansion is then inserted in the KS equations H KS c G ′ c G ′ ψ ∗ k , ν ψ G ′ ( k , r ) = λ k ν ψ ∗ G ( k , r ) ∑ ˆ G ( k , r ) ∑ k , ν ψ G ′ ( k , r ) , G ′ G ′ and, by defining the matrix entries for the left and right hand side respectively as Hamiltonian A k and overlap matrices B k , � ˆ � [ A k B k ] GG ′ = ∑ d r ψ ∗ H KS ˆ � G ( k , r ) ψ G ′ ( k , r ) 1 α one arrives at generalized eigenvalue equations parametrized by k ( A k ) GG ′ c G ′ ( B k ) GG ′ c G ′ ∑ k ν = λ k ν ∑ ≡ A k x i = λ i B k x i . P k : k ν . G ′ G ′ Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 8
Discretized Kohn-Sham scheme Self-consistent cycle Solve a set of Initial guess Compute KS potential eigenproblems n start ( r ) = ⇒ V 0 ( r )[ n ] − → P k 1 ... P k N ↑ No ↓ Converged? Compute new density OUTPUT Yes | n ( ℓ + 1 ) − n ( ℓ ) | < η n ( r ) = ∑ k , ν | φ k , ν ( r ) | 2 ⇐ = ← − Energy, ... Observations: 1 A and B are respectively hermitian and hermitian positive definite 2 eigenproblems across k index have different size and we consider them independent from each other (for the moment) 3 eigenvectors of problems of same k are seemingly uncorrelated across iterations i 4 k = 1:10-100 ; i = 1:20-50 Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 9
Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution A LGORITHM ⇐ S IM – Exploiting eigenvector collinearity: block iterative eigensolvers Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 10
Sequences of eigs: an A LGORITHM ⇐ S IM case Sequences of eigenproblems Consider the set of generalized eigenproblems P ( 1 ) ... P ( ℓ ) P ( ℓ + 1 ) ... P ( N ) � = ( P ) N � P ( ℓ ) � Could this sequence of eigenproblems evolve following a convergence pattern in line with the convergence of n ( r ) ? Numerical study: studied the evolutions of the angles b/w eigenvectors of successive iterations developed a method that establishes systematically a one-to-one correspondence b/w eigenvectors collected data on eigenvectors deviation angles 1 analyzed deviation angles at fixed λ for all k s 2 analyzed deviation angles at fixed k for all λ s below Fermi Energy Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 11
Angle evolution fixed k Example: a metallic compound at fixed k Evolution of subspace angle for eigenvectors of k − point 1 and lowest 75 eigs 0 10 AuAg Angle b/w eigenvectors of adjacent iterations − 2 10 − 4 10 − 6 10 − 8 10 − 10 10 2 6 10 14 18 22 Iterations (2 − > 22) Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 12
Correlation and its exploitation ∃ correlation between successive eigenvectors x ( ℓ − 1 ) and x ( ℓ ) Angles decrease monotonically with some oscillation Majority of angles are small after the first few iterations Note: Mathematical model � correlation. Correlation ⇐ systematic analysis of the simulation . A LGORITHM ⇐ S IM The stage is favorable to an iterative eigensolver where the eigenvectors of P ( ℓ − 1 ) are fed to the solve P ( ℓ ) Next stages of the investigation: 1 Establish which eigensolvers can exploit the evolution (Implicit Restarted Arnoldi, Krylov-Schur, Subspace Iteration, Davidson-like, etc.) 2 Investigate if approximate eigenvectors can speed-up iterative solvers 3 Understand if iterative methods be competitive with direct methods for dense problems Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 13
Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution A LGORITHM ⇐ S IM – Exploiting eigenvector collinearity: block iterative eigensolvers Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 14
Block Iterative Eigensolvers Two essential properties the iterative algorithms have to comply with: 1 the ability to receive as input a sizable set of approximate eigenvectors; 2 the capacity to solve simultaneously for a substantial portion of eigenpairs. Block iterative methods constitutes the natural choice: they accept a variable set of multiple starting vectors; these methods have a faster convergence rate and avoid stalling when facing small clusters of eigenvalues; when augmented with polynomial accelerators their performance is further improved. ALGORITHMS: Block Krylov-Schur – Block Chebyshev-Davidson – Chebyshev Subspace Iteration – LOBPCG Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 15
Experimental tests setup Matrix sizes: 2,000 ÷ 5,700 Num of fixed-point iterations: 15 ÷ 30 Num of k-points: 6 ÷ 27 B ill-conditioned B is in general almost singular. Examples: size ( A ) = 50 → κ ( A ) ≈ 10 4 size ( A ) = 500 → κ ( A ) ≈ 10 7 We used the standard form for the problem A ′ = L − 1 AL − T A ′ y = λ y y = L T x Ax = λ Bx − → with and Numerical Study: Approx. vs Random solutions against Iteration Index Approx. vs Random solutions against Spectrum Fraction Numerical tests were performed using Matlab version R2011b (7.13.0.564) running on an Intel i7 CPU with 8 cores at 2.93 GHz. Four cores and 8 Gb of RAM were fully dedicated to computations Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 16
Speed − up vs. iterations for CaFe 2 As 2 (n=2612) 100 90 ChSI 80 70 BChDav 3% BChDav 7% ChSI 3% Speed − up (%) 60 Lobpcg ChSI 7% Lobpcg 3% 50 Lobpcg 7% 40 30 20 BChDav 10 0 3 7 11 14 17 20 23 Iteration Index Figure: Comparison between the 3 most effective block iterative eigensolvers for CaFe 2 As 2 with respect to the outer-iteration index Birkbeck University, London, June the 29th 2012 Edoardo Di Napoli Folie 17
Recommend
More recommend
