An optimized subspace iteration eigensolver applied to sequences of - PowerPoint PPT Presentation

Mitglied der Helmholtz-Gemeinschaft An optimized subspace iteration eigensolver applied to sequences of dense eigenproblems in ab initio simulations PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa

Motivation: electronic structure computation Full-potential Linearized Augmented Plane Waves (FLAPW) self-consistent field cycle Solve a set of Initial guess Compute discretized eigenproblems for charge density Kohn-Sham P ( ℓ ) k 1 ... P ( ℓ ) n start ( r ) equations k N No Compute new OUTPUT Yes Converged? Electronic charge density | n ( ℓ ) − n ( ℓ − 1 ) | < η structure, n ( ℓ ) ( r ) ... 1 every P ( ℓ ) : A ( ℓ ) k x = B ( ℓ ) k λ x is a generalized eigenvalue problem; k 2 A and B are DENSE and hermitian (B is positive definite); 3 required: lower 2 ÷ 10 % of eigenpairs; 4 momentum vector index: k = 1 : 10 ÷ 100 ; 5 iteration cycle index: ℓ = 1 : 20 ÷ 50 . PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 2

Outline Sequences of correlated eigenproblems The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) ChFSI parallelization and numerical tests PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 3

Outline Sequences of correlated eigenproblems The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) ChFSI parallelization and numerical tests PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 4

Sequences of eigenproblems Definitions and solving strategies Definition: Eigenproblem sequence A sequence of eigenproblems is a finite and index-ordered set of problems { P } N . = P ( 1 ) ··· P ( ℓ ) ··· P ( N ) with same size = n such that the eigenpairs of P ( ℓ ) are used (directly or indirectly) to initialize P ( ℓ + 1 ) . PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 5

Sequences of eigenproblems Definitions and solving strategies Definition: Eigenproblem sequence A sequence of eigenproblems is a finite and index-ordered set of problems { P } N . = P ( 1 ) ··· P ( ℓ ) ··· P ( N ) with same size = n such that the eigenpairs of P ( ℓ ) are used (directly or indirectly) to initialize P ( ℓ + 1 ) . Current solving strategy The set of generalized eigenproblems P ( 1 ) ... P ( ℓ ) P ( ℓ + 1 ) ... P ( N ) is handled as a set of disjoint problems N × P ; Each problem P ( ℓ ) is solved independently using a direct solver as a black-box from a standard library (i.e. ScaLAPACK). PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 5

Sequences of Eigenproblems Adjacent iteration cycles ITERATION ( ℓ ) ITERATION ( ℓ + 1 ) direct direct P ( ℓ ) ( X ( ℓ ) k 1 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 1 ) ) k 1 k 1 k 1 k 1 solver solver direct direct P ( ℓ ) ( X ( ℓ ) k 2 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 2 ) ) k 2 k 2 k 2 k 2 solver solver Next cycle direct direct P ( ℓ ) ( X ( ℓ ) k N , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k N ) ) k N k N k N k N solver solver X ≡ { x 1 ,..., x n } Λ ≡ diag ( λ 1 ,..., λ n ) PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 6

Sequences of Eigenproblems Adjacent iteration cycles ITERATION ( ℓ ) ITERATION ( ℓ + 1 ) direct direct P ( ℓ ) ( X ( ℓ ) k 1 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 1 ) ) k 1 k 1 k 1 k 1 solver solver direct direct P ( ℓ ) ( X ( ℓ ) k 2 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 2 ) ) k 2 k 2 k 2 k 2 solver solver Next cycle direct direct P ( ℓ ) ( X ( ℓ ) k N , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k N ) ) k N k N k N k N solver solver X ≡ { x 1 ,..., x n } Λ ≡ diag ( λ 1 ,..., λ n ) PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 6

Sequences of Eigenproblems Adjacent iteration cycles ITERATION ( ℓ ) ITERATION ( ℓ + 1 ) direct direct P ( ℓ ) ( X ( ℓ ) k 1 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 1 ) ) k 1 k 1 k 1 k 1 solver solver direct direct P ( ℓ ) ( X ( ℓ ) k 2 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 2 ) ) k 2 k 2 k 2 k 2 solver solver Next cycle direct direct P ( ℓ ) ( X ( ℓ ) k N , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k N ) ) k N k N k N k N solver solver X ≡ { x 1 ,..., x n } Λ ≡ diag ( λ 1 ,..., λ n ) PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 6

Correlation between eigenproblems Definition and solving strategies Definition: Correlation Two adjacent problems P ( ℓ + 1 ) and P ( ℓ ) are said to be correlated when the eigenpairs ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) ) are in some relation with the eigenpairs ( X ( ℓ ) , Λ ( ℓ ) ) . PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 7

