The approximate Dirichlet Domain Decomposition method. Part II: Applications to 2nd-order Elliptic B.V.P.s

In the first part of this article series, we had derived Domain Decomposition (DD) preconditioners containing three block matrices which must be specified for specific applications. In the present paper, we consider finite element equations arising from the DD discretization of plane, symmetric, 2nd-order, elliptic b.v.p.s and specify the matrices involved in the preconditioner via multigrid and hierarchical techniques. The resulting DD-PCCG methods are asymptotically almost optimal with respect to the operation count and well suited for parallel computations on MIMD computers with local memory and message passing. The numerical experiments performed on a transputer hypercube confirm the efficiency of the DD preconditioners proposed.     

Domain decomposition preconditioners for the conjugate gradient method

The aim of this paper is to compare different techniques for constructing incomplete Domain Decomposition preconditioners which can be used on parallel computers with a large number of processors (typically from 1 to 100) in connection with the Conjugate Gradient method for solving symmetric linear systems that arise from finite difference discretization of two dimensional elliptic partial differential equations. We summarize some methods which have been presented during the last years and we introduce new ideas to improve the degree of parallelism and the efficiency of our preconditioners. Finally, we show some numerical results that demonstrate the usefulness of the preconditioners described in this paper.     

New preconditioners for systems of linear equations with Toeplitz structure

In this paper, we consider applying the preconditioned conjugate gradient (PCG) method to solve system of linear equations $T x = \mathbf b $ where $T$ is a block Toeplitz matrix with Toeplitz blocks (BTTB). We first consider Level-2 circulant preconditioners based on generalized Jackson kernels. Then, BTTB preconditioners based on a splitting of BTTB matrices are proposed. We show that the BTTB preconditioners based on splitting are special cases of embedding-based BTTB preconditioners, which are also good BTTB preconditioners. As an application, we apply the proposed preconditioners to solve BTTB least squares problems. Our preconditioners work for BTTB systems with nonnegative generating functions. The implementations of the construction of the preconditioners and the relevant matrix-vector multiplications are also presented. Finally, Numerical examples, including image restoration problems, are presented to demonstrate the efficiency of our proposed preconditioners.     

Preconditioning of the p-version of the finite element method

The p-version finite element method for linear, second-order elliptic equations in an arbitrary, sufficiently smooth (incl. polygonal), bounded domain is studied in the framework of the Domain Decomposition (DD) method. Two types of square reference elements are used with coordinate functions given by the products of the integrated Legendre polynomials. Estimates for the condition numbers and some useful inequalities are given. We consider preconditioning of the problems arising on subdomains and of the Schur complement, as well as the derivation and analysis of the DD preconditioner for the entire system. This is done for a class of curvilinear finite elements. We obtain several DD preconditioners for which the generalized condition numbers vary from ((log p)³) to (1).

This paper is based on [19–21,27]. We have omitted most of the proofs in order to shorten it and have described instead what could be done as well as outlined some additional ideas. The full proofs omitted can in most cases be found in [19,20,27].     

Efficient Nonlinear Solvers for Nodal High-Order Finite Elements in 3D

Conventional high-order finite element methods are rarely used for industrial problems because the Jacobian rapidly loses sparsity as the order is increased, leading to unaffordable solve times and memory requirements. This effect typically limits order to at most quadratic, despite the favorable accuracy and stability properties offered by quadratic and higher order discretizations. We present a method in which the action of the Jacobian is applied matrix-free exploiting a tensor product basis on hexahedral elements, while much sparser matrices based on Q ₁ sub-elements on the nodes of the high-order basis are assembled for preconditioning. With this “dual-order” scheme, storage is independent of spectral order and a natural taping scheme is available to update a full-accuracy matrix-free Jacobian during residual evaluation. Matrix-free Jacobian application circumvents the memory bandwidth bottleneck typical of sparse matrix operations, providing several times greater floating point performance and better use of multiple cores with shared memory bus. Computational results for the p-Laplacian and Stokes problem, using block preconditioners and AMG, demonstrate mesh-independent convergence rates and weak (bounded) dependence on order, even for highly deformed meshes and nonlinear systems with several orders of magnitude dynamic range in coefficients. For spectral orders around 5, the dual-order scheme requires half the memory and similar time to assembled quadratic (Q ₂) elements, making it very affordable for general use.     

On the robustness and optimality of algebraic multilevel methods for reaction–diffusion type problems

