Preconditioned conjugate gradient methods for large-scale fluid flow applications

The simulation of large-scale fluid flow applications often requires the efficient solution of extremely large nonsymmetric linear and nonlinear sparse systems of equations arising from the discretization of systems of partial differential equations. While preconditioned conjugate gradient methods work well for symmetric, positive-definite matrices, other methods are necessary to treat large, nonsymmetric matrices. The applications may also involve highly localized phenomena which can be addressed via local and adaptive grid refinement techniques. These local refinement methods usually cause non-standard grid connections which destroy the bandedness of the matrices and the associated ease of solution and vectorization of the algorithms. The use of preconditioned conjugate gradient or conjugate-gradient-like iterative methods in large-scale reservoir simulation applications is briefly surveyed. Then, some block preconditioning methods for adaptive grid refinement via domain decomposition techniques are presented and compared. These techniques are being used efficiently in existing large-scale simulation codes.     

Numerical integration using sparse grids

We present new and review existing algorithms for the numerical integration of multivariate functions defined over d-dimensional cubes using several variants of the sparse grid method first introduced by Smolyak [49]. In this approach, multivariate quadrature formulas are constructed using combinations of tensor products of suitable one-dimensional formulas. The computing cost is almost independent of the dimension of the problem if the function under consideration has bounded mixed derivatives. We suggest the usage of extended Gauss (Patterson) quadrature formulas as the one‐dimensional basis of the construction and show their superiority in comparison to previously used sparse grid approaches based on the trapezoidal, Clenshaw–Curtis and Gauss rules in several numerical experiments and applications. For the computation of path integrals further improvements can be obtained by combining generalized Smolyak quadrature with the Brownian bridge construction. This revised version was published online in August 2006 with corrections to the Cover Date.     

Four algorithms for the the efficient computation of truncated pivoted QR approximations to a sparse matrix

Summary. In this paper we propose four algorithms to compute truncated pivoted QR approximations to a sparse matrix. Three are based on the Gram–Schmidt algorithm and the other on Householder triangularization. All four algorithms leave the original matrix unchanged, and the only additional storage requirements are arrays to contain the factorization itself. Thus, the algorithms are particularly suited to determining low-rank approximations to a sparse matrix. Received February 23, 1998 / Revised version received April 16, 1998     

A bidiagonalization algorithm for solving large and sparse ill-posed systems of linear equations

An iterative method based on Lanczos bidiagonalization is developed for computing regularized solutions of large and sparse linear systems, which arise from discretizations of ill-posed problems in partial differential or integral equations. Determination of the regularization parameter and termination criteria are discussed. Comments are given on the computational implementation of the algorithm.Dedicated to Peter Naur on the occasion of his 60th birthday     

A sparse counterpart of Reichel and Gragg’s package QRUP

We describe how to maintain the triangular factor of a sparse QR factorization when columns are added and deleted and Q cannot be stored for sparsity reasons. The updating procedures could be thought of as a sparse counterpart of Reichel and Gragg’s package QRUP. They allow us to solve a sequence of sparse linear least squares subproblems in which each matrix B_k is an independent subset of the columns of a fixed matrix A, and B_k+1 is obtained by adding or deleting one column. Like Coleman and Hulbert [T. Coleman, L. Hulbert, A direct active set algorithm for large sparse quadratic programs with simple bounds, Math. Program. 45 (1989) 373-406], we adapt the sparse direct methodology of Björck and Oreborn of the late 1980s, but without forming A^TA, which may be only positive semidefinite. Our Matlab 5 implementation works with a suitable row and column numbering within a static triangular sparsity pattern that is computed in advance by symbolic factorization of A^TA and preserved with placeholders.     

Solving sparse linear least-squares problems on some supercomputers by using large dense blocks

