Sparse Systems Of Linear Equations Research Articles

Monte Carlo method via a numerical algorithm to solve a parabolic problem

This paper is intended to provide a numerical algorithm consisted of the combined use of the finite difference method and Monte Carlo method to solve a one-dimensional parabolic partial differential equation. The numerical algorithm is based on the discretize governing equations by finite difference method. Due to the application of the finite difference method, a large sparse system of linear algebraic equations is obtained. An approach of Monte Carlo method is employed to solve the linear system. Numerical tests are performed in order to show the efficiency and accuracy of the present work.

On inexact hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems

We study theoretical properties of two inexact Hermitian/skew-Hermitian splitting (IHSS) iteration methods for the large sparse non-Hermitian positive definite system of linear equations. In the inner iteration processes, we employ the conjugate gradient (CG) method to solve the linear systems associated with the Hermitian part, and the Lanczos or conjugate gradient for normal equations (CGNE) method to solve the linear systems associated with the skew-Hermitian part, respectively, resulting in IHSS(CG, Lanczos) and IHSS(CG, CGNE) iteration methods, correspondingly. Theoretical analyses show that both IHSS(CG, Lanczos) and IHSS(CG, CGNE) converge unconditionally to the exact solution of the non-Hermitian positive definite linear system. Moreover, their contraction factors and asymptotic convergence rates are dominantly dependent on the spectrum of the Hermitian part, but are less dependent on the spectrum of the skew-Hermitian part, and are independent of the eigenvectors of the matrices involved. Optimal choices of the inner iteration steps in the IHSS(CG, Lanczos) and IHSS(CG, CGNE) iterations are discussed in detail by considering both global convergence speed and overall computation workload, and computational efficiencies of both inexact iterations are analyzed and compared deliberately.

Sheafhom: Software for Sparse Integer Matrices

Sheafhom is a free software package for large-scale computations in the category of finitely-generated modules over the integers and related rings. Its front end is a language for problems in algebraic topology and geometry. These problems come down to sparse systems of linear equations over the integers. Sheafhom's back end solves the sparse systems, with emphasis on avoiding fill-in and integer explosion. We survey and compare algorithms for integer sparse matrices, and we present implementation techniques in Common Lisp. The final section finds the quotient of the free abelian group on 26 letters by sets of words in the dictionary, in the spirit of [8].

Parallel imaging reconstruction for arbitrary trajectories usingk-space sparse matrices (kSPA)

Although the concept of receiving MR signal using multiple coils simultaneously has been known for over two decades, the technique has only recently become clinically available as a result of the development of several effective parallel imaging reconstruction algorithms. Despite the success of these algorithms, it remains a challenge in many applications to rapidly and reliably reconstruct an image from partially-acquired general non-Cartesian k-space data. Such applications include, for example, three-dimensional (3D) imaging, functional MRI (fMRI), perfusion-weighted imaging, and diffusion tensor imaging (DTI), in which a large number of images have to be reconstructed. In this work, a systematic k-space-based reconstruction algorithm based on k-space sparse matrices (kSPA) is introduced. This algorithm formulates the image reconstruction problem as a system of sparse linear equations in k-space. The inversion of this system of equations is achieved by computing a sparse approximate inverse matrix. The algorithm is demonstrated using both simulated and in vivo data, and the resulting image quality is comparable to that of the iterative sensitivity encoding (SENSE) algorithm. The kSPA algorithm is noniterative and the computed sparse approximate inverse can be applied repetitively to reconstruct all subsequent images. This algorithm, therefore, is particularly suitable for the aforementioned applications.

Comparison of iterative pressure solvers for turbulent flow over hills

AbstractThis article discusses the effectiveness of a number of different iterative solvers for the solution of large (sparse) systems of linear algebraic equations. Specific interest is placed upon the behaviour of the methods when applied to the problem of a neutral atmospheric boundary layer encountering a three‐dimensional hill. Here the behaviour of the iterative solver as the flow develops is of primary interest since it can significantly affect the amount of time required to perform a computer simulation. Copyright © 2007 Royal Meteorological Society

A numerical experimental study of inverse preconditioning for the parallel iterative solution to 3D finite element flow equations

