Parallel block Chebyshev subspace iteration algorithm optimized for - PowerPoint PPT Presentation

Mitglied der Helmholtz-Gemeinschaft Parallel block Chebyshev subspace iteration algorithm optimized for sequences of correlated dense eigenproblems ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. 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 ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 2

Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) Exploiting approximate eigenvectors: numerical results ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 3

Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) Exploiting approximate eigenvectors: numerical results ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. 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 ) ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. 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 ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 6

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, ... FLAPW details Observations: 1 every P k : Ax = B λ x is a generalized eigenvalue problem; 2 A and B are DENSE and hermitian (B is also pos. def.); 3 P k s with different k index have different size and are independent from each other. 4 k = 1:10-100 ; i = 1:20-50 ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 7

Algorithmic digression Direct solvers. Iterative solvers. ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 8

Algorithmic digression Direct solvers. Iterative solvers.  ∗ ∗ ∗ ∗ ∗ ∗  ∗ ∗ ∗ ∗ ∗ ∗   ∗ ∗ ∗ ∗ ∗ ∗     ∗ ∗ ∗ ∗ ∗ ∗     ∗ ∗ ∗ ∗ ∗ ∗   ∗ ∗ ∗ ∗ ∗ ∗  ∗ ∗  ∗ ∗ ∗   ∗ ∗ ∗     ∗ ∗ ∗     ∗ ∗ ∗   ∗ ∗ ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 8

Algorithmic digression Direct solvers. Iterative solvers. | λ 1 | > | λ 2 | > | λ 3 | > ...  ∗ ∗ ∗ ∗ ∗ ∗  ∗ ∗ ∗ ∗ ∗ ∗   Ax j = λ j x j ∗ ∗ ∗ ∗ ∗ ∗     ∗ ∗ ∗ ∗ ∗ ∗     v = ∑ j γ j x j ∗ ∗ ∗ ∗ ∗ ∗   ∗ ∗ ∗ ∗ ∗ ∗ Av = ∑ j λ j γ j x j ⇒ A k v = ∑ j λ k j γ j x j  ∗ ∗  � � λ 1 � � ∗ ∗ ∗ Rate of convergence → magnitude of � �   ∗ ∗ ∗ � λ j �     � � ∗ ∗ ∗     ∗ ∗ ∗   ∗ ∗ ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 8

Algorithmic digression Direct solvers. Iterative solvers. Sparse matrices. Dense matrices. ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 8

Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) Exploiting approximate eigenvectors: numerical results ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 9

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 ) not a a set of disjoint problems ( P ) N , but as a sequence; � P ( ℓ ) � Could this sequence of eigenproblems evolve following a convergence pattern in line with the convergence of n ( r ) ? R EVERSE S IMULATION method: numerical simulations analyzed employing a parameter-based “inverse” problem method; collected data on deviation angles b/w eigenvectors of adjacent eigenproblems; identified “evolutions” of eigenvectors along the sequence. ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 10

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) ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 11

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 ⇐ numerical 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 Development of a block iterative eigensolver that can exploit the correlation 2 Investigate if approximate eigenvectors can speed-up the iterative solver of choice 3 Understand if such an iterative method be competitive with direct methods for dense problems ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 12

Algorithmic choice Direct solvers. Iterative solvers. Sparse matrices. Dense matrices. ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 13

Outline Stating the problem: how sequences of generalized eigenproblems arise in all-electron computations Eigenvectors angle evolution The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) Exploiting approximate eigenvectors: numerical results ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 14

Chebyshev Filtered Sub-space Iteration method Two essential properties the iterative algorithm has 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. ChFSI constitutes the natural choice: it accepts the full set of multiple starting vectors; it avoids stalling when facing small clusters of eigenvalues; when augmented with polynomial accelerators it has a much faster convergence rate; converged eigenvectors can be easily locked. ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 15

Pseudocode I NPUT : Hamiltonian, approximation for the eigenpairs – ( Λ , W ) , TOL , DEG . O UTPUT : Wanted eigenpairs. 1 Lanczos step. Identify the bounds for the interval to be filtered out. R EPEAT U NTIL CONVERGENCE : 2 Chebyshev filter. Filter a block of vectors W . 3 QR decomposition. Re-orthogonalize the vectors outputted by the filter. 4 Compute the Rayleigh quotient G = W H HW . 5 Compute the primitive Ritz pairs ( Λ , Q ) . 6 Compute the approximate Ritz pairs ( Λ , WQ ) . 7 Check which one among the Ritz vectors converged . 8 Deflate and lock the converged vectors. E ND R EPEAT ERCIM 2012 Oviedo, Spain , Dec. 2nd M. Berljafa and E. Di Napoli Folie 16