Large Systems Of Linear Equations Research Articles

Wavelet-In-Time Multigrid-In-Space Preconditioning of Parabolic Evolution Equations

Two space-time variational formulations of linear parabolic evolution equations are considered: one is symmetric and elliptic on the trial space, while the other is not. In each case, a space-time Petrov--Galerkin discretization using suitable tensor product trial and test functions leads to a large linear system of equations. The well-posedness of this system with respect to parabolic norms induces a canonical preconditioner for the algebraic equations that arise after a choice of basis. For the iterative resolution of this algebraic system with parallelization in the temporal direction, we propose a sparse algebraic wavelet-in-time transformation on possibly nonuniform temporal meshes. This transformation approximately block-diagonalizes the preconditioner, and the individual spatial blocks can then be inverted, for instance, by standard spatial multigrid methods in parallel. The performance of the preconditioner is documented in a series of numerical experiments.

Restarted Pulay mixing for efficient and robust acceleration of fixed-point iterations

We present a variant of the restarted Pulay's Direct Inversion in the Iterative Subspace (DIIS) method for efficiently and robustly accelerating the convergence of fixed-point iterations. Specifically, we propose a simple modification of DIIS without any additional parameters, which we refer to as the r-Pulay method. We demonstrate the efficacy of r-Pulay in the context of the Jacobi iteration for solving large linear systems of equations, as well as in the Self Consistent Field (SCF) approach for Density Functional Theory (DFT) calculations. Overall, we find r-Pulay to be an attractive version of the restarted DIIS method.

Some properties of range restricted GMRES methods

The GMRES method is one of the most popular iterative schemes for the solution of large linear systems of equations with a square nonsingular matrix. GMRES-type methods also have been applied to the solution of linear discrete ill-posed problems. Computational experience indicates that for the latter problems variants of the standard GMRES method, that require the solution to live in the range of a positive power of the matrix of the linear system of equations to be solved, generally yield more accurate approximations of the desired solution than standard GMRES. This paper investigates properties of these variants of GMRES.

An Efficient Parallel Gauss-Seidel Algorithm on a 3D Torus Network-on-Chip

Network-on-chip (NoC) multi-core architectures with a large number of processing elements are becoming a reality with the recent developments in technology. In these modern systems the processing elements are interconnected with regular NoC topologies such as meshes and tori. In this paper we propose a parallel Gauss-Seidel (GS) iterative algorithm for solving large systems of linear equations on a 3-dimensional torus NoC architecture. The proposed parallel algorithm is O(Nn2/k3) time complexity for solving a system with a matrix of order n on a k×k×k 3D torus NoC architecture with N iterations assuming n and N are large compared to k. We show that under these conditions the proposed parallel GS algorithm has near optimal speedup.

An effective method to suppress numerical dispersion in 2D acoustic and elastic modelling using a high-order Padé approximation

In this paper, we propose a method that can effectively reduce the numerical dispersion for solving the acoustic and elastic-wave equations. The method is a fourth-order Pade approximation scheme, in which the time difference operator is a rational function and a block tridiagonal system needs to be solved at each step. On the one hand, to efficiently solve a large linear system of equations we propose an explicit method for this implicit algrithom. On the other hand, to approximate the high-order spatial derivatives we use an eighth-order stereo-modelling method using wavefield displacements and their gradients simultaneously. For this new method, we investigate some mathematical properties including the stability, errors and the numerical dispersion relationship for 1D and 2D cases. We also present numerical results computed by the Pade approximation and compare them with the eighth-order Lax–Wendroff correction method and the eighth-order staggered-grid method. Numerical results show that the high-order Pade approximation scheme can effectively suppress the numerical dispersion caused by discretizing the wave equations when coarse spatial grids are used or models have strong velocity contrasts between adjacent grids. In contrast to other high-order finite-difference methods, the new method takes substantially less computational time and requires less memory because large spatial and time increments can be used. Thus the high-order Pade approximation method can potentially be used to solve large-scale wave propagation problems and seismic tomography based on the wave equations.

Tikhonov regularization via flexible Arnoldi reduction

Flexible GMRES, introduced by Saad, is a generalization of the standard GMRES method for the solution of large linear systems of equations. It is based on the flexible Arnoldi process for reducing a large square matrix to a small matrix. We describe how the flexible Arnoldi process can be applied to implement one-parameter and multi-parameter Tikhonov regularization of linear discrete ill-posed problems. The method proposed is well suited for large-scale problems. Moreover, computed examples show that our method can give approximate solutions of higher accuracy than available direct methods for small-scale problems.

