A Chebyshev-based two-stage iterative method as an alternative to - PowerPoint PPT Presentation

A Chebyshev-based two-stage iterative method as an alternative to the direct solution of linear systems Mario Arioli m.arioli@rl.ac.uk CCLRC-Rutherford Appleton Laboratory with Daniel Ruiz (E.N.S.E.E.I.H.T) Cortona, 20-24 September, 2004 – p.1/30

Outline Problems with Preconditioning Techniques for improving the solvers The two phase algorithm Error Analysis Influence of the different parameters Numerical experiments Cortona, 20-24 September, 2004 – p.2/30

Problems with Preconditioning Au = b R N × N SPD A ∈ I Cortona, 20-24 September, 2004 – p.3/30

Problems with Preconditioning Au = b R N × N SPD A ∈ I M = R T R a SPD preconditioner R − T AR − 1 Ru = R − T b A = R − T AR − 1 , x = Ru , ˆ � b = R − T b . Cortona, 20-24 September, 2004 – p.3/30

Problems with Preconditioning Au = b R N × N SPD A ∈ I M = R T R a SPD preconditioner R − T AR − 1 Ru = R − T b A = R − T AR − 1 , x = Ru , ˆ � b = R − T b . Good clusterization but problem still ill-conditioned Cortona, 20-24 September, 2004 – p.3/30

A PDE example R 2 simply connected bounded polygonal. Let Ω ⊂ I ∀ u , v ∈ H 1 a ( u , v ) 0 (Ω) be a continuous and coercive bilinear form, and L ( v ) ∈ H − 1 (Ω) . The problem � Find u ∈ H 1 0 (Ω) such that a ( u , v ) = L ( v ) , ∀ v ∈ H 1 0 (Ω) , has a unique solution. Cortona, 20-24 September, 2004 – p.4/30

A PDE example � ∀ u , v ∈ H 1 a ( u , v ) = K ( x ) ∇ u · ∇ v d x , 0 (Ω) Ω � ∀ v ∈ H 1 L ( v ) = Ω 10 v d x , 0 (Ω) * Cortona, 20-24 September, 2004 – p.5/30

✁ � ✆ ✂ ✁ ✄ �✝ ✂✆ � ✂✆ ✁ ✁ A PDE example � ∀ u , v ∈ H 1 a ( u , v ) = K ( x ) ∇ u · ∇ v d x , 0 (Ω) Ω � ∀ v ∈ H 1 L ( v ) = Ω 10 v d x , 0 (Ω)  ✂☎✄ (1,1)  1 x ∈ Ω \ { Ω 1 ∪ Ω 2 } ,  10 6 K ( x ) = x ∈ Ω 1 , Ω Ω   1 10 4 x ∈ Ω 2 . (0,0) (0,1) Ω 2 ✂☎✄ * Cortona, 20-24 September, 2004 – p.5/30

A PDE example Using pdetool c � , we generated a mesh satisfying the usual regularity conditions of Ciarlet and we computed a finite-element approximation of the problem with the use of continuous piece-wise linear elements. The approximated problem is equivalent to the following system of linear equations: Au = b . In our mesh, the largest triangle has an area of 3 . 123 × 10 − 4 , therefore, the resulting linear system has 16256 triangles, 8289 nodes, and 7969 degrees of freedom. Cortona, 20-24 September, 2004 – p.6/30

A PDE example We used three kinds of preconditioners: the classical Jacobi diagonal matrix, M = diag ( A ) , the incomplete Cholesky decomposition of A with zero fill-in, and the incomplete Cholesky decomposition of A with drop tolerance 10 − 2 . Using the incomplete Cholesky decompositions, we computed the upper triangular matrix R such that M = R T R . Cortona, 20-24 September, 2004 – p.7/30

A PDE example We used three kinds of preconditioners: the classical Jacobi diagonal matrix, M = diag ( A ) , the incomplete Cholesky decomposition of A with zero fill-in, and the incomplete Cholesky decomposition of A with drop tolerance 10 − 2 . Using the incomplete Cholesky decompositions, we computed the upper triangular matrix R such that M = R T R . κ ( M − 1 A ) λ min λ max M 2.6 10 9 3.7 10 − 3 9.6 10 6 I 6.8 10 8 3.1 10 − 9 Jacobi 2.08 9.4 10 7 1.7 10 − 8 Inc. Cholesky(0) 1.6 Inc. Cholesky(10 − 2 ) 6.2 10 6 1.8 10 − 7 1.1 Cortona, 20-24 September, 2004 – p.7/30

A PDE example Preconditioner M Inc. Cholesky(10 − 2 ) µ Jacobi Inc. Cholesky(0) λ max / 10 3 3 λ max / 500 5 λ max / 200 18 λ max / 100 43 3 λ max / 50 11 λ max / 20 32 λ max / 10 > 200 68 3 λ max / 5 9 λ max / 2 40 Number of eigenvalues in [ λ min , µ ] . Cortona, 20-24 September, 2004 – p.8/30

Problems with Preconditioning Good clusterization but problem still ill-conditioned Cortona, 20-24 September, 2004 – p.9/30

Problems with Preconditioning Good clusterization but problem still ill-conditioned Several right-hand sides in succession Cortona, 20-24 September, 2004 – p.9/30

