Approximate eigen-decomposition preconditioners for solving numerical PDE problems

https://doi.org/10.1016/j.amc.2004.11.018Get rights and content

Abstract

Approximate matrix inverse has played key role in constructing preconditioners since last 1990s, because of their solid theoretical background. Sparse approximate inverse preconditioners (SAIPs) have attracted much attention recently, because their potential usefulness in a parallel environment. As a general matrix approach, in this paper, we propose an approximate eigen-decomposition preconditioners by combining a FFT-like multiplication algorithm. Some numerical tests are given to show this algorithm is more effective than the traditional method such as ILU with PETSc for solving a wide class of discrete elliptic problems.

Introduction

Since the past more than 10 years, we have been involved in parallel solving for oil-reservoir simulation system. This work stimulate us to do a deep research on parallel preconditioning technique for implicit iteration schemes arisen from numerical PDE. In this paper we propose an approximate eigen-decomposition preconditioning algorithm to get more efficient algorithms.

Let us consider a general linear systemAu=f,where sparse matrix A arises from numerical PDE and usually is assumed to be positive, but not necessary to be symmetry in general. To solve (1.1) efficiently, in the most current commercial software, one adopts iteration methods if the order of the system exceeds 3000–5000. There are many iterative methods, see [3]. Usually CG-like method and GMRES-like methods are chosen for symmetry and non-symmetry case, respectively. The convergent rate depends mainly on the eigenvalue distribution of the matrix A. With more closer and clustering of all eigenvalues, convergence rate is more fast. Since the eigenvalue distribution of (1.1) is getting worse as the step of space or time becomes smaller, it is necessary to find a so-called preconditioning algorithm to improve the condition number, i.e. to reform the original linear system toBAu=Bf,orAB(B-1u)=f,where B is called a left or right preconditioner which can be taken as an approximation inverse of A in (1.1). The idea of approximate matrix inverse has played key role in constructing preconditioners since last 1990s [3], because of their solid theoretical background. Recently, as a general and simple approach, matrix approach sparse approximate inverse preconditioners (SAIPs) [1] have attracted much attention, because their potential usefulness in a parallel environment. However, it is worth to point that because (1.1) is a partial difference equation system arisen from a PDE problem and quite different from a general sparse matrix system. Based on some special PDE properties we have proposed a number of efficient preconditioning algorithms, such as approximate Green function preconditioning [4], presymmetry iteration algorithms [5], diagonal-block preconditioning for enhancing the elliptic dominant part [7], etc. In this paper, we propose an approximate eigen-decomposition preconditioners by combining a FFT-like multiplication algorithm. Though the computational complexity for each matrix multiplication in our case is O(N Log N) instead of O(N) in SAIPs. However, once approximate eigen-decomposition can be applied, e.g. in some PDE with special partitions in 2-D and 3-D, the iteration rate will be greatly reduced.

Some examples are given for solving discrete Helmholtz equation on equilateral triangle with zero boundary condition and on 2-D hexagon and 3-D dodecahedron domains with periodic boundary conditions. Numerical tests are listed to show this algorithm is more efficient, comparing with the usual ILU, e.g. in PETSc (Portable, Extensible Toolkit for Scientific Computation) [2].

Section snippets

Approximate matrix eigen-decomposition preconditioning

As is well known, an approximate sparse inverse may be a good preconditioning [1]. Now we propose that a reasonable approximate eigen-decomposition, based on fast algorithms, also can be taken as a preconditioner. In this case, it need not require the preconditioner B to be sparse, but the working amount of Bu must less O(N2), e.g. O(N Log N), where N is the matrix order.

Suppose A0 is an approximation of A. Let E=I-A0-1A, if ρ(E) < 1, thenA-1={I-E}-1A0-1=I+∑k=1∞EkA0-1.We may define various levels

Approximate eigen-decomposition for solving 2-D Helmholtz equation over a triangle and a hexagon domain

Consider the following 2-D Helmholtz equation-Δu+qu=finΩwithu|∂Ω=0oru|∂Ω=periodicfor a triangle or a parallel hexagon domain, respectively.

First we replace the usual 2-D Cartesian coordinates by the following three direction coordinates. Let e1 and e2 be two independent vectors of the given domain, a bi-orthogonal vector set n1 and n2 is defined asn1=(e1,e1)e2-(e1,e2)e1(e1,e1)(e2-e2)-(e1,e2)2,n2=(e1,e1)e1-(e1,e2)e2(e1,e1)(e2-e2)-(e1,e2)2.Thus, for any 2-D point P, we define its three direction

3-D dodecahedron partition case

In 3-D case we consider parallel dodecahedron domain with six direction partition. For a given three linear independent vectors: e1, e2, e3, e4 = −(e1 + e2 + e3), there are an bi-orthogonal vectors n1, n2, n3, such that(ej,nk)=δj,k,(j,k=1,2,3).Setn4=n1-n2,n5=n2-n3,n6=n3-n1,and 6-D affine coordinates for a 3-D point P are defined ast1=(P,n1),t2=(P,n2),t3=(P,n3),t4=(P,n4),t5=(P,n5),t6=(P,n6),where t4 = t1 − t2, t5 = t2 − t3, t6 = t3 − t1.

A basic 3-D domain Ω is a parallel dodecahedron with 14 vertex and 24 edges defined

Generalized discrete Fourier transformation (HDFT)

To explain our idea more clear, here we describe the fast algorithm in 2-D case only. Let 1 − N ⩽ k1 ⩽ N, −N ⩽ k2, k3 ⩽ N − 1, k1 + k2 + k3 = 0, and 3N2 given number (real or complex) fk1,k2,k3. DefineFn1,n2,n3=∑k1=1-NN∑k2,k3=-Nk2+k3=-k1N-1fk1,k2,k3e2πi3(k1n1+k2n2+k3n3)/Nas the generalized discrete Fourier transformation (HDFT).

The transform and the related inverse transform can be rewritten in terms of matrix formF=GNf,f=(GN)-1F,where matrix GN is a complex dense matrix of order 3N2.

It is obvious that the

Numerical tests

Test 1

Fast Laplace solver on a unit regular hexagon domain with three direction periodic boundary conditions.

Table 1, Table 2, Table 3, Table 4 list three term comparisons on iteration counts, Mflops and CPU time among four different PCG solver: CG, LU, the best ILU(k) running on PETSc and our HFT algorithm, where mesh side h=1256. Fig. 9, Fig. 10, Fig. 11, Fig. 12, Fig. 13 show CPU time comparison among the above four algorithms for four refining mesh step size 1h=32,64,128,256.

Test 2

Fast Helmholtz solver

Acknowledgement

The figures and part computing are done by Yao Jifeng and Xie Ligang.

References (9)

  • L.Yu. Kolotilina et al.

    Factorized sparse approximate inverse preconditionings I: theory

    SIAM Journal on Matrix Analysis and Applications

    (1993)
  • PETc home page, http://www-unix.mcs.anl.gov/petsd/petsc-2/,...
  • Y. Saad, ILUs and Factorized Approximate Inverses are Strongly, Bollhöfer,...
  • J. Sun

    A class of local green-like parallel preconditioner algorithms for elliptic discrete equations: I

    Mathetitica Numerica Sinica

    (1995)
There are more references available in the full text version of this article.

Cited by (3)

1

Project supported by the Major Basic Project of China (No.G19990328) and National Natural Science Foundation of China (No. 60173021).

View full text
Copyright © 2004 Elsevier Inc. All rights reserved.