Space-Time Regularization for Video Decompression

We consider the problem of reconstructing frames from a video which has been compressed using the video compressive sensing (VCS) method. In VCS data, each frame comes from first subsampling the original video data in space and then averaging the subsampled sequence in time. This results in a large linear system of equations whose inversion is ill-posed. We introduce a convex regularizer to invert the system, where the spatial component is regularized by the total variation seminorm, and the temporal component is regularized by enforcing sparsity on the difference between the spatial gradients of each frame. Since the regularizers are $L^1$-like norms, the model can be written in the form of an easy-to-solve saddle point problem. The saddle point problem is solved by the primal-dual algorithm, whose implementation calls for nearly pointwise operations (i.e., no direct linear inversion) and has a simple parallel version. Results show that our model decompresses videos more accurately than other popular mod...

Accurate mirror-based camera pose estimation with explicit geometric meanings

An accurate closed form solution is proposed to estimate camera pose by several mirrored reference object images acquired via a planar mirror under different unknown poses. Compared with state-of-the-art methods, our method is more accurate when there are more than three images and has explicit geometric meanings. This method also properly handles cases in which some of the mirror poses are parallel. The central idea is to minimize an error metric based on all reflections of rotation, which enables the camera rotation to be estimated directly by SVD of sum of mirrored camera rotations. After that, the camera translation is computed by solving a large system of linear equations to minimize object space collinearity error. Both synthesized data and real data experiments show the advantages of our approach.

On facilitating the use of HARDI in population studies by creating rotation-invariant markers.

We design and evaluate a novel method to compute rotationally invariant features using High Angular Resolution Diffusion Imaging (HARDI) data. These measures quantify the complexity of the angular diffusion profile modeled using a higher order model, thereby giving more information than classical diffusion tensor-derived parameters. The method is based on the spherical harmonic (SH) representation of the angular diffusion information, and is generalizable to a range of HARDI reconstruction models. These scalars are obtained as homogeneous polynomials of the SH representation of a HARDI reconstruction model. We show that finding such polynomials is equivalent to solving a large linear system of equations, and present a numerical method based on sparse matrices to efficiently solve this system. Among the solutions, we only keep a subset of algebraically independent polynomials, using an algorithm based on a numerical implementation of the Jacobian criterion. We compute a set of 12 or 25 rotationally invariant measures representative of the underlying white matter for the rank-4 or rank-6 spherical harmonics (SHs) representation of the apparent diffusion coefficient (ADC) profile, respectively. Synthetic data was used to investigate and quantify the difference in contrast. Real data acquired with multiple repetitions showed that within subject variation in the invariants was less than the difference across subjects - facilitating their use to study population differences. These results demonstrate that our measures are able to characterize white matter, especially complex white matter found in regions of fiber crossings and hence can be used to derive new biomarkers for HARDI and can be used for HARDI-based population analysis.

Large Scale Finite Element Analysis Via Assembly-Free Deflated Conjugate Gradient

Large-scale finite element analysis (FEA) with millions of degrees of freedom (DOF) is becoming commonplace in solid mechanics. The primary computational bottleneck in such problems is the solution of large linear systems of equations. In this paper, we propose an assembly-free version of the deflated conjugate gradient (DCG) for solving such equations, where neither the stiffness matrix nor the deflation matrix is assembled. While assembly-free FEA is a well-known concept, the novelty pursued in this paper is the use of assembly-free deflation. The resulting implementation is particularly well suited for large-scale problems and can be easily ported to multicore central processing unit (CPU) and graphics-programmable unit (GPU) architectures. For demonstration, we show that one can solve a 50 × 106 degree of freedom system on a single GPU card, equipped with 3 GB of memory. The second contribution is an extension of the “rigid-body agglomeration” concept used in DCG to a “curvature-sensitive agglomeration.” The latter exploits classic plate and beam theories for efficient deflation of highly ill-conditioned problems arising from thin structures.

Polynomial-based mean weighted residuals methods for elliptic problems with nonlocal boundary conditions in the rectangle