Integration of the subsurface flow equation by finite elements (FE) in space and finite differences (FD) in time requires the repeated solution to sparse symmetric positive definite systems of linear equations. Iterative techniques based on preconditioned conjugate gradients (PCG) are one of the most attractive tool to solve the problem on sequential computers. A present challenge is to make PCG attractive in a parallel computing environment as well. To this aim a key factor is the development of an efficient parallel preconditioner. FSAI (factorized sparse approximate inverse) and enlarged FSAI relying on the approximate inverse of the coefficient matrix appears to be a most promising parallel preconditioner. In the present paper PCG using FSAI, diagonal and pARMS (parallel algebraic recursive multilevel solvers) preconditioners is implemented on the IBM SP4/512 and CLX/768 supercomputers with up to 32 processors to solve underground flow problems of a large size. The results show that FSAI may allow for a parallel relative efficiency E p * larger than 50% on the largest problems with p = 32 processors. Moreover, FSAI turns out to be significantly less expensive and more robust than pARMS. Finally, it is shown that E p * for p in the upper range may be much improved if PCG–FSAI is implemented on CLX.

Visually Mapping the RMS Titanic: Conservative Covariance Estimates for SLAM Information Filters

This paper describes a vision-based, large-area, simultaneous localization and mapping (SLAM) algorithm that respects the low-overlap imagery constraints typical of underwater vehicles while exploiting the inertial sensor information that is routinely available on such platforms. We present a novel strategy for efficiently accessing and maintaining consistent covariance bounds within a SLAM information filter, thereby greatly increasing the reliability of data association. The technique is based upon solving a sparse system of linear equations coupled with the application of constant-time Kalman updates. The method is shown to produce consistent covariance estimates suitable for robot planning and data association. Real-world results are reported for a vision-based, six degree of freedom SLAM implementation using data from a recent survey of the wreck of the RMS Titanic.

DC resistivity modelling with a Gaussian quadrature grid

We present a new numerical scheme for 2D/3D direct current resistivity modelling. This method co-operatively combines the solution of the variational principle of the partial differential equation, Gaussian global quadrature abscissae and local cardinal functions so that it has the main advantages of the finite element method and the spectral method. The formulation shows that the method is close to the spectral element method, but it does not require the element mesh or the element integrations, and it makes it much easier to deal with geological models having a 2D/3D complex topography than the traditional numerical methods. It can achieve a similar convergence rate to the spectral element method. We show it transforms the 2D/3D resistivity modelling problem into a sparse and symmetric linear equation system, which can be solved by an iterative or matrix inversion method.

Linear system solution by null‐space approximation and projection (SNAP)

AbstractSolutions of large sparse linear systems of equations are usually obtained iteratively by constructing a smaller dimensional subspace such as a Krylov subspace. The convergence of these methods is sometimes hampered by the presence of small eigenvalues, in which case, some form of deflation can help improve convergence. The method presented in this paper enables the solution to be approximated by focusing the attention directly on the ‘small’ eigenspace (‘singular vector’ space). It is based on embedding the solution of the linear system within the eigenvalue problem (singular value problem) in order to facilitate the direct use of methods such as implicitly restarted Arnoldi or Jacobi–Davidson for the linear system solution. The proposed method, called ‘solution by null‐space approximation and projection’ (SNAP), differs from other similar approaches in that it converts the non‐homogeneous system into a homogeneous one by constructing an annihilator of the right‐hand side. The solution then lies in the null space of the resulting matrix. We examine the construction of a sequence of approximate null spaces using a Jacobi–Davidson style singular value decomposition method, called restarted SNAP‐JD, from which an approximate solution can be obtained. Relevant theory is discussed and the method is illustrated by numerical examples where SNAP is compared with both GMRES and GMRES‐IR. Copyright © 2006 John Wiley & Sons, Ltd.

Deconvolution‐interpolation gridding (DING): Accurate reconstruction for arbitrary k‐space trajectories