Problems with Preconditioning Good clusterization but problem still ill-conditioned Several right-hand sides in succession Implicit matrix Cortona, 20-24 September, 2004 – p.9/30

Techniques for improving the solvers Cortona, 20-24 September, 2004 – p.10/30

Techniques for improving the solvers Deflation: compute the invariant space corresponding to the smallest eigenvalues Cortona, 20-24 September, 2004 – p.10/30

Techniques for improving the solvers Deflation: compute the invariant space corresponding to the smallest eigenvalues Filtering + Lanczos Cortona, 20-24 September, 2004 – p.10/30

The two-phase iterative approach Fix µ such that λ min ( � A ) < µ < λ max ( � A ) Cortona, 20-24 September, 2004 – p.11/30

The two-phase iterative approach Fix µ such that λ min ( � A ) < µ < λ max ( � A ) Start with an initial set of s randomly generated vectors ( s is the block-size), and “ damp ”, in these starting vectors, the eigenfrequencies associated with all the eigenvalues in [ µ, λ max ( � A )] . Cortona, 20-24 September, 2004 – p.11/30

The two-phase iterative approach Fix µ such that λ min ( � A ) < µ < λ max ( � A ) Start with an initial set of s randomly generated vectors ( s is the block-size), and “ damp ”, in these starting vectors, the eigenfrequencies associated with all the eigenvalues in [ µ, λ max ( � A )] . Compute all the eigenvectors associated with all the eigenvalues in the range [ λ min ( � A ) , µ ] . Cortona, 20-24 September, 2004 – p.11/30

The two-phase iterative approach Fix µ such that λ min ( � A ) < µ < λ max ( � A ) Start with an initial set of s randomly generated vectors ( s is the block-size), and “ damp ”, in these starting vectors, the eigenfrequencies associated with all the eigenvalues in [ µ, λ max ( � A )] . Compute all the eigenvectors associated with all the eigenvalues in the range [ λ min ( � A ) , µ ] . If the eigenvalues are well clustered the number of remaining eigen- values in [ λ min ( � A ) , µ ] , with reasonable µ ( λ max / 100 , or λ max / 10 ), should be small compared to the size of the linear system. Cortona, 20-24 September, 2004 – p.11/30

The two-phase iterative approach 2 1.8 26 eigs. are below the threshold ( µ = λ max /10) ← 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 −12 −10 −8 −6 −4 −2 0 10 10 10 10 10 10 10 Smallest eig. = 1.0897e−13; Largest eig. = 2.5923 Cortona, 20-24 September, 2004 – p.12/30

The two-phase iterative approach We can use Chebyshev polynomials to damp � T 0 ( ω ) = 1 , T 1 ( ω ) = ω n ≥ 1 . . T n +1 ( ω ) = 2 ωT n ( ω ) − T n − 1 ( ω ) If we consider d > 1 and H n ( ω ) = T n ( ω ) T n ( d ) , then H n has minimum l ∞ norm on the interval [ − 1 , 1] over all polynomials Q n of degree less than of equal to n and satisfying the condition Q n ( d ) = 1 , and we have 1 ω ∈ [ − 1 , 1] H n ( ω ) = max T n ( d ) . Cortona, 20-24 September, 2004 – p.13/30

The two-phase iterative approach → ω µ ( λ ) = λ max ( � A ) + µ − 2 λ λ ∈ R �− , λ max ( � A ) − µ with λ min ( � A ) < µ < λ max ( � A ) given above. This transformation maps λ max ( � A ) to − 1 , µ to 1 , and 0 to ω µ (0) = d µ = λ max ( � A ) + µ > 1 . λ max ( � A ) − µ Then, P n ( λ ) = T n ( ω µ ( λ )) , T n ( d µ ) Cortona, 20-24 September, 2004 – p.14/30

The two-phase iterative approach � A = UΛU T the eigendecomposition of the SPD matrix � A , with Λ = diag( λ i ) Cortona, 20-24 September, 2004 – p.15/30

The two-phase iterative approach � A = UΛU T the eigendecomposition of the SPD matrix � A , with Λ = diag( λ i ) n � z = u i ξ i , i =1 with ξ i = u T i z , i = 1 , . . . , n , Cortona, 20-24 September, 2004 – p.15/30

The two-phase iterative approach � A = UΛU T the eigendecomposition of the SPD matrix � A , with Λ = diag( λ i ) n � z = u i ξ i , i =1 with ξ i = u T i z , i = 1 , . . . , n , n � v = P n ( � A ) z = u i ( P n ( λ i ) ξ i ) , i =1 Cortona, 20-24 September, 2004 – p.15/30

The two-phase iterative approach � A = UΛU T the eigendecomposition of the SPD matrix � A , with Λ = diag( λ i ) n � z = u i ξ i , i =1 with ξ i = u T i z , i = 1 , . . . , n , n � v = P n ( � A ) z = u i ( P n ( λ i ) ξ i ) , i =1 The eigencomponents of the resulting vector v are close to that of the initial vector z for all i such that λ i is close to 0 and relatively much smaller for large enough degree n and all i such that λ i ∈ [ µ, λ max ( � A )] . Cortona, 20-24 September, 2004 – p.15/30