In this paper, polynomial-based mean weighted residuals methods for the solution of elliptic problems with nonlocal boundary conditions, in rectangular domains, are numerically studied. When these methods are employed, the nonclassical boundary conditions involve the solution of large systems of linear equations or least squares problems, hence some numerical techniques for these solvers are compared to show the importance of using efficient algorithms for this purpose. Different kinds of nodes are used to demonstrate how they can be employed to solve different numerical problems when large derivatives of the solution appear. We will also study how using extra precision and/or oversampling can often reduce the computational effort.These methods can also be combined with others (as, for example, finite difference or spline methods) to solve linear and nonlinear parabolic or hyperbolic partial differential equations. The numerical study of some techniques as those explained above can help to obtain significantly better numerical approximations with a smaller computational effort.

Hydraulic Analysis of Water Supply Networks Using a Modified Hardy Cross Method

There are different methods for the hydraulic analysis of water supply networks. In the solution process of most of these methods, a large system of linear equations is solved in each iteration. This usually requires a high computational effort. Hardy Cross method is one of the approaches that do not need such a process and may converge to the solution through scalar divisions. However, this method has two shortcomings: first, initial discharges should satisfy continuity equation at each node; second a large number of iterations are required to converge to solution. In this article an algorithm is suggested for initial discharges that are close to the final results while the continuity equations are automatically established. This algorithm may be directly implemented in the Hardy Cross method. To reduce the number of iterations the Hardy Cross method is combined with third-order and sixteenth-order methods. The results of some numerical examples demonstrate that the use of the combined approach with the suggested initial guess reduces the number of iterations and hydraulic analysis time and the solutions converge with a high accuracy.

A coarse space for heterogeneous Helmholtz problems based on the Dirichlet-to-Neumann operator

The Helmholtz equation governing wave propagation and scattering phenomena is difficult to solve numerically. Its discretization with piecewise linear finite elements results in typically large linear systems of equations. The inherently parallel domain decomposition methods constitute hence a promising class of preconditioners. An essential element of these methods is a good coarse space. Here, the Helmholtz equation presents a particular challenge, as even slight deviations from the optimal choice can be devastating.In this paper, we present a coarse space that is based on local eigenproblems involving the Dirichlet-to-Neumann operator. Our construction is completely automatic, ensuring good convergence rates without the need for parameter tuning. Moreover, it naturally respects local variations in the wave number and is hence suited also for heterogeneous Helmholtz problems. The resulting method is parallel by design and its efficiency is demonstrated on 2D homogeneous and heterogeneous numerical examples.

A Multiple-FPGA parallel computing architecture for real-time simulation of soft-object deformation

Hardware-based parallel computing is proposed for acceleration of finite-element (FE) analysis of linear elastic deformation models. An implementation of the Preconditioned Conjugate Gradient algorithm on N Field Programmable Gate Array (FPGA) devices solves the large linear system of equations arising from the FE discretization. The system employs a large number of customized fixed-point computing units with a high-throughput memory architecture. An implementation of this scalable architecture on four Altera EP3SE110 FPGA devices yields a peak performance of 604 Giga Operations per second. This enables haptic simulation of a 3-dimensional deformable object of 21000 elements at an update rate of 400Hz.

Vertex-Ball Spring Smoothing: An efficient method for unstructured dynamic hybrid meshes

Spring analogy approach is one of the most popular dynamic mesh deformation methods. In the ball-vertex method, perpendicular linear springs are introduced to deal with the element collapse problem. However, it is not very efficient as a large system of linear equations has to be resolved. In order to overcome this difficulty, the Vertex-Ball Spring Smoothing algorithm (VerBSS) is proposed in this paper. Following the mesh smoothing concept, a sub-spring system derived from the “ball-vertex” model is built and solved on a node by node basis using an LDLT solver. Interior nodes are smoothed layer by layer in an iterative manner to achieve the best result. Parallel scheme is also introduced in the smoothing process for further improvement of the efficiency. Numerical examples in two and three dimensions show that VerBSS is much more efficient than the ball-vertex method, and is capable of dealing with practical engineering objects with complex geometries subject to large deformations. VerBSS can be applied to complicated mesh topologies as well, not only to 2D/3D dynamic mesh, but also to the hybrid dynamic mesh.

Performance evaluation of kernel fusion BLAS routines on the GPU: iterative solvers as case study