This paper is on preconditioners for reaction–diffusion problems that are both, uniform with respect to the reaction–diffusion coefficients, and optimal in terms of computational complexity. The considered preconditioners belong to the class of so-called algebraic multilevel iteration (AMLI) methods, which are based on a multilevel block factorization and polynomial stabilization. The main focus of this work is on the construction and on the analysis of a hierarchical splitting of the conforming finite element space of piecewise linear functions that allows to meet the optimality conditions for the related AMLI preconditioner in case of second-order elliptic problems with non-vanishing zero-order term. The finite element method (FEM) then leads to a system of linear equations with a system matrix that is a weighted sum of stiffness and mass matrices. Bounds for the constant \(\gamma \) in the strengthened Cauchy–Bunyakowski–Schwarz inequality are computed for both mass and stiffness matrices in case of a general \(m\) -refinement. Moreover, an additive preconditioner is presented for the pivot blocks that arise in the course of the multilevel block factorization. Its optimality is proven for the case \(m=3\) . Together with the estimates for \(\gamma \) this shows that the construction of a uniformly convergent AMLI method with optimal complexity is possible (for \(m \ge 3\) ). Finally, we discuss the practical application of this preconditioning technique in the context of time-periodic parabolic optimal control problems.     

Spectral Numerical Exterior Calculus Methods for Differential Equations on Radial Manifolds

We analyse a class of nonoverlapping domain decomposition preconditioners for nonsymmetric linear systems arising from discontinuous Galerkin finite element approximations of fully nonlinear Hamilton–Jacobi–Bellman (HJB) partial differential equations. These nonsymmetric linear systems are uniformly bounded and coercive with respect to a related symmetric bilinear form, that is associated to a matrix \(\mathbf {A}\). In this work, we construct a nonoverlapping domain decomposition preconditioner \(\mathbf {P}\), that is based on \(\mathbf {A}\), and we then show that the effectiveness of the preconditioner for solving the nonsymmetric problems can be studied in terms of the condition number \(\kappa (\mathbf {P}^{-1}\mathbf {A})\). In particular, we establish the bound \(\kappa (\mathbf {P}^{-1}\mathbf {A})\lesssim 1+ p^6 H^3 /q^3 h^3\), where H and h are respectively the coarse and fine mesh sizes, and q and p are respectively the coarse and fine mesh polynomial degrees. This represents the first such result for this class of methods that explicitly accounts for the dependence of the condition number on q; our analysis is founded upon an original optimal order approximation result between fine and coarse discontinuous finite element spaces. Numerical experiments demonstrate the sharpness of this bound. Although the preconditioners are not robust with respect to the polynomial degree, our bounds quantify the effect of the coarse and fine space polynomial degrees. Furthermore, we show computationally that these methods are effective in practical applications to nonsymmetric, fully nonlinear HJB equations under h-refinement for moderate polynomial degrees.     

New block triangular preconditioner for linear systems arising from the discretized time-harmonic Maxwell equations

In this paper, based on the preconditioners presented by Rees and Greif [T. Rees, C. Greif, A preconditioner for linear systems arising from interior point optimization methods, SIAM J. Sci. Comput. 29 (2007) 1992-2007], we present a new block triangular preconditioner applied to the problem of solving linear systems arising from finite element discretization of the mixed formulation of the time-harmonic Maxwell equations (k=0) in electromagnetic problems, since linear systems arising from the corresponding equations and methods have the same matrix block structure. Similar to spectral distribution of the preconditioners presented by Rees and Greif, this paper analyzes the corresponding spectral distribution of the new preconditioners considered in this paper. From the views of theories and applications, the presented preconditioners are as efficient as the preconditioners presented by Rees and Greif to apply. Moreover, numerical experiments are also reported to illustrate the efficiency of the presented preconditioners.     

The non-overlapping domain decomposition multiplicative schwarz method

Once a decomposition of the finite element space V into two or more subspaces is given, e.g., via domain decomposition, a specific Multiplicative Schwarz Method (MSM) and Additive Schwarz Method (ASM) is defined. In this paper, we analyse the MSM for the decomposition induced by the approximate discrete harmonic finite element basis which was introduced in a joint paper of the authors with A. Meyer (1990). The main theorem of the present paper states that a special symmetric version of the MSM with approximate orthoprojections is equivalent to some ASM with specially chosen basic transformation and block preconditioners. From this observation we can benefit twice. Indeed, the MSM-DD-preconditioner can be analysed in the MSM framework and implemented as specific ASM-DD-preconditioner in the parallel PCG method studied previously. Emphasis that we look at the ASM and MSM as techniques for defining and analysing parallel DD preconditioners used then in a parallelized version of the PCG-method which is well suited for computations on MIMD computers with local memory and message passing principle.     

Numerical and computational aspects of some block-preconditioners for saddle point systems