A simple iterative algorithm, termed deconvolution-interpolation gridding (DING), is presented to address the problem of reconstructing images from arbitrarily-sampled k-space. The new algorithm solves a sparse system of linear equations that is equivalent to a deconvolution of the k-space with a small window. The deconvolution operation results in increased reconstruction accuracy without grid subsampling, at some cost to computational load. By avoiding grid oversampling, the new solution saves memory, which is critical for 3D trajectories. The DING algorithm does not require the calculation of a sampling density compensation function, which is often problematic. DING's sparse linear system is inverted efficiently using the conjugate gradient (CG) method. The reconstruction of the gridding system matrix is simple and fast, and no regularization is needed. This feature renders DING suitable for situations where the k-space trajectory is changed often or is not known a priori, such as when patient motion occurs during the scan. DING was compared with conventional gridding and an iterative reconstruction method in computer simulations and in vivo spiral MRI experiments. The results demonstrate a stable performance and reduced root mean square (RMS) error for DING in different k-space trajectories.

Static load balancing applied to Schur complement method

A finite element method often leads to large sparse symmetric and positive definite systems of linear equations. We consider parallel solvers based on the Schur complement method on homogeneous parallel machines with distributed memory. A finite element mesh is partitioned by graph partitioning. Such partitioning results in submeshes with similar numbers of elements and, consequently, submatrices of similar sizes. The submatrices are partially factorised. The time spent on the partial factorisation can be different, i.e., disbalanced, because methods exploiting the sparsity of submatrices are used. This paper proposes a Quality Balancing heuristic that modifies classic mesh partitioning so that the partial factorisation times are balanced, which saves overall computation time, especially for time dependent mechanical and nonstationary transport problems.

Regularized conjugate gradient method for skew-symmetric indefinite system of linear equations and applications

A class of regularized conjugate gradient methods is presented for solving large-scale sparse system of linear equations of which the coefficient matrix is an ill-conditioned skew-symmetric indefinite matrix. The convergence is proved and the possible choices of the parameters involved in the new methods are discussed in detail. Preliminary numerical computations show that the numerical behaviors of the new methods are superior to those of some standard Krylov subspace methods, such as CGNE, CGS, GMRES etc.

Parallel ART for image reconstruction in CT using processor arrays

Algebraic Reconstruction Technique (ART) is a widely-used iterative method for solving sparse systems of linear equations. This method (originally due to Kaczmarz) is inherently sequential according to its mathematical definition since, at each step, the current iterate is projected toward one of the hyperplanes defined by the equations. The main advantages of ART are its robustness, its cyclic convergence on inconsistent systems, and its relatively good initial convergence. ART is widely used as an iterative solution to the problem of image reconstruction from projections in computerized tomography (CT), where its implementation with a small relaxation parameter produces excellent results. It is shown that for this particular problem, ART can be implemented in parallel on a linear processor array. Reconstructing an image of n pixels from Θ(n) equations can be done on a linear array of processors with optimal efficiency (linear speedup) and O(n/p) memory for each processor. The parallel technique can be applied to various geometric models of image reconstruction, as well as to 3D reconstruction with spherically symmetric volume elements, using a 2D rectangular mesh-connected array of processors.

New subclasses of block H-matrices with applications to parallel decomposition-type relaxation methods

The parallel decomposition-type relaxation methods for solving large sparse systems of linear equations on SIMD multiprocessor systems have been proposed in [3] and [2]. In case when the coefficient matrix of the linear system is a block \(H\)-matrix, sufficient conditions for the convergence of methods given in [2], [3] have been further improved in [5] and [4]. From the practical point of view, the convergence area obtained there is not always suitable for computation, so we propose new, easily computable ones, for some special subclasses of block \(H\)-matrices. Furthermore, this approach improves the already known convergence area for the class of block strictly diagonally dominant (SDD) matrices.

The SSOR‐preconditioned inner outer flexible GMRES method for the FEM analysis of EM problems

AbstractWhen solving electromagnetic problems by finite element method, the resulting linear systems are ill‐conditioned when the PML absorbing layers are adopted. Traditional Krylov iterative algorithms usually converge slowly or may even stagnate. In this paper, the combination of the inner–outer Flexible Generalized Minimum Residual Method (GMRES) and the symmetric successive overrelaxation (SSOR) preconditioning scheme is used to accelerate the iteration for solving a large sparse system of linear equations resulting from the analysis of EM problems by use of edge‐based finite element method (FEM). The numerical results show that tens of times convergence improvement can be achieved compared with the conventional GMRES solver. © 2006 Wiley Periodicals, Inc. Microwave Opt Technol Lett 48: 1708–1712, 2006; Published online in Wiley InterScience (www.interscience.wiley.com). DOI 10.1002/mop.21758