Programmers usually implement iterative methods that solve partial differential equations by expressing them using a sequence of basic kernels from libraries optimized for the graphics processing unit (GPU). The global runtime of the resulting combination is often penalized by the smallest and most inefficient vector operations. To improve the GPU exploitation, we identify and analyze the potential kernels to be fused according to the data dependence, data type and size, and GPU resources. This paper provides an extensive analysis of the impact of fusing vector operations [level 1 of Basic Linear Algebra Subprograms (BLAS)] on the performance of the GPU. The experimental evaluation shows that this optimization provides noticeable improvement especially for kernels with lower memory requirements and on more modern GPUs. It is worth noting that the fused BLAS operations can be very useful to help programmers efficiently code iterative methods to solve large linear systems of equations for the GPU. Iterative methods such as biconjugate gradient method (BCG) are one of the examples that can benefit from this optimization strategy. Indeed, kernel fusion of vector routines makes the most efficient GPU implementation of BCG run between $$1.09\times $$ 1.09 × and $$1.27\times $$ 1.27 × faster on three GPUs of different characteristics.

Comment on “Plasmonic nanoparticle monomers and dimers: from nanoantennas to chiral metamaterials” by Chigrin et al. in Appl Phys B (2011) 105:81–97

Chigrin et al. (Appl Phys B 105:81, 2011) have recently presented a general direct solution of the integral equation for the electric field inside the scatterer, which does not require inversion of integral operator or, equivalently, solution of large system of linear equations. This comment points out that the rigorous applicability of the considered derivation is limited to ellipsoids much smaller than the wavelength.

FGMRES for linear discrete ill-posed problems

GMRES is one of the most popular iterative methods for the solution of large linear systems of equations. However, GMRES does not always perform well when applied to the solution of linear systems of equations that arise from the discretization of linear ill-posed problems with error-contaminated data represented by the right-hand side. Such linear systems are commonly referred to as linear discrete ill-posed problems. The FGMRES method, proposed by Saad, is a generalization of GMRES that allows larger flexibility in the choice of solution subspace than GMRES. This paper explores application of FGMRES to the solution of linear discrete ill-posed problems. Numerical examples illustrate that FGMRES with a suitably chosen solution subspace may determine approximate solutions of higher quality than commonly applied iterative methods.

Asynchronous collision integrators: Explicit treatment of unilateral contact with friction and nodal restraints

This article presents asynchronous collision integrators and a simple asynchronous method treating nodal restraints. Asynchronous discretizations allow individual time step sizes for each spatial region, improving the efficiency of explicit time stepping for finite element meshes with heterogeneous element sizes. The article first introduces asynchronous variational integration being expressed by drift and kick operators. Linear nodal restraint conditions are solved by a simple projection of the forces that is shown to be equivalent to RATTLE. Unilateral contact is solved by an asynchronous variant of decomposition contact response. Therein, velocities are modified avoiding penetrations. Although decomposition contact response is solving a large system of linear equations (being critical for the numerical efficiency of explicit time stepping schemes) and is needing special treatment regarding overconstraint and linear dependency of the contact constraints (for example from double-sided node-to-surface contact or self-contact), the asynchronous strategy handles these situations efficiently and robust. Only a single constraint involving a very small number of degrees of freedom is considered at once leading to a very efficient solution. The treatment of friction is exemplified for the Coulomb model. Special care needs the contact of nodes that are subject to restraints. Together with the aforementioned projection for restraints, a novel efficient solution scheme can be presented. The collision integrator does not influence the critical time step. Hence, the time step can be chosen independently from the underlying time-stepping scheme. The time step may be fixed or time-adaptive. New demands on global collision detection are discussed exemplified by position codes and node-to-segment integration. Numerical examples illustrate convergence and efficiency of the new contact algorithm. Copyright © 2013 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons, Ltd.

Communication Avoiding ILU0 Preconditioner

In this paper we present a communication avoiding ILU0 preconditioner for solving large linear systems of equations by using iterative Krylov subspace methods. Recent research has focused on communication avoiding Krylov subspace methods based on so-called s-step methods. However, there are not many communication avoiding preconditioners yet, and this represents a serious limitation of these methods. Our preconditioner allows us to perform $s$ iterations of the iterative method with no communication, through ghosting some of the input data and performing redundant computation. To avoid communication, an alternating reordering algorithm is introduced for structured and well partitioned unstructured matrices, which requires the input matrix to be ordered by using a graph partitioning technique such as k-way or nested dissection. We show that the reordering does not affect the convergence rate of the ILU0 preconditioned system as compared to k-way or nested dissection ordering, while it reduces data movement and is expected to reduce the time needed to solve a linear system. In addition to communication avoiding Krylov subspace methods, our preconditioner can be used with classical methods such as GMRES to reduce communication.