Linear systems with two-by-two block matrices are usually preconditioned by block lower- or upper-triangular systems that require an approximation of the related Schur complement. In this work, in the finite element framework, we consider one special such approximation, namely, the element-wise Schur complement. It is sparse and its construction is perfectly parallelizable, making it an appropriate ingredient when building preconditioners for iterative solvers executed on both distributed and shared memory computer architectures. For saddle point matrices with symmetric positive (semi-)definite blocks we show that the Schur complement is spectrally equivalent to the so-constructed approximation and derive spectral equivalence bounds. We also illustrate the quality of the approximation for nonsymmetric problems, where we observe the same good numerical efficiency.Furthermore, we demonstrate the computational and numerical performance of the corresponding preconditioned iterative solution method on a large scale model benchmark problem originating from the elastic glacial isostatic adjustment model discretized using the finite element method.     

Hierarchical matrix preconditioners for the Oseen equations

Hierarchical matrices provide a technique for the data-sparse approximation and matrix arithmetic of large, fully populated matrices. In particular, approximate inverses as well as approximate LU factorizations of finite element stiffness matrices may be computed and stored in nearly optimal complexity. In this paper, we develop efficient -matrix preconditioners for the Oseen equations. In particular, -matrices will provide efficient preconditioners for the auxiliary (scalar) discrete convection–diffusion and pressure Schur complement problems. We will provide various numerical tests comparing the resulting preconditioners with each other. This work was supported in part by the US Department of Energy under Grant No. DE-FG02-04ER25649 and in part by the National Science Foundation under grant No. DMS-0408950.     

Computation of minimum eigenvalue through minimization of rayleigh's quotient for large sparse matrices using vector computer:

Rayleigh Quotient Minimization methods for the calculation of minimal eigenpair achieve great efficiency when ad hoc preconditioners are employed. The performance of different preconditioners on a vector computer is compared and analyzed from a computational standpoint. The numerical experiments are carried out on large symmetric sparse positive definite matrices arising from finite element discretizations of practical problems. Diagonal, polynomial and Kershaw preconditioners are considered for the generalized and the classical eigenvalue problems. Speed-up factors and MFLOPS are calculated. The results show that a good vectorization level of the computational code is achieved. The speed-up factors obtained with the best schemes are generally high and uniformly distributed. Numerical experiments suggest that the highest efficiency for the solution of the classical problem is achieved in vector computers by diagonal preconditioning. This conclusion differs from the results that can be obtained in scalar computers. All the computations were performed on a CRAY X-MP/48 supercomputer.     

Circumventing the ill-conditioning problem with multiquadric radial basis functions: Applications to elliptic partial differential equations

Madych and Nelson [1] proved multiquadric (MQ) mesh-independent radial basis functions (RBFs) enjoy exponential convergence. The primary disadvantage of the MQ scheme is that it is global, hence, the coefficient matrices obtained from this discretization scheme are full. Full matrices tend to become progressively more ill-conditioned as the rank increases.In this paper, we explore several techniques, each of which improves the conditioning of the coefficient matrix and the solution accuracy. The methods that were investigated are

• 1.(1) replacement of global solvers by block partitioning, LU decomposition schemes,
• 2.(2) matrix preconditioners,
• 3.(3) variable MQ shape parameters based upon the local radius of curvature of the function being solved,
• 4.(4) a truncated MQ basis function having a finite, rather than a full band-width,
• 5.(5) multizone methods for large simulation problems, and
• 6.(6) knot adaptivity that minimizes the total number of knots required in a simulation problem.

The hybrid combination of these methods contribute to very accurate solutions.Even though FEM gives rise to sparse coefficient matrices, these matrices in practice can become very ill-conditioned. We recommend using what has been learned from the FEM practitioners and combining their methods with what has been learned in RBF simulations to form a flexible, hybrid approach to solve complex multidimensional problems.     

All-floating coupled data-sparse boundary and interface-concentrated finite element tearing and interconnecting methods

Efficient and robust tearing and interconnecting solvers for large scale systems of coupled boundary and finite element domain decomposition equations are the main topic of this paper. In order to reduce the complexity of the finite element part from ${\mathcal{O}}((H/h)^{d})$ to ${\mathcal{O}}((H/h)^{d-1})$ , we use an interface-concentrated hp finite element approximation. The complexity of the boundary element part is reduced by data-sparse approximations of the boundary element matrices. Finally, we arrive at a parallel solver whose complexity behaves like ${\mathcal{O}}((H/h)^{d-1})$ up to some polylogarithmic factor, where H, h, and d denote the usual scaling parameters of the subdomains, the minimal discretization parameter of the subdomain boundaries, and the spatial dimension, respectively.     

