An Algorithm for Solving Non-Symmetric Saddle-Point Systems - PowerPoint PPT Presentation

An Algorithm for Solving Non-Symmetric Saddle-Point Systems Jaroslav Haslinger, Charles University, Prague s Kozubek, Vˇ Tom´ aˇ SB–TU Ostrava cera, Vˇ Radek Kuˇ SB–TU Ostrava HARRACHOV, August 2007 • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

OUTLINE Motivation: Fictitious domain method Algorithm PSCM: Schur complement method + Null-space method Inner solver: Projected BiCGSTAB Preconditioning: Hierarchical multigrid Singular matrices: Poisson-like solver based on circulants Numerical experiments • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

OUTLINE Motivation: Fictitious domain method Algorithm PSCM: Schur complement method + Null-space method Inner solver: Projected BiCGSTAB Preconditioning: Hierarchical multigrid Singular matrices: Poisson-like solver based on circulants Numerical experiments • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

MODEL PROBLEM 1: Dirichlet problem − ∆ u = f on ω (1) u = g in γ ≡ ∂ω (2) Fictitious domain method (FDM): PDE (1) is solved on the fictitious domain Ω , ω ⊂ Ω , with a simple geome- try. The corresponding stiffness matrix A is structured. The original boundary conditions (2) on γ are enforced by Lagrange multipliers or control variables. Ω Γ γ δ ω Ξ Classical FDM with Γ = γ Smooth FDM with Γ � = γ � � � � � � � � � � � � B ⊤ B ⊤ A u f A u f γ Γ = = B γ 0 λ γ g B γ 0 λ Γ g • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

MODEL PROBLEM 2: Signorini problem − ∆ u = f on ω (3) ∂u ≥ 0 , ( u − g ) ∂u u − g ≥ 0 , in γ ≡ ∂ω = 0 (4) ∂n γ ∂n γ FDM formulation uses the non-differentiable max -function to express BC (4): � Au + B Γ λ Γ = f (5) C γ,i u = max { 0 , C γ,i u − ρ ( B γ,i u − g i ) } , i = 1 , . . . , m where B γ,i , B Γ ,i and C γ,i are rows of Dirichlet and Neumann trace matrices, respectively. The equations (5) can be solved by the semi-smooth Newton method, in which � � B ⊤ A Γ Jacobian = ∂G ( u ) 0 is determined by the generalized derivative ∂G ( u ) . • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

MODEL PROBLEM 2: Newton method = Active set algorithm (0) Set k := 1 , ρ > 0 , ε u > 0 , u (0) ∈ R n , λ (0) ∈ R m . (1) Define the inactive and active sets by: I k := { i : C γ,i u k − 1 − ρ ( B γ,i u k − 1 − g i ) ≤ 0 } A k := { i : C γ,i u k − 1 − ρ ( B γ,i u k − 1 − g i ) > 0 } (2) Solve:     � � B ⊤ f A u k Γ     = B γ, A k 0 g A k λ k Γ C γ, I k 0 0 (3) If � u k − u k − 1 � / � u k � ≤ ε u , return u := u k . (4) Set k := k + 1 , and go to step (1). Remark: The mixed Dirichlet-Neumann problem is solved in each Newton step, that is described by the non-symmetric saddle-point system. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

FORMULATION: Non-symmetric sadle-point system � � � � � � A B ⊤ u f 1 = B 2 0 λ g General assumptions A . . . non-symmetric ( n × n ) –matrix . . . singular with p = dim Ker A B 1 , B 2 . . . full rank ( m × n ) –matrices . . . B 1 � = B 2 Special FDM assumptions • n is large ( n = 4198401 ) • m ≪ n ( m = 360 ) • p ≪ m ( p = 1 ) • A is structured so that actions of A † or ( A − 1 ) are ”cheap” • B 1 , B 2 are highly sparse so that their actions are ”cheap” • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

OUTLINE Motivation: Fictitious domain method Algorithm PSCM: Schur complement method + Null-space method Inner solver: Projected BiCGSTAB Preconditioning: Hierarchical multigrid Singular matrices: Poisson-like solver based on circulants Numerical experiments • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