Efficient subroutines for dense matrix computations have recently been developed and are available on many high-speed computers. On some computers the speed of many dense matrix operations is near to the peak-performance. For sparse matrices storage and operations can be saved by operating only and storing only nonzero elements. However, the price is a great degradation of the speed of computations on supercomputers (due to the use of indirect addresses, to the need to insert new nonzeros in the sparse storage scheme, to the lack of data locality, etc.). On many high-speed computers a dense matrix technique is preferable to sparse matrix technique when the matrices are not large, because the high computational speed compensates fully the disadvantages of using more arithmetic operations and more storage. For very large matrices the computations must be organized as a sequence of tasks in each of which a dense block is treated. The blocks must be large enough to achieve a high computational speed, but not too large, because this will lead to a large increase in both the computing time and the storage. A special “locally optimized reordering algorithm” (LORA) is described, which reorders the matrix so that dense blocks can be constructed and treated with some standard software, say LAPACK or NAG. These ideas are implemented for linear least-squares problems. The rectangular matrices (that appear in such problems) are decomposed by an orthogonal method. Results obtained on a CRAY C92A computer demonstrate the efficiency of using large dense blocks.     

Fouriertransformation auf dünnen Gittern mit hierarchischen Basen

Zusammenfassung Es wird ein auf der schnellen Fouriertransformation beruhender Algorithmus zur trigonometrischen Interpolation von multivariaten Funktionen auf einem dünnen Gitter beschrieben. Bei nur geringfügig schlechterer Approximationsgenauigkeit für Funktionen aus Korobovräumen hat dieser Algorithmus eine viel kleinere Komplexität als die Algorithmen, die auf gewöhnlichen Gittern arbeiten. Die Transformation ist auf einfache Weise umkehrbar.

---

Fouriertransform on sparse grids with hierarchical bases

Summary An algorithm for the trigonometric interpolation of functions ofn variables on a sparse grid is described. This discrete Fourier transform based on FFT has a greatly reduced complexity in comparison to the fourier transform on a regular grid whereas the approximation quality is only slightly reduced if the function belongs to a Korobov space. The transformation is easily invertible.

     

Incomplete Gröbner basis as a preconditioner for polynomial systems

Precondition plays a critical role in the numerical methods for large and sparse linear systems. It is also true for nonlinear algebraic systems. In this paper incomplete Gröbner basis (IGB) is proposed as a preconditioner of homotopy methods for polynomial systems of equations, which transforms a deficient system into a system with the same finite solutions, but smaller degree. The reduced system can thus be solved faster. Numerical results show the efficiency of the preconditioner.     

Estimating the largest singular values of large sparse matrices via modified moments

We describe a procedure for determining a few of the largest singular values of a large sparse matrix. The method by Golub and Kent which uses the method of modified moments for estimating the eigenvalues of operators used in iterative methods for the solution of linear systems of equations is appropriately modified in order to generate a sequence of bidiagonal matrices whose singular values approximate those of the original sparse matrix. A simple Lanczos recursion is proposed for determining the corresponding left and right singular vectors. The potential asynchronous computation of the bidiagonal matrices using modified moments with the iterations of an adapted Chebyshev semi-iterative (CSI) method is an attractive feature for parallel computers. Comparisons in efficiency and accuracy with an appropriate Lanczos algorithm (with selective re-orthogonalization) are presented on large sparse (rectangular) matrices arising from applications such as information retrieval and seismic reflection tomography. This procedure is essentially motivated by the theory of moments and Gauss quadrature.This author's work was supported by the National Science Foundation under grants NSF CCR-8717492 and CCR-910000N (NCSA), the U.S. Department of Energy under grant DOE DE-FG02-85ER25001, and the Air Force Office of Scientific Research under grant AFOSR-90-0044 while at the University of Illinois at Urbana-Champaign Center for Supercomputing Research and Development.This author's work was supported by the U.S. Army Research Office under grant DAAL03-90-G-0105, and the National Science Foundation under grant NSF DCR-8412314.     

Sparse grid spaces for the numerical solution of the electronic Schrödinger equation

This article complements the author’s recent work [Numer. Math. 98, 731–759 (2004)] on the regularity of the electronic Schrödinger equation in Hilbert spaces of mixed derivatives. It has been shown there that the solutions of this equation are surprisingly smooth and possess square integrable mixed weak derivatives of order up to N+1 with N the number of electrons across the singularities of the interaction potentials, and it has been claimed that this result can help to break the complexity barriers in computational quantum mechanics using correspondingly antisymmetrized sparse grid trial functions. A construction of this kind that can be interpreted as a sparse grid sampling theorem is sketched here.     