Correlation between eigenproblems Definition and solving strategies Definition: Correlation Two adjacent problems P ( ℓ + 1 ) and P ( ℓ ) are said to be correlated when the eigenpairs ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) ) are in some relation with the eigenpairs ( X ( ℓ ) , Λ ( ℓ ) ) . Alternative solving strategy Consider the entire sequence of eigenproblems P ( 1 ) ,..., P ( N ) as one self-contained problem; Exploit information provided from the correlation; In FLAPW, correlation appears in the way angles b/w eigenvectors of adjacent eigenproblems evolve [EDN, Blügel and Bientinesi 2012] Θ ( ℓ ) � 1 −� X ( ℓ − 1 ) X ( ℓ ) � , ˜ k i ≡ { θ 1 ,..., θ n } = diag k i � k i θ ( 2 ) � θ ( 3 ) � ··· � θ ( N ) θ ( 2 ) ≫ θ ( N ) for fixed k i : j j j j j PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 7

An alternative solving strategy Adjacent cycles ITERATION ( ℓ ) ITERATION ( ℓ + 1 ) iterative iterative P ( ℓ ) ( X ( ℓ ) k 1 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 1 ) ) k 1 k 1 k 1 k 1 solver solver iterative iterative P ( ℓ ) ( X ( ℓ ) k 2 , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k 2 ) ) k 2 k 2 k 2 k 2 solver solver Next cycle iterative iterative P ( ℓ ) ( X ( ℓ ) k N , Λ ( ℓ ) P ( ℓ + 1 ) ( X ( ℓ + 1 ) , Λ ( ℓ + 1 ) k N ) ) k N k N k N k N solver solver X ≡ { x 1 ,..., x n } Λ ≡ diag ( λ 1 ,..., λ n ) PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 8

Outline Sequences of correlated eigenproblems The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) ChFSI parallelization and numerical tests PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 9

Chebyshev Filtered Subspace Iteration method (ChFSI) Properties and algorithm evolution Iterative solver musts input: the full set of multiple starting vectors Z 0 ≡ X ( ℓ − 1 ) ( : , 1 : NEV ) ; k i needed: it can efficiently use dense linear algebra kernels (i.e. xGEMM ); needed: it avoids stalling when facing small clusters of eigenvalues; Chebyshev Subspace Iteration Firstly introduced in [Rutishauser 1969] A version (called CheFSI) tailored to electronic structure computation in [Zhou, Saad, Tiago and Chelikowski 2006] for sparse eigenvalue problems. Our ChFSI: 1) is tailored for dense eigenproblem sequences, 2) introduces a locking mechanism, 3) contains a refining inner loop, and 4) optimizes the polynomial degree. PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 10

ChFSI pseudocode I NPUT : Hamiltonian, approximate eigenvectors – Z 0 , extreme eigenvalues { λ 1 , λ NEV } , TOL , DEG . O UTPUT : NEV wanted eigenpairs ( Λ , W ) . PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 11

ChFSI pseudocode I NPUT : Hamiltonian, approximate eigenvectors – Z 0 , extreme eigenvalues { λ 1 , λ NEV } , TOL , DEG . O UTPUT : NEV wanted eigenpairs ( Λ , W ) . 1 Lanczos step. Identify the bounds for the eigenspectrum interval corresponding to the wanted eigenspace. PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 11

ChFSI pseudocode I NPUT : Hamiltonian, approximate eigenvectors – Z 0 , extreme eigenvalues { λ 1 , λ NEV } , TOL , DEG . O UTPUT : NEV wanted eigenpairs ( Λ , W ) . 1 Lanczos step. Identify the bounds for the eigenspectrum interval corresponding to the wanted eigenspace. R EPEAT U NTIL CONVERGENCE : 2 Chebyshev filter. Filter a block of vectors W ← − Z 0 . 3 Re-orthogonalize the vectors outputted by the filter; W = QR . 4 Compute the Rayleigh quotient G = Q † HQ . 5 Compute the primitive Ritz pairs ( Λ , Y ) by solving for GY = Y Λ . 6 Compute the approximate Ritz pairs ( Λ , W ← QY ) . 7 Check which one among the Ritz vectors converged . 8 Deflate and lock the converged vectors. E ND R EPEAT PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 11

The core of the algorithm: Chebyshev filter The basic principle Theorem Let | γ | > 1 and P m denote the set of polynomials of degree smaller or equal to m . Then the extremum p ∈ P m , p ( γ )= 1 max min t ∈ [ − 1 , 1 ] | p ( t ) | is reached by = C m ( t ) p m ( t ) . C m ( γ ) . where C m is the Chebyshev polynomial of the first kind of order m , defined as � cos ( m arccos ( t )) , t ∈ [ − 1 , 1 ] ; C m ( t ) = cosh ( m arccosh ( t )) , | t | > 1 . PMAA14, Lugano, Switzerland. July 4th E. Di Napoli and M. Berljafa Folie 12