ALGORITHMS based on the Schur complement reduction Case 1: A non-singular, symmetric Case 2: A non-singular, non-symmetric Case 3: A singular, symmetric Case 4: A singular, non-symmetric • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Case 1: A non-singular, symmetric � A B ⊤ � � u � � f � u = A − 1 ( f − B ⊤ λ ) = ⇒ = BA − 1 B ⊤ � λ = BA − 1 f − g ⇒ = λ g B 0 � �� negative Schur complement S Algorithm Assemble d := BA − 1 f − g . 1 ◦ 2 ◦ Solve iteratively S λ = d with S := BA − 1 B ⊤ . 3 ◦ Assemble u := A − 1 ( f − B ⊤ λ ) . If A is positive defined, then CGM can be used. Matrix-vector products S µ are performed by: � � A − 1 � ��� B ⊤ µ S µ := B • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Case 2: A non-singular, non-symmetric � A B ⊤ � � u � � f � u = A − 1 ( f − B ⊤ = ⇒ 1 λ ) 1 = B 2 A − 1 B ⊤ � λ = B 2 A − 1 f − g = ⇒ λ g B 2 0 � �� 1 negative Schur complement S Algorithm is analogous. • an iterative method for non-symmetric matrices is required (GMRES, BiCG, BiCGSTAB, ...) � A B ⊤ � � � � A B ⊤ � I 0 1 1 A := = B 2 A − 1 I 0 − S B 2 0 Theorem 1 Let A be non-singular. Then A is invertible iff S is invertible. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Case 3: A singular, symmetric • a generalized inverse A † satisfying A = AA † A • an ( n × p ) –matrix N whose columns span Ker A Au + B ⊤ λ = f f − B ⊤ λ ∈ Im A ⊥ Ker A ⇐ ⇒ � � u = A † ( f − B ⊤ λ ) + N α N ⊤ ( f − B ⊤ λ ) = 0 & The reduced system: � � � � � � BA † B ⊤ BA † f − g − BN λ Bu = g = − N ⊤ B ⊤ − N ⊤ f α 0 ⇓ BA † B ⊤ λ − BN α = BA † f − g If A is positive semidefinite, then it corresponds to the algebra in FETI DDM. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Case 4: A singular, non-symmetric • a generalized inverse A † • columns of ( n × p ) –matrices N , M span Ker A , Ker A ⊤ , respectively Au + B ⊤ f − B ⊤ 1 λ ∈ Im A ⊥ Ker A ⊤ ⇐ ⇒ 1 λ = f � � u = A † ( f − B ⊤ M ⊤ ( f − B ⊤ 1 λ ) + N α 1 λ ) = 0 & The reduced system: � � � � � � B 2 A † B ⊤ B 2 A † f − g − B 2 N λ 1 B 2 u = g = − M ⊤ B ⊤ − M ⊤ f α 0 1 ⇓ B 2 A † B ⊤ 1 λ − B 2 N α = B 2 A † f − g • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

� A B ⊤ � 1 Theorem 2 The saddle-point matrix A := is invertible iff B 2 0   B 1 has full row-rank   Ker A ∩ Ker B 2 = { 0 } (NSC)    A Ker B 2 ∩ Im B ⊤ 1 = { 0 } Remark: The third equality is equivalent to the MinMax condition that is well- known in the continuous setting: v ⊤ Au ∃ C > 0 : min max � v �� u � ≥ C u ∈ Ker B 2 , u � =0 v ∈ Ker B 1 , v � =0 • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

The generalized Schur complement: the matrix of the reduced system � − B 2 A † B ⊤ � B 2 N 1 S := M ⊤ B ⊤ 0 1 Theorem 3 The following three statements are equivalent: • The necessary and sufficient condition (NSC) holds. • A is invertible. • S is invertible. Remark: The generalized Schur complement S is not defined uniquely. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

First step of the algorithm = Schur complement reduction:  � � � � � �  G ⊤ F λ d � � � � � �   1 = A B ⊤ u f G 2 0 α e 1 ⇐ ⇒ = B 2 0 λ g    u = A † ( f − B ⊤ 1 λ ) + N α How to solve the reduced system again with the saddle-point structure? � � A † � ��� B ⊤ • matrix-vector products via F µ := B 2 1 µ • G 1 , G 2 , d , e may be assembled   1) Again the Schur complement reduction (the second elimination)     E α = r with E := G 2 F − 1 G ⊤ 1 , then λ := F − 1 ( d − G ⊤  1 α ) and u  (R.K., Appl. Math. 50(2005))     2) Null-space method    (Farhat, Mandel, Roux: FETI DDM, 1994) • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Second step of the algorithm = Null-space method: � � � � � � G ⊤ λ d F 1 = α e G 2 0 Two orthogonal projectors P 1 and P 2 onto Ker G 1 and Ker G 2 : P k : R m �→ Ker G k , P k := I − G ⊤ k ( G k G ⊤ k ) − 1 G k , k = 1 , 2 Ker P k = Im G ⊤ P k G ⊤ ⇐ ⇒ k = 0 Property: k P 1 F λ + P 1 G ⊤ • P 1 splits the saddle-point structure: 1 α = P 1 d α := ( G 1 G ⊤ 1 ) − 1 ( G 1 d − G 1 F λ ) P 1 F λ = P 1 d , G 2 λ = e , λ Im ∈ Im G ⊤ • P 2 decomposes λ = λ Im + λ Ker , λ Ker ∈ Ker G 2 2 , λ Im := G ⊤ 2 ( G 2 G ⊤ 2 ) − 1 e G 2 λ = G 2 λ Im = e = ⇒ At first: P 1 F λ Ker = P 1 ( d − F λ Im ) on Ker G 2 At second: • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit