Large Sparse Linear Systems Research Articles

A third order gas-kinetic scheme for unstructured grid

In order to study the high order gas-kinetic scheme on unstructured grid, we combine the third order gas-kinetic flux solver with the compact least-square reconstruction (CLS) method together. The CLS method can achieve arbitrary high order compact reconstruction using the stencil from the whole computational domain implicitly. A large sparse linear system resulted from the CLS reconstruction method is solved by applying the generalized minimal residual algorithm (GMRES). To accelerate the convergence, the Reverse-Cuthill–McKee (RCM) algorithm and the incomplete lower–upper (ILU) factorization method are implemented. Different from the traditional flux solver, the gas-kinetic scheme developed by Xu is of high spatial and temporal accuracy. Applying the second order Taylor expansion of the initial gas distribution function and equilibrium state at cell interface, the flux solver can be extended to third order accuracy directly. The accuracy of present method is validated by several numerical cases such as the advection of density perturbation problem, isotropic vortex propagation problem, Sod shock–wave problem, Lax shock tube test case, Shu–Osher problem, shock–vortex interaction, lid-driven cavity flow, and flat plate boundary layer. The advantages of this high order gas-kinetic scheme are exhibited in some benchmarks including incompressible flow and supersonic compressible flow, inviscid flow and viscous flow.

Robust Preconditioners via Generalized Eigenproblems for Hybrid Sparse Linear Solvers

The solution of large sparse linear systems is one of the most time consuming kernels in many 4 numerical simulations. The domain decomposition community has developed many efficient and robust methods in 5 the last decades. While many of these solvers fall into the abstract Schwarz (aS) framework, their robustness has 6 originally been demonstrated on a case-by-case basis. In this paper, we propose a bound for the condition number 7 of all deflated aS methods provided that the coarse grid consists of the assembly of local components that contain 8 the kernel of some local operators. We show that classical results from the literature on particular instances of 9 aS methods can be retrieved from this bound. We then show that such a coarse grid correction can be explicitly 10 obtained algebraically via generalized eigenproblems, leading to a condition number independent of the number of 11 domains. This result can be readily applied to retrieve or improve the bounds previously obtained via generalized 12 eigenproblems in the particular cases of Neumann-Neumann (NN), Additive Schwarz (AS) and optimized Robin but 13 also generalizes them when applied with approximate local solvers. Interestingly, the proposed methodology turns 14 out to be a comparison of the considered particular aS method with generalized versions of both NN and AS for 15 tackling the lower and upper part of the spectrum, respectively. We furthermore show that the application of the 16 considered grid corrections in an additive fashion is robust in the AS case although it is not robust for aS methods in 17 general. In particular, the proposed framework allows for ensuring the robustness of the AS method applied on the 18 Schur complement (AS/S), either with deflation or additively, and with the freedom of relying on an approximate 19 local Schur complement. Numerical experiments illustrate these statements.

Convergence in Norm of Nonsymmetric Algebraic Multigrid

Algebraic multigrid (AMG) is one of the fastest numerical methods for solving large sparse linear systems. For SPD matrices, convergence of AMG is well motivated in the $A$-norm, and AMG has proven...

A variational reconstructed discontinuous Galerkin method for the steady-state compressible flows on unstructured grids

In this paper, a third-order reconstructed DG(p1p2) method based on the variational reconstruction (VR) (Wang et al., 2017 [24]) is developed for simulating the two dimensional steady-state compressible flows on unstructured and hybrid grids. The proposed method combines the advantages of the DG discretization with the flexibility of the variational reconstruction, which exhibits its superior potential in enhancing the level of accuracy compared to the underlying DG method. In this variational rDG(p1p2) method, the low order degrees of freedom are evolved through the underlying DG(p1) method, while the high order degrees of freedom are reconstructed through the variational reconstruction, in which the constitutive relations are built by minimizing the so-called ‘cost function’. The cost function is defined by the total interfacial jump integration in the computational domain using the variational method. The large sparse linear system resulted by the variational reconstruction is solved in an efficient way coupled with the temporal discretization for the steady-state simulations. A number of test cases are presented to assess the performance of the new high order variational rDG(p1p2) method.

An Operator Smoothing with ILU(0) for Aggregation-based Algebraic Multigrid

Due to the potential optimal convergence, algebraic multigrid is widely used in solving large sparse linear systems. The aggregation-based version is one of the widely used methods for its cheap cost and easy implementation. But its convergence is often slow compared to other versions. In this paper, a smoothing technique based on incomplete LU factorization without fill-in is presented to the operators on each level. Each operator is approximately factorized, and the derived lower and upper triangular factors are approximately inverted, which are applied to the operator from both sides to improve the diagonal dominance, and then the effectiveness of the smoothing process and the accuracy of the grid transfer operators. The numerical results show that when incorporated into the preconditioned conjugate gradient iteration, the convergence rate is greatly improved, and though the time used for setup is larger, the time used for iteration and the overall time can also be reduced for large-scale systems.

A transformation approach that makes SPAI, PSAI and RSAI procedures efficient for large double irregular nonsymmetric sparse linear systems

A sparse matrix is called double irregular sparse if it has at least one relatively dense column and row, and it is double regular sparse if all the columns and rows of it are sparse. The sparse approximate inverse preconditioning procedures SPAI, PSAI( t o l ) and RSAI( t o l ) are costly and even impractical to construct preconditioners for a large sparse nonsymmetric linear system with the coefficient matrix being double irregular sparse, but they are efficient for double regular sparse problems. Double irregular sparse linear systems have a wide range of applications, and 24.4% of the nonsymmetric matrices in the Florida University collection are double irregular sparse. For this class of problems, we propose a transformation approach, which consists of four steps: (i) transform a given double irregular sparse problem into a small number of double regular sparse ones with the same coefficient matrix A ˆ , (ii) use SPAI, PSAI( t o l ) and RSAI( t o l ) to construct sparse approximate inverses M of A ˆ , (iii) solve the preconditioned double regular sparse linear systems by Krylov solvers, and (iv) recover an approximate solution of the original problem with a prescribed accuracy from those of the double regular sparse ones. A number of theoretical and practical issues are considered on the transformation approach. Numerical experiments on a number of real-world problems confirm the very sharp superiority of the transformation approach to the standard approach that preconditions the original double irregular sparse problem by SPAI, PSAI( t o l ) or RSAI( t o l ) and solves the resulting preconditioned system by Krylov solvers.

Direct solutions of 3-D magnetotelluric fields using edge-based finite element

This paper utilizes the 3-array variation of the compressed sparse row (CSR) format to store upper triangular elements of coefficient matrix, and applies the Intel MKL PARDISO solver to solve large sparse linear systems of equations, when calculating the 3-D magnetotelluric (MT) fields using edge-based finite element (FE). By default, Intel MKL PARDISO stores matrices in RAM, namely in-core (IC) mode, however, it is quite memory consumption. One can specify that Intel MKL PARDISO stores matrices on disk, it is referred to as out-of-core (OOC) mode. But the OOC mode consumes longer time than the IC mode. We compare both time and memory consumption of IC and OOC mode via COMMEMI 3D-1A and 3D-2A models. As shown in the examples, the OOC mode can solve very large problems since the amount of main memory required is significantly reduced, in the acceptable computing time, which makes frequency loop parallelization of 3-D MT modeling possible on a small memory PC.

A hybrid parallel algorithm for large sparse linear systems

SummaryLarge sparse linear systems arise in many areas of scientific computing, and the solution of these systems is the most time‐consuming part in many large‐scale problems. We present a hybrid parallel algorithm, named incomplete WZ parallel solver (IWZPS), for the solution of large sparse nonsingular diagonally dominant linear systems on distributed memory architectures. The method is a combination of both direct and iterative techniques. We compare the present hybrid parallel sparse algorithm IWZPS with the direct and iterative sparse solvers, namely, MUMPS and ILUPACK, respectively. In addition, we compare it with a hybrid parallel solver, DDPS.

A modulus iteration method for retinex problem

SummaryRetinex theory was proposed by Land and McCann, which manifests that the perception of object colors by human eyes is only dependent on the reflectance of the object and unrelated to the illumination amount on the object. Image intensity recorded by digital cameras is the product of reflectance and illumination. The main purpose of the retinex problem is to recover the reflectance from the recorded image just like the human vision system. In this paper, we propose a variational minimization model with physical constraints imposed on reflectance values. We show that the proposed model is equivalent to a linear complementarity problem, and a modulus iteration method is applied to solve it. A large sparse linear system of equations arises in the modulus iteration method. By utilizing the special structure of the coefficient matrix, the solution of the linear system is obtained by solving a smaller linear system of only half of the unknowns. The convergence of the modulus iteration method for solving the linear complementarity problem for the proposed model is also demonstrated. The experiments show that the convergence of the proposed method is much faster than the existing efficient methods for the retinex problem, and the proposed method is also competitive to the existing methods for the retinex problem when considering recovered reflectance quality.

A single-step iteration method for non-Hermitian positive definite linear systems

An efficient single-step iteration method is presented for solving the large sparse non-Hermitian positive definite linear systems. We theoretically prove that this method converges to the unique solution of the system of linear equations under suitable restrictions. Moreover, we derive an upper bound for the spectral radius of the new iteration matrix. Furthermore, we consider acceleration of the new iteration by Krylov subspace methods and some special properties of the new preconditioned matrix are proposed. Numerical experiments on a few model problems are presented to further examine the effectiveness of our new method.

An EM-based iterative method for solving large sparse linear systems

ABSTRACTWe propose a novel iterative algorithm for solving a large sparse linear system. The method is based on the EM algorithm. If the system has a unique solution, the algorithm guarantees convergence with a geometric rate. Otherwise, convergence to a minimal Kullback–Leibler divergence point is guaranteed. The algorithm is easy to code and competitive with other iterative algorithms.

Scalable Linear Solvers Based on Enlarged Krylov Subspaces with Dynamic Reduction of Search Directions

Krylov methods are widely used for solving large sparse linear systems of equations. On distributed architectures, their performance is limited by the communication needed at each iteration of the ...

Preconditioned Dirichlet-Dirichlet Methods for Optimal Control of Elliptic PDE

Abstract The discretization of optimal control of elliptic partial differential equations problems yields optimality conditions in the form of large sparse linear systems with block structure. Correspondingly, when the solution method is a Dirichlet-Dirichlet non-overlapping domain decomposition method, we need to solve interface problems which inherit the block structure. It is therefore natural to consider block preconditioners acting on the interface variables for the acceleration of Krylov methods with substructuring preconditioners. In this paper we describe a generic technique which employs a preconditioner block structure based on the fractional Sobolev norms corresponding to the domains of the boundary operators arising in the matrix interface problem, some of which may include a dependence on the control regularization parameter. We illustrate our approach on standard linear elliptic control problems. We present analysis which shows that the resulting iterative method converges independently of the size of the problem. We include numerical results which indicate that performance is also independent of the control regularization parameter and exhibits only a mild dependence on the number of the subdomains.

Using Jacobi iterations and blocking for solving sparse triangular systems in incomplete factorization preconditioning

When using incomplete factorization preconditioners with an iterative method to solve large sparse linear systems, each application of the preconditioner involves solving two sparse triangular systems. These triangular systems are challenging to solve efficiently on computers with high levels of concurrency. On such computers, it has recently been proposed to use Jacobi iterations, which are highly parallel, to approximately solve the triangular systems from incomplete factorizations. The effectiveness of this approach, however, is problem-dependent: the Jacobi iterations may not always converge quickly enough for all problems. Thus, as a necessary and important step to evaluate this approach, we experimentally test the approach on a large number of realistic symmetric positive definite problems. We also show that by using block Jacobi iterations, we can extend the range of problems for which such an approach can be effective. For block Jacobi iterations, it is essential for the blocking to be cognizant of the matrix structure.

A BDDC algorithm for the Stokes problem with weak Galerkin discretizations

The BDDC (balancing domain decomposition by constraints) methods have been applied successfully to solve the large sparse linear algebraic systems arising from conforming finite element discretizations of second order elliptic and Stokes problems. In this paper, the Stokes equations are discretized using the weak Galerkin method, a newly developed nonconforming finite element method. A BDDC algorithm is designed to solve the linear system such obtained. Edge/face velocity interface average and mean subdomain pressure are selected for the coarse problem. The condition number bounds of the BDDC preconditioned operator are analyzed, and the same rate of convergence is obtained as for conforming finite element methods. Numerical experiments are conducted to verify the theoretical results.

On relaxed greedy randomized Kaczmarz methods for solving large sparse linear systems

For solving large sparse systems of linear equations by iteration methods, we further generalize the greedy randomized Kaczmarz method by introducing a relaxation parameter in the involved probability criterion, obtaining a class of relaxed greedy randomized Kaczmarz methods. We prove the convergence of these methods when the linear system is consistent, and show that these methods can be more efficient than the greedy randomized Kaczmarz method if the relaxation parameter is chosen appropriately.

A high order compact least-squares reconstructed discontinuous Galerkin method for the steady-state compressible flows on hybrid grids

In this paper, a class of new high order reconstructed DG (rDG) methods based on the compact least-squares (CLS) reconstruction [23,24] is developed for simulating the two dimensional steady-state compressible flows on hybrid grids. The proposed method combines the advantages of the DG discretization with the flexibility of the compact least-squares reconstruction, which exhibits its superior potential in enhancing the level of accuracy and reducing the computational cost compared to the underlying DG methods with respect to the same number of degrees of freedom. To be specific, a third-order compact least-squares rDG(p1p2) method and a fourth-order compact least-squares rDG(p2p3) method are developed and investigated in this work. In this compact least-squares rDG method, the low order degrees of freedom are evolved through the underlying DG(p1) method and DG(p2) method, respectively, while the high order degrees of freedom are reconstructed through the compact least-squares reconstruction, in which the constitutive relations are built by requiring the reconstructed polynomial and its spatial derivatives on the target cell to conserve the cell averages and the corresponding spatial derivatives on the face-neighboring cells. The large sparse linear system resulted by the compact least-squares reconstruction can be solved relatively efficient when it is coupled with the temporal discretization in the steady-state simulations. A number of test cases are presented to assess the performance of the high order compact least-squares rDG methods, which demonstrates their potential to be an alternative approach for the high order numerical simulations of steady-state compressible flows.

Parallel Schur Complement Techniques Based on Multiprojection Methods

The simulation of several physical phenomena, arising from a wide class of scientific fields, requires solving large sparse linear systems effectively. In order to exploit appropriately the available supercomputing infrastructures, designing efficient and scalable solvers for large sparse linear systems is necessary. A hybrid parallel algebraic iterative method based on the Schur complement method and multiprojection techniques that utilizes semiaggregated subdomains is proposed. The graph corresponding to the coefficient matrix of the sparse linear system is partitioned by a graph partitioning algorithm and the Schur complement problem is solved by a preconditioned Krylov subspace method. The preconditioning scheme is based on domain decomposition and the resulting subdomains consist of fine and aggregated components from the inner boundaries. Moreover, two variants for approximating the local Schur complements, derived from the semiaggregation procedure, are proposed in order to improve memory requirements and preprocessing time. The proposed method is suitable for distributed memory systems with multicore nodes and has improved convergence behavior for a large number of workstations. Numerical and comparative results concerning the applicability, the convergence behavior, and the parallel performance of the proposed scheme are given.

A Max-Plus Approach to Incomplete Cholesky Factorization Preconditioners

We present a new method for constructing incomplete Cholesky factorization preconditioners for use in solving large sparse symmetric positive-definite linear systems. This method uses max-plus algebra to predict the positions of the largest entries in the Cholesky factor and then uses these positions as the sparsity pattern for the preconditioner. Our method builds on the max-plus incomplete LU factorization preconditioner recently proposed in [J. Hook and F. Tisseur, Incomplete LU preconditioner based on max-plus approximation of LU factorization, MIMS Eprint 2016.47, Manchester, 2016] but applied to symmetric positive-definite matrices, which comprise an important special case for the method and its application. A attractive feature of our approach is that the sparsity pattern of each column of the preconditioner can be computed in parallel. Numerical comparisons are made with other incomplete Cholesky factorization preconditioners using problems from a range of practical applications. We demonstrate that the new preconditioner can outperform traditional level-based preconditioners and offer a parallel alternative to a serial limited-memory based approach.

Varying the s in your s-step GMRES

Krylov subspace methods are commonly used iterative methods for solving large sparse linear systems, however they suffer from communication bottlenecks on parallel computers. Therefore, $s$-step methods have been developed where the Krylov subspace is built block by block, so that $s$ matrix-vector multiplications can be done before orthonormalizing the block. Then Communication-Avoiding algorithms can be used for both kernels. This paper introduces a new variation on $s$-step GMRES in order to reduce the number of iterations necessary to ensure convergence, with a small overhead in the number of communications. Namely, we develop a $s$-step GMRES algorithm, where the block size is variable and increases gradually. Our numerical experiments show a good agreement with our analysis of condition numbers and demonstrate the efficiency of our variable $s$-step approach.