A cascadic multigrid algorithm for semilinear elliptic problems

Summary. We propose a cascadic multigrid algorithm for a semilinear elliptic problem. The nonlinear equations arising from linear finite element discretizations are solved by Newton's method. Given an approximate solution on the coarsest grid on each finer grid we perform exactly one Newton step taking the approximate solution from the previous grid as initial guess. The Newton systems are solved iteratively by an appropriate smoothing method. We prove that the algorithm yields an approximate solution within the discretization error on the finest grid provided that the start approximation is sufficiently accurate and that the initial grid size is sufficiently small. Moreover, we show that the method has multigrid complexity. Received February 12, 1998 / Revised version received July 22, 1999 / Published online June 8, 2000     

Stability and accuracy of adapted finite element methods for singularly perturbed problems

The stability and accuracy of a standard finite element method (FEM) and a new streamline diffusion finite element method (SDFEM) are studied in this paper for a one dimensional singularly perturbed connvection-diffusion problem discretized on arbitrary grids. Both schemes are proven to produce stable and accurate approximations provided that the underlying grid is properly adapted to capture the singularity (often in the form of boundary layers) of the solution. Surprisingly the accuracy of the standard FEM is shown to depend crucially on the uniformity of the grid away from the singularity. In other words, the accuracy of the adapted approximation is very sensitive to the perturbation of grid points in the region where the solution is smooth but, in contrast, it is robust with respect to perturbation of properly adapted grid inside the boundary layer. Motivated by this discovery, a new SDFEM is developed based on a special choice of the stabilization bubble function. The new method is shown to have an optimal maximum norm stability and approximation property in the sense that where u _N is the SDFEM approximation in linear finite element space V ^N of the exact solution u. Finally several optimal convergence results for the standard FEM and the new SDFEM are obtained and an open question about the optimal choice of the monitor function for the moving grid method is answered. This work was supported in part by NSF DMS-0209497 and NSF DMS-0215392 and the Changjiang Professorship through Peking University.     

Patch-adaptive multilevel iteration

Themultilevel adaptive iteration is an attempt to improve both the robustness and efficiency of iterative sparse system solvers. Unlike in most other iterative methods, the order of processing and sequence of operations is not determined a priori. The method consists of a relaxation scheme with an active set strategy and can be viewed as an efficient implementation of the Gauß-Southwell relaxation. With this strategy, computational work is focused on where it can efficiently improve the solution quality. To obtain full efficiency, the algorithm must be used on a multilevel structure. This algorithm is then closely related to multigrid or multilevel preconditioning algorithms, and can be shown to have asymptotically optimal convergence. In this paper the focus is on a variant that uses data structures with a locally uniform grid refinement. The resulting grid system consists of a collection of patches where each patch is a uniform rectangular grid and where adaptive refinement is accomplished by arranging the patches flexibly in space. This construction permits improved implementations that better exploit high performance computer designs. This will be demonstrated by numerical examples.     

A generalized successive overrelaxation method for least squares problems

In this paper a new iterative method is given for solving large sparse least squares problems and computing the minimum norm solution to underdetermined consistent linear systems. The new scheme is called the generalized successive overrelaxation (GSOR) method and is shown to be convergent ifA is full column rank. The GSOR method involves a parameter ρ and an auxiliary matrixP. One can choose matrix P so that the GSOR method only involves matrix and vector operations; therefore the GSOR method is suitable for parallel computations. Besides, the GSOR method can be combined with preconditioning techniques, and therefore can be expected to be more effective. This author's work was supported by Natural Science Foundation of Liaoning Province, China.     

A survey of multilevel preconditioned iterative methods

We survey multilevel iterative methods applied for solving large sparse systems with matrices, which depend on a level parameter, such as arise by the discretization of boundary value problems for partial differential equations when successive refinements of an initial discretization mesh is used to construct a sequence of nested difference or finite element meshes.We discuss various two-level (two-grid) preconditioning techniques, including some for nonsymmetric problems. The generalization of these techniques to the multilevel case is a nontrivial task. We emphasize several ways this can be done including classical multigrid methods and a recently proposed algebraic multilevel preconditioning method. Conditions for which the methods have an optimal order of computational complexity are presented.On leave from the Institute of Mathematics, and Center for Informatics and Computer Technology, Bulgarian Academy of Sciences, Sofia, Bulgaria. The research of the second author reported here was partly supported by the Stichting Mathematisch Centrum, Amsterdam.     