A general class of shift-splitting preconditioners for non-Hermitian saddle point problems with applications to time-harmonic eddy current models

Based on a general splitting of the (1,1) leading block matrix, we first construct a general class of shift-splitting (GCSS) preconditioners for non-Hermitian saddle point problems. Convergence conditions of the corresponding matrix splitting iteration methods and preconditioning properties of the GCSS preconditioned saddle point matrices are analyzed. Then the GCSS preconditioner is specifically applied to the non-Hermitian saddle point problems arising from the finite element discretizations of the hybrid formulations of the time-harmonic eddy current models. With suitable choices of the splittings, the new GCSS preconditioners are easier to implement and have faster convergence rates than the existing shift-splitting preconditioner and its modified variant. Two numerical examples are presented to verify the theoretical results and show effectiveness of the new proposed preconditioners.     

Adaptive AMG with coarsening based on compatible weighted matching

We introduce a new composite adaptive Algebraic Multigrid (composite \(\alpha \) AMG) method to solve systems of linear equations without a-priori knowledge or assumption on characteristics of near-null components of the AMG preconditioned problem referred to as algebraic smoothness. Our version of \(\alpha \) AMG is a composite solver built through a bootstrap strategy aimed to obtain a desired convergence rate. The coarsening process employed to build each new solver component relies on a pairwise aggregation scheme based on weighted matching in a graph, successfully exploited for reordering algorithms in sparse direct methods to enhance diagonal dominance, and compatible relaxation. The proposed compatible matching process replaces the commonly used characterization of strength of connection in both the coarse space selection and in the interpolation scheme. The goal is to design a method leading to scalable AMG for a wide class of problems that go beyond the standard elliptic Partial Differential Equations (PDEs). In the present work, we introduce the method and demonstrate its potential when applied to symmetric positive definite linear systems arising from finite element discretization of highly anisotropic elliptic PDEs on structured and unstructured meshes. We also report on some preliminary tests for 2D and 3D elasticity problems as well as on problems from the University of Florida Sparse Matrix Collection.     

Some preconditioners for elliptic PDE-constrained optimization problems

For the structured systems of linear equations arising from the Galerkin finite element discretizations of elliptic PDE-constrained optimization problems, some preconditioners are proposed to accelerate the convergence rate of Krylov subspace methods such as GMRES for both cases of the Tikhonov parameter $\beta$ not very small (equal or greater than 1e?6) and sufficiently small (less than 1e?6), respectively. We derive the explicit expressions for the eigenvalues and eigenvectors of the corresponding preconditioned matrices. Numerical results show that the corresponding preconditioned GMRES methods perform and match well with the theoretical results.     

Parallel black box \mathcal {H}-LU preconditioning for elliptic boundary value problems

Hierarchical ( $\mathcal {H}$ -) matrices provide a data-sparse way to approximate fully populated matrices. The two basic steps in the construction of an $\mathcal {H}$ -matrix are (a) the hierarchical construction of a matrix block partition, and (b) the blockwise approximation of matrix data by low rank matrices. In the context of finite element discretisations of elliptic boundary value problems, $\mathcal {H}$ -matrices can be used for the construction of preconditioners such as approximate $\mathcal {H}$ -LU factors. In this paper, we develop a new black box approach to construct the necessary partition. This new approach is based on the matrix graph of the sparse stiffness matrix and no longer requires geometric data associated with the indices like the standard clustering algorithms. The black box clustering and a subsequent $\mathcal {H}$ -LU factorisation have been implemented in parallel, and we provide numerical results in which the resulting black box $\mathcal {H}$ -LU factorisation is used as a preconditioner in the iterative solution of the discrete (three-dimensional) convection-diffusion equation.     

Block preconditioners for saddle point systems arising from liquid crystal directors modeling

We analyze two types of block preconditioners for a class of saddle point problems arising from the modeling of liquid crystal directors using finite elements. Spectral properties of the preconditioned matrices are investigated, and numerical experiments are performed to assess the behavior of preconditioned iterations using both exact and inexact versions of the preconditioners.     

On some ways of approximating inverses of banded matrices in connection with deriving preconditioners based on incomplete block factorizations

A unified approach of deriving band approximate inverses of band symmetric positive definite matrices is considered. Such band approximations to the inverses of successive Schur complements are required throughout incomplete block factorizations of block-tridiagonal matrices. Such block-tridiagonal matrices arise, for example, in finite element solution of second order elliptic differential equations. A sharp decay rate estimate for inverses of blocktridiagonal symmetric positive definite matrices is given in addition. Numerical tests on a number of model elliptic boundary value problems are presented comparing thus derived preconditioning matrices.