We introduce block versions of the multielimination incomplete LU (ILUM) factorization preconditioning technique for solving general sparse unstructured linear systems. These preconditioners have a multilevel structure and, for certain types of problems, may exhibit properties that are typically enjoyed by multigrid methods. Several heuristic strategies for forming blocks of independent sets are introduced and their relative merits are discussed. The advantages of block ILUM over point ILUM include increased robustness and efficiency. We compare several versions of the block ILUM, point ILUM, and the dual-threshold-based ILUT preconditioners. In particular, tests with some convection-diffusion problems show that it may be possible to obtain convergence that is nearly independent of the Reynolds number as well as of the grid size.
Recently a great deal of attention has been focused on quantum computation following a sequence of results [Bernstein and Vazirani, in Proc. 25th Annual ACM Symposium Theory Comput., 1993, pp. 11--20, SIAM J. Comput., 26 (1997), pp. 1277--1339], [Simon, in Proc. 35th Annual IEEE Symposium Foundations Comput. Sci., 1994, pp. 116--123, SIAM J. Comput., 26 (1997), pp. 1340--1349], [Shor, in Proc. 35th Annual IEEE Symposium Foundations Comput. Sci., 1994, pp. 124--134] suggesting that quantum computers are more powerful than classical probabilistic computers. Following Shor's result that factoring and the extraction of discrete logarithms are both solvable in quantum polynomial time, it is natural to ask whether all of $\NP$ can be efficiently solved in quantum polynomial time. In this paper, we address this question by proving that relative to an oracle chosen uniformly at random with probability 1 the class $\NP$ cannot be solved on a quantum Turing machine (QTM) in time $o(2^{n/2})$. We also show that relative to a permutation oracle chosen uniformly at random with probability 1 the class $\NP \cap \coNP$ cannot be solved on a QTM in time $o(2^{n/3})$. The former bound is tight since recent work of Grover [in {\it Proc.\ $28$th Annual ACM Symposium Theory Comput.}, 1996] shows how to accept the class $\NP$ relative to any oracle on a quantum computer in time $O(2^{n/2})$.
An algorithm for solving large nonlinear optimization problems with simple bounds is described. It is based on the gradient projection method and uses a limited memory BFGS matrix to approximate the Hessian of the objective function. It is shown how to take advantage of the form of the limited memory approximation to implement the algorithm efficiently. The results of numerical tests on a set of large problems are reported.
Motivated by a recent method of Freund [SIAM J. Sci. Comput., 14 (1993), pp. 470–482], who introduced a quasi-minimal residual (QMR) version of the conjugate gradients squared (CGS) algorithm, a QMR variant of the biconjugate gradient stabilized (Bi-CGSTAB) algorithm of van der Vorst that is called QMRCGSTAB, is proposed for solving nonsymmetric linear systems. The motivation for both QMR variants is to obtain smoother convergence behavior of the underlying method. The authors illustrate this by numerical experiments that also show that for problems on which Bi-CGSTAB performs better than CGS, the same advantage carries over to QMRCGSTAB.
We study the eigenvalue problem for a class $\mathcal{H}$ of band matrices which includes as a proper subclass all band matrices with Toeplitz inverses. Toeplitz matrices of this kind occur, for example, as autocorrelation matrices of purely autoregressive stationary time series. A formula is given for the characteristic polynomial $p_n ( \lambda )$ of an nth order matrix $H_n $ in $\mathcal{H}$, with bandwidth $k + 1\leqq n$, as the ratio of $k \times k$ determinants whose entries are polynomials in the zeros of a certain kth degree polynomial which is independent of n and has one coefficient which depends upon $\lambda $. The formula permits the evaluation of $p_n ( \lambda )$ by means of a computation with complexity independent of n. Also given is a formula for the eigenvectors in terms of these zeros and k coefficients which can be obtained by solving a $k \times k$ homogeneous system.
An implementation is presented of the fast multipole method, which uses approximations based on Poisson's formula. Details for the implementation in both two and three dimensions are given. Also discussed is how the multigrid aspect of the fast multipole method can be exploited to yield efficient programming procedures. The issue of the selection of an appropriate refinement level for the method is addressed. Computational results are given that show the importance of good level selection. An efficient technique that can be used to determine an optimal level to choose for the method is presented.
A group of parallel algorithms, and their implementation for solving a special class of nonlinear equations, are discussed. The type of sparsity occurring in these problems, which arise in VLSI design, structural engineering, and many other areas, is called a block bordered structure. The explicit method and several implicit methods are described, and the new corrected implicit method for solving block bordered nonlinear problems is presented. The relationship between the two types of methods is analyzed, and some computational comparisons are performed. Several variations and globally convergent modifications of the implicit method are also described. Parallel implementations of these algorithms for solving block bordered nonlinear equations are described, and experimental results on the Intel hypercube that show the effectiveness of the parallel implicit algorithms are presented. These experiments include a fairly large circuit simulation that leads to a multilevel block bordered system of nonlinear equations.
The rank revealing QR factorization of a rectangular matrix can sometimes be used as a reliable and efficient computational alternative to the singular value decomposition for problems that involve rank determination. This is illustrated by showing how the rank revealing QR factorization can be used to compute solutions to rank deficient least squares problems, to perform subset selection, to compute matrix approximations of given rank, and to solve total least squares problems.
The paper deals with the numerical solution of time-dependent nonisothermal flow problems, governed by the axisymmetric incompressible Boussinesq equations, on array processors. Using finite difference methods for discretization and a pressure correction method in combination with a successive iteration process for linearization and decoupling of variables, the problem is approximated by a sequence of sparse linear systems. The parallel solution of these systems by preconditioned conjugate gradient (PCG) methods and multigrid (MG) methods is discussed. Numerical experiments on an array processor for a number of test problems, derived from a practical application in the field of crystal growth, are reported. Included are comparisons of the implicit Euler scheme and the Crank–Nicolson scheme for time discretization for different flows, as well as comparisons of PCG methods and MG methods.
Consider the method of independent replications with initial transient deletion for generating confidence intervals for 'steady-state' quantities. To produce intervals with good convergence characteristics, the relative growth rates of the number of replications, the length of each replication, and the deletion period must be controlled. Critical rates for these parameters are determined. The applicability of these results to simultaneously running multiple replications on a highly parallel computer is discussed.
In [SIAM J. Numer. Anal., 21 (1984), pp. 285–299], a method was introduced for solving Poisson’s or the biharmonic equation on an irregular region by making use of an integral equation formulation. Because fast solvers were used to extend the solution to an enclosing rectangle, this method avoided many of the standard problems associated with integral equations. The equations that arose were Fredholm integral equations of the second kind with bounded kernels. In this paper iterative methods are used to solve the dense nonsymmetric linear systems arising from the integral equations. Because the matrices are very well conditioned, conjugate gradient-like methods can be used and will converge very rapidly. The methods are very amenable to vectorization and parallelization, and parallel and vector implementations are described on shared memory multiprocessors. Numerical experiments are described and results presented for a three-dimensional interface problem for the Laplacian on a recording head geometry.
Consideration of an abstract improvement algorithm leads to the following principle, which is similar to that underlying iterative refinement: By making judicious use of relatively few high accuracy computations, high accuracy solutions can be obtained very efficiently by the algorithm. This principle is applied specifically to ${\text{GMRES}}(m)$ here; it can be similarly applied to a number of other “restarted” iterative linear methods as well. Results are given for numerical experiments in solving a discretized linear elliptic boundary value problem and in computing a step of an inexact Newton method using finite differences for a discretized nonlinear elliptic boundary value problem.
Preconditioners based on domain decomposition appear natural for the Krylov solution of implicitly discretized partial differential equations (PDEs) on parallel computers. Two-scale preconditioners (involving a global coarse-grid solve, independent solves over interfaces connecting the coarse-grid points, and independent subdomain solves) have been known since the early 1980s to be “near optimal” in the sense of ensuring a bounded, or at most logarithmically growing, iteration count as the mesh is refined. As a result, the refinement of the mesh can be chosen locally on the basis of truncation error, and the granularity of the domain decomposition can be chosen globally on the basis of parallel computing considerations with only mild effects on the convergence rate of the algorithm. However, overall computational complexity depends not only on the algebraic convergence rate, but also on the operation counts of the components of the preconditioner that must be applied at each iteration. The costs of solving the subdomain systems and the crosspoint system show superlinear growth in their respective (and inversely related) sizes. On the subdomains, the superlinear terms arise from arithmetic only; in the crosspoint system the cost of nonlocal data exchange is also superlinear. These factors make the preconditioner granularity and the choice of its components problem- and machine-dependent compromises. The tradeoffs involved are illustrated through numerical experiments on both shared- and distributed-memory computers for convection-diffusion problems. Because of the development of boundary layers, these problems benefit from local mesh refinement, which is straightforward to accommodate within the domain decomposition framework in a locally uniform sense, but which introduces load balancing as a further consideration in selecting the granularity of the preconditioner. In spite of the tradeoffs, cumulative speedups are obtainable out to at least medium-scale granularity (up to 64 processors in our tests). The largest problems involve $\mathcal{O}(10^5 )$ unknowns partitioned into $\mathcal{O}(10^3 )$ subdomains and converge in $\mathcal{O}(10)$ iterations requiring $\mathcal{O}(1)$ seconds on the Intel iPSC/860.
This paper presents a modification of the block Householder method based on the compact WY representation [R. Schreiber and C. Van Loan, SIAM J. Sci. Statist. Comput., 10 (1989), pp. 52–57]. It is modified in order to introduce more matrix–matrix operations.
Conjugate gradient-type methods for the solution of large sparse linear systems Ax = b with complex symmetric coefficient matrices A = A(T) are considered. Such linear systems arise in important applications, such as the numerical solution of the complex Helmholtz equation. Furthermore, most complex non-Hermitian linear systems which occur in practice are actually complex symmetric. Conjugate gradient-type iterations which are based on a variant of the nonsymmetric Lanczos algorithm for complex symmetric matrices are investigated. In particular, a new approach with iterates defined by a quasi-minimal residual property is proposed. The resulting algorithm presents several advantages over the standard biconjugate gradient method. Some remarks are also included on the obvious approach to general complex linear systems by solving equivalent real linear systems for the real and imaginary parts of x. Finally, numerical experiments for linear systems arising from the complex Helmholtz equation are reported.
Many elliptic partial differential equations can be solved numerically with near optimal efficiency through the uses of adaptive refinement and multigrid solution techniques. This paper presents a more unified approach to the combined process of adaptive refinement and multigrid solution which can be used with high order finite elements. Refinement is achieved by the bisection of pairs of triangles, corresponding to the addition of one or more basis functions to the approximation space. An approximation of the resulting change in the solution is used as an error indicator. The multigrid iteration uses red-black Gauss-Seidel relaxation with local black relaxations. The grid transfers use the change between the nodal and hierarchical bases. This multigrid iteration requires only O(N) operations, even for highly nonuniform grids, and is defined for any finite element space. The full multigrid method is an optimal blending of the processes of adaptive refinement and multigrid iteration. To minimize the number of operations required, the duration of the refinement phase is based on increasing the dimension of the approximation space by the largest possible factor, given the error reduction of the multigrid iteration. The algorithm (i) uses only O(N) operations, (ii) solves the discrete system to the accuracy of the discretization error, and (iii) achieves optimal convergence of the discretization error in the presence of singularities. Numerical experiments confirm this for linear, quadratic, and cubic elements.
The truncated singular value decomposition (SVD) method is useful for solving the standard-form regularization problem: $\min ||{\bf x}||_2 $ subject to $\min ||A{\bf x} - {\bf b}||_2 $. This paper presents a modification of the truncated SVD method, which solves the more general problem: $\min ||L{\bf x}||_2 $ subject to $\min ||A{\bf x} - {\bf b}||_2 $, where L is a general matrix with full row rank. The extra work, associated with the introduction of the matrix L, is dominated by a QR-factorization of a matrix with dimensions smaller than those of L. In order to determine the optimal solution, it is often necessary to compute a sequence of regularized solutions, and it is shown how this can be accomplished with little extra computational effort. Finally, the new method is illustrated with an example from helioseismology.
A domain decomposition algorithm based on a hybrid variational principle is developed for the parallel finite element solution of selfadjoint elliptic partial differential equations. The spatial domain is partitioned into a set of totally disconnected subdomains, each assigned to an individual processor. Lagrange multipliers are introduced to enforce compatibility at the interface points. Within each subdomain, the singularity due to the disconnection is resolved in a two-step procedure. First, the null space component of each local operator is eliminated from the local problem. Next, its contribution to the local solution is related to the Lagrange multipliers through an orthogonality condition. Finally, a conjugate projected gradient algorithm is developed for the solution of the coupled system of local null space components and Lagrange multipliers. When implemented on local memory multiprocessors, the proposed hybrid method requires fewer interprocessor communications than conventional Schur methods. It is also suitable for parallel/vector computers with shared memory. Moreover, unlike parallel direct solvers, it exhibits a degree of parallelism that is not limited by the bandwidth of the finite element system of equations. In this paper, it is applied to the solution of large-scale structural and solid mechanics problems.