On additive Schwarz preconditioners for sparse grid discretizations

Summary Based on the framework of subspace splitting and the additive Schwarz scheme, we give bounds for the condition number of multilevel preconditioners for sparse grid discretizations of elliptic model problems. For a BXP-like preconditioner we derive an estimate of the optimal orderO(1) and for a HB-like variant we obtain an estimate of the orderO(k ² ·2 ^k/2 ), wherek denotes the number of levels employed. Furthermore, we confirm these results by numerically computed condition numbers.     

Stable three-point wavelet bases on general meshes

Summary. Recently, we introduced a wavelet basis on general, possibly locally refined linear finite element spaces. Each wavelet is a linear combination of three nodal basis functions, independently of the number of space dimensions. In the present paper, we show -stability of this basis for a range of , that in any case includes , which means that the corresponding additive Schwarz preconditioner is optimal for second order problems. Furthermore, we generalize the construction of the wavelet basis to manifolds. We show that the wavelets have at least one-, and in areas where the manifold is smooth and the mesh is uniform even two vanishing moments. Because of these vanishing moments, apart from preconditioning, the basis can be used for compression purposes: For a class of integral equation problems, the stiffness matrix with respect to the wavelet basis will be close to a sparse one, in the sense that, a priori, it can be compressed to a sparse matrix without the order of convergence being reduced. Received November 6, 1996 / Revised version received June 30, 1997     

Robust iterative methods for elliptic problems with highly varying coefficients in thin substructures

Summary. In this paper we introduce a class of robust multilevel interface solvers for two-dimensional finite element discrete elliptic problems with highly varying coefficients corresponding to geometric decompositions by a tensor product of strongly non-uniform meshes. The global iterations convergence rate is shown to be of the order with respect to the number of degrees of freedom on the single subdomain boundaries, uniformly upon the coarse and fine mesh sizes, jumps in the coefficients and aspect ratios of substructures. As the first approach, we adapt the frequency filtering techniques [28] to construct robust smoothers on the highly non-uniform coarse grid. As an alternative, a multilevel averaging procedure for successive coarse grid correction is proposed and analyzed. The resultant multilevel coarse grid preconditioner is shown to have (in a two level case) the condition number independent of the coarse mesh grading and jumps in the coefficients related to the coarsest refinement level. The proposed technique exhibited high serial and parallel performance in the skin diffusion processes modelling [20] where the high dimensional coarse mesh problem inherits a strong geometrical and coefficients anisotropy. The approach may be also applied to magnetostatics problems as well as in some composite materials simulation. Received December 27, 1994     

Fourier analysis of a robust multigrid method for convection-diffusion equations

Summary. We consider a two-grid method for solving 2D convection-diffusion problems. The coarse grid correction is based on approximation of the Schur complement. As a preconditioner of the Schur complement we use the exact Schur complement of modified fine grid equations. We assume constant coefficients and periodic boundary conditions and apply Fourier analysis. We prove an upper bound for the spectral radius of the two-grid iteration matrix that is smaller than one and independent of the mesh size, the convection/diffusion ratio and the flow direction; i.e. we have a (strong) robustness result. Numerical results illustrating the robustness of the corresponding multigrid -cycle are given. Received October 14, 1994     

On the numerical solution of diffusion–reaction equations with singular source terms

A numerical study is presented of reaction–diffusion problems having singular reaction source terms, singular in the sense that within the spatial domain the source is defined by a Dirac delta function expression on a lower dimensional surface. A consequence is that solutions will be continuous, but not continuously differentiable. This lack of smoothness and the lower dimensional surface form an obstacle for numerical discretization, including amongst others order reduction. In this paper the standard finite volume approach is studied for which reduction from order two to order one occurs. A local grid refinement technique is discussed which overcomes the reduction.