Equivalent Conditions for Convergence of Splittings of Non-Hermitian Indefinite Matrices

We present some sufficient and necessary conditions for convergent splitting of a non-Hermitian indefinite matrix. Some sufficient conditions to determinate a matrix with a (strongly) dominant symmetric part for a class of boundary value problem are also obtained. These results are applicable to identify the convergence of iterative methods for solving large sparse systems of linear equations.

A frontal solver for the 21st century

SUMMARY In recent years there have been a number of important developments in frontal algorithms for solving the large sparse linear systems of equations that arise fromnite-element problems. We report on the design of a new fully portable and ecient frontal solver for large-scale real and complex unsymmetric linear systems fromnite-element problems that incorporates these developments. The new package oers both aexible reverse communication interface and a simple to use all-in-one interface, which is designed to make the package more accessible to new users. Other key features include automatic element ordering using a state-of-the-art hybrid multilevel spectral algorithm, minimal main memory requirements, the use of high-level BLAS, and facilities to allow the solver to be used as part of a parallel multiple front solver. The performance of the new solver, which is written in Fortran 95, is illustrated using a range of problems from practical applications. The solver is available as package HSL MA42 ELEMENT within the HSL mathematical software library and, for element problems, supersedes the well-known MA42 package. Copyright ? 2006 John Wiley & Sons, Ltd.

Performance Aspects of Sparse Matrix-Vector Multiplication

Sparse matrix-vector multiplication (shortly SpM×V) is an important building block in algorithms solving sparse systems of linear equations, e.g., FEM. Due to matrix sparsity, the memory access patterns are irregular and utilization of the cache can suffer from low spatial or temporal locality. Approaches to improve the performance of SpM×V are based on matrix reordering and register blocking [1, 2], sometimes combined with software-pipelining [3]. Due to its overhead, register blocking achieves good speedups only for a large number of executions of SpM×V with the same matrix A.We have investigated the impact of two simple SW transformation techniques (software-pipelining and loop unrolling) on the performance of SpM×V, and have compared it with several implementation modifications aimed at reducing computational and memory complexity and improving the spatial locality. We investigate performance gains of these modifications on four CPU platforms.

THE IMPROVED PARALLEL ICGS METHOD FOR LARGE AND SPARSE UNSYMMETRIC LINEAR SYSTEMS

For the solutions of large and sparse linear systems of equations with unsymmetric coefficient matrices, we propose an improved version of the Conjugate Gradient Squared method (ICGS) method. The algorithm is derived such that all inner products, matrix-vector multiplications and vector updates of a single iteration step are independent and communication time required for inner product can be overlapped efficiently with computation time of vector updates. Therefore, the cost of global communication on parallel distributed memory computers can be significantly reduced. The resulting ICGS algorithm maintains the favorable properties of the algorithm while not increasing computational costs. Data distribution suitable for both irregularly and regularly structured matrices based on the analysis of the non-zero matrix elements is also presented. Communication scheme is supported by overlapping execution of computation and communication to reduce waiting times. The efficiency of this method is demonstrated by numerical experimental results carried out on a parallel distributed memory system.

Structured preconditioners for nonsingular matrices of block two-by-two structures

For the large sparse block two-by-two real nonsingular matrices, we establish a general framework of practical and efficient structured preconditioners through matrix transformation and matrix approximations. For the specific versions such as modified block Jacobi-type, modified block Gauss-Seidel-type, and modified block unsymmetric (symmetric) Gauss-Seidel-type preconditioners, we precisely describe their concrete expressions and deliberately analyze eigenvalue distributions and positive definiteness of the preconditioned matrices. Also, we show that when these structured preconditioners are employed to precondition the Krylov subspace methods such as GMRES and restarted GMRES, fast and effective iteration solvers can be obtained for the large sparse systems of linear equations with block two-by-two coefficient matrices. In particular, these structured preconditioners can lead to efficient and high-quality preconditioning matrices for some typical matrices from the real-world applications.