Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Solving Nonlinear Eigenvalue Problems with SLEPc Jose E. Roman D. Sistemes Inform` atics i Computaci´ o Universitat Polit` ecnica de Val` encia, Spain Joint work with Carmen Campos NLA-HPC 2013, Hsinchu (Taiwan) 1/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Nonlinear Eigenproblems Increasing interest in nonlinear eigenvalue problems arising in many application domains ◮ Structural analysis with damping effects ◮ Vibro-acoustics (fluid-structure interaction) ◮ Linear stability of fluid flows Problem types ◮ QEP: quadratic eigenproblem, ( λ 2 M + λC + K ) x = 0 ◮ PEP: polynomial eigenproblem, P ( λ ) x = 0 ◮ REP: rational eigenproblem, P ( λ ) Q ( λ ) − 1 x = 0 ◮ NEP: general nonlinear eigenproblem, T ( λ ) x = 0 Test cases available in the NLEVP collection [Betcke et al. 2013] 2/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Example from NLEVP Collection Connected damped mass-spring system Second-order differential equation: M d 2 dt 2 x + C d dtx + Kx = 0 3/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Solution Strategy We are interested in large-scale nonlinear eigenvalue problems, with methods appropriate for sparse matrices on parallel computers ◮ Only a small part of the spectrum Alternatives for the PEP: 1. Projection method such as Jacobi-Davidson, P ( θ ) u ⊥ K 2. Linearization: solve via a linear eigenproblem Linearization can leverage existing linear eigensolvers, but... ◮ Dimension of linearized problem is dn → Need memory-efficient variants capable of exploiting structure The NEP requires completely different approaches 4/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers SLEPc : Scalable Library for Eigenvalue Problem Computations A general library for solving large-scale sparse eigenproblems on parallel computers ◮ For standard and generalized eigenproblems ◮ For real and complex arithmetic ◮ For Hermitian or non-Hermitian problems ◮ Also support for SVD, QEP, and more ( λ 2 M + λC + K ) x = 0 Ax = λx Ax = λBx Av i = σ i u i Co-authors: C. Campos, E. Romero, A. Tomas http://www.grycap.upv.es/slepc Current version: 3.4 (released July 2013) 5/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers PETSc/SLEPc Numerical Components PETSc Nonlinear Systems Time Steppers Trust Pseudo Line Backward Other Euler Other Search Region Euler Time Step Krylov Subspace Methods GMRES CG CGS Bi-CGStab TFQMR Richardson Chebychev Other Preconditioners Additive Block Jacobi ILU ICC LU Other Schwarz Jacobi Matrices Compressed Block Symmetric Dense CUSP Other Sparse Row CSR Block CSR Vectors Index Sets Standard CUSP Indices Block Stride Other 6/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers PETSc/SLEPc Numerical Components PETSc SLEPc Nonlinear Systems Time Steppers Quadratic Eigensolver Nonlinear Eigensolver Trust Pseudo Line Backward Linear- Other Euler Other ization Q-ArnoldiQ-Lanczos SLP RII N-Arnoldi Search Region Euler Time Step Krylov Subspace Methods SVD Solver M. Function Cross Matrix Lanczos Thick R. Cyclic GMRES CG CGS Bi-CGStab TFQMR Richardson Chebychev Other Krylov Product Lanczos Linear Eigensolver Preconditioners Additive Block Jacobi ILU ICC LU Other Krylov-Schur GD JD RQCG CISS Other Schwarz Jacobi Matrices Spectral Transformation Compressed Block Symmetric Shift Shift-and-invert Cayley Dense CUSP Other Fold Preconditioner Sparse Row CSR Block CSR Vectors Index Sets IP DS FN Standard CUSP Indices Block Stride Other 6/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Outline Introduction 1 PEP: Polynomial Eigensolvers 2 Q-Arnoldi and TOAR for QEP Extension to PEP Usage and Numerical Results NEP: General Nonlinear Eigensolvers 3 User Interface Examples 7/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers PEP: Polynomial Eigensolvers 8/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Quadratic Arnoldi ( λ 2 M + λC + K ) x = 0 QEP: linearized as � N � � � 0 N 0 A − λB = − λ − K − C 0 M with either N = I or N = − K . The eigenvector is [ x ∗ , λx ∗ ] ∗ Q-Arnoldi [Meerbergen 2008] builds an Arnoldi relation for B − 1 A � � V 0 � V 0 v 0 � � � 0 I k k = H − M − 1 K − M − 1 C V 1 V 1 v 1 k k V 1 k = [ V 0 k , v 0 ] H From the 1st row block: Hence, only { V 0 k , H, v } needs to be stored 9/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Q-Arnoldi For each i = 1 , . . . , k 1. Expand: w 0 = v 1 , w 1 = − M − 1 ( Kv 0 + Cv 1 ) 2. Gram-Schmidt coefficients � V 0 � V 0 � ∗ ∗ w 0 + V 1 ∗ w 1 v 0 � h j = k w = k k V 1 v 1 u ∗ w k ∗ w 1 = H ∗ [ V 0 with V 1 k , v 0 ] ∗ w 1 k 3. Gram-Schmidt update w 0 − [ V 0 v 0 k , v 0 ] h j ˜ = w 1 − v 1 [ V 0 k , v 0 ] H, v 1 � � ˜ = h j 4. Normalize with � ˜ v � 2 10/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Stabilized Variant Q-Arnoldi presents instability when � H � gets large ◮ Problem comes from representing the Krylov basis with V 0 k , v 0 � � matrices and H , whose columns are not orthogonal TOAR: Two-level Orthogonalization Arnoldi [Su et al. 2008] ◮ Main idea: use only orthogonal matrices � � G 0 � U k +1 � V k = k G 1 U k +1 k Arnoldi relation: � � U k +1 G 0 � U k +2 G 0 � 0 I � � k k +1 = H k − M − 1 K − M − 1 C U k +1 G 1 U k +2 G 1 k +1 k 11/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Main Steps in TOAR At step j , the new vector w is computed similarly to Q-Arnoldi. Then the bottom part w 1 is used to extend U j +1 to obtain a basis of span([ V 0 j , V 1 j , w 0 , w 1 ]) The new basis vector u j +2 ∈ U ⊥ j +1 is obtained via Gram-Schmidt, � � � [ U j +1 u j +2 ] g 0 � g 0 together with g = so that w = g 1 [ U j +1 u j +2 ] g 1 Gram-Schmidt coefficients: � ∗ � U ∗ � G 0 U j +1 g 1 � g 1 � � � � h j = V ∗ j j +1 = G ∗ j j j w = G 1 U ∗ j U j +1 ˆ w + u j +2 α w ˆ j j +1 The orthogonalization of w can be expressed as � �� g 0 � G 0 �� � � � � U j +1 u j +2 j w = w − V j h j = ˜ − h j G 1 � � g 1 U j +1 u j +2 j 12/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Linearizations P ( λ ) x = 0 PEP: P ( λ ) = A 0 f 0 ( λ ) + · · · + A d f d ( λ ) Monomial basis: f j ( λ ) = λ j ◮ Companion linearization L 0 − λ L 1 , with I I x ... ... λx L 0 = L 1 = y = . . I I . λ d − 1 x − A 0 − A 1 · · · − A d − 1 A d Chebyshev basis: f j ( λ ) = τ j ( λ ) with τ j ( λ ) = 2 λτ j − 1 ( λ ) − τ j − 2 ( λ ) ◮ Linearization contains more nonzero blocks Other bases/lineariz. [Amiraslani et al. 2009, De Ter´ an et al. 2010] 13/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Shift-and-Invert Transformation For computing eigenvalues closest to a target σ ◮ Can be applied to the linearization or the original problem ( θ 2 M σ + θC σ + K σ ) x = 0 Transformed QEP: M σ = σ 2 M + σC + K C σ = C + 2 σM K σ = M The relation with the original eigenvalue is θ = ( λ − σ ) − 1 In the polynomial case, use the reverse of the shifted polynomial ˜ P ( θ ) with coefficients d − k � j + k � � σ j A j + k T k = k j =0 14/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Polynomial TOAR Arnoldi decomposition of order j for S := L − 1 1 L 0 v 0 . � � . SV j = V j v H j , vectors split as v = . v d − 1 Due to the structure of S we have the relations V i +1 � V i v i � = H j i = 0 , . . . , d − 2 j j If U j + d is an orthogonal basis of span([ V 0 j , . . . , V d − 1 , v 0 , . . . , v d − 1 ]) = span([ V 0 j , v 0 , . . . , v d − 1 ]) j then it is possible to represent � V i v i � � G i g i � = U j + d , i = 0 , . . . , d − 1 j j � � Arnoldi: S ( I d ⊗ U j + d − 1 ) G j = ( I d ⊗ U j + d ) G j g j +1 H j Can also be adapted to Chebyshev basis [Kressner/Roman 2013] 15/30
Introduction PEP: Polynomial Eigensolvers NEP: General Nonlinear Eigensolvers Scaling Scaling the QEP can improve the backward error of the solutions obtained via linearization [Tisseur 2000] The FLV scheme [Fan et al. 2004] replaces the QEP with the scaled problem ( µ 2 γ 2 δM + µγδC + δK ) x = 0 � where λ = γµ with γ = � K � / � M � , δ = 2 / ( � K � + γ � C � ) Betcke [2008] extends this to polynomials of arbitrary degree as well as more general diagonal scalings D 1 P ( λ ) D 2 ◮ Scaling is optimal for a given λ 16/30
Recommend
More recommend
