Sparse Systems Of Linear Equations Research Articles

Minimum residual NDSS iteration method for a class of complex symmetric linear systems

The double-step scale splitting and the minimum residual modified HSS iteration method are effective methods for a class of large sparse complex symmetric systems of linear equations. In this paper, we give an efficient method to solve the linear equations related to complex symmetric systems. Firstly, we modify the double-step scale splitting (DSS) iteration methods, obtaining the new double-step scale splitting (NDSS) iteration methods. Furthermore, we improve the algorithm by taking use of minimum residual technique. We call the modified algorithm minimum residual NDSS (MRNDSS) iteration method. For the new presented method, we study some properties of the algorithm carefully. Compared with the NDSS method, MRNDSS has two different parameters in each alternating iteration step and two additional parameters depending on the residual error from previous step. In our numerical experiments, our new method performs more efficient than four other iteration methods which are usually applied to solve the same problems.

Parallel optimization and application of unstructured sparse triangular solver on new generation of Sunway architecture

Large-scale sparse linear equation solver plays an important role in both numerical simulation and artificial intelligence, and sparse triangular equation solver is a key step in solving sparse linear equation systems. Its parallel optimization can effectively improve the efficiency of solving sparse linear equation systems. In this paper, we design and implement a parallel algorithm for solving sparse triangular equations in combination with the features of the new generation of Sunway architecture, and optimize the access and communication respectively for 949 real equations and 32 complex equations in the SuiteSparse collection. The solution efficiency of the algorithm presented in this paper outperforms the cuSparse algorithm on NVIDIA V100 GPU platforms in more than 71% of the cases, and the speedup is even better in solving larger cases (matrix size greater than 10,000): our method increases the speedup from 1.29 time of the previous version to an average speedup of 5.54 and the best speedup of 32.18 over the sequential method on the next generation of Sunway architecture when using 64 slave cores.

Method for calculating a thermal expansion induced mechanical stress in three-dimensional solid-state structures using mathematical modeling

At the end of the 20th century, the demand for more efficient methods for solving large sparse unstructured linear systems of equations increased dramatically. Classical single-level methods had already reached their limits, and new hierarchical algorithms had to be developed to provide efficient solutions to even larger problems. Efficient numerical solution of large systems of discrete elliptic PDEs requires hierarchical algorithms that provide a fast reduction of both shortwave and longwave components in the error vector expansion. The breakthrough, and certainly one of the most important advances of the last three decades, was due to the multigrid principle. Any appropriate method works with a grid hierarchy specified a priori by coarsening a given sampling grid in a geometrically natural way (a "geometric" multigrid method). However, defining a natural hierarchy can become very difficult for very complex, unstructured meshes, if possible at all. The article proposes an algorithm for calculating the deformation that occurs under the action of a thermal expansion force in three-dimensional solid models based on a grid approximation of the problem by hexagonal 8-node cells. The operation of the algorithm is illustrated by solving three problems.

Automated tuning for the parameters of linear solvers

Robust iterative methods for solving large sparse systems of linear algebraic equations often suffer from the problem of optimizing the corresponding tuning parameters. To improve the performance of the problem of interest, specific parameter tuning is required, which in practice can be a time-consuming and tedious task. This paper proposes an optimization algorithm for tuning the numerical method parameters. The algorithm combines the evolution strategy with the pre-trained neural network used to filter the individuals when constructing the new generation. The proposed coupling of two optimization approaches allows to integrate the adaptivity properties of the evolution strategy with a priori knowledge realized by the neural network. The use of the neural network as a preliminary filter allows for significant weakening of the prediction accuracy requirements and reusing the pre-trained network with a wide range of linear systems.The detailed algorithm efficiency evaluation is performed for a set of model linear systems, including the ones from the SuiteSparse Matrix Collection and the systems from the turbulent flow simulations. The obtained results show that the pre-trained neural network can be effectively reused to optimize parameters for various linear systems, and a significant speedup in the calculations can be achieved at the cost of about 100 trial solves. The hybrid evolution strategy decreases the calculation time by more than 6 times for the black box matrices from the SuiteSparse Matrix Collection and by a factor of 1.4–2 for the sequence of linear systems when modeling turbulent flows. This results in a speedup of up to 1.8 times for the turbulent flow simulations performed in the paper.

Narrow escape in composite domains forming heterogeneous networks

Cellular networks are often composed of thin tubules connecting much larger node compartments. These structures serve for active or diffusion transport of proteins. Examples are glial networks in the brain, the endoplasmic reticulum in cells or dendritic spines. In this latter case, a large ball forming the head is connected to the dendrite by a narrow passage. In all cases, how the transport of molecules, ions or proteins is regulated determines the time scale of chemical reactions or signal transduction. In the present study, based on modeling diffusion in three dimensions, we compute the mean time for a Brownian particle to reach a narrow target inside such a composite network made of tubules connected to spherical nodes. We derive asymptotic formulas by solving a mixed Neumann–Dirichlet boundary value problem with small Dirichlet part. We first consider the general case of a network domain organized in a 2-D lattice structure that consists of spherical ball compartments connected via narrow cylindrical passages. For a single target located on the boundary of one of the spherical domains, we derive a sparse linear system of equations for each Mean First Passage Time (MFPT) averaged over the different compartments. We then consider a composite domain consisting of a spherical head-like domain, with narrow absorbing targets on its boundary, that is connected to a large cylinder via a narrow cylindrical neck. For Brownian particles starting within the neck, we derive an asymptotic formula for the MFPT. When diffusing particles can be absorbed upon hitting additional absorbing boundaries of the large cylindrical compartment, we derive asymptotic formulas for the splitting probability and conditional MFPT to reach a target on the spherical head. We compare these formulas with numerical solutions of the mixed boundary value problem and with Brownian simulations, allowing to explore geometrical parameter ranges. To conclude, the present analysis reveals that the mean arrival time, driven by diffusion in heterogeneous networks, is controlled by the targets and tubules sizes, as well as the size of the nodal compartments.

On multi-step randomized extended Kaczmarz method for solving large sparse inconsistent linear systems

In order to improve the convergence property and computational behavior of the randomized extended Kaczmarz method, we propose a multi-step randomized extended Kaczmarz method, in which we repeatedly update the iterate several times at each iteration step, obtaining a nonstationary inner-outer iteration scheme for solving large-scale, sparse, and inconsistent system of linear equations. For this multi-step randomized extended Kaczmarz method, we prove its convergence, derive an upper bound for its convergence rate, and demonstrate that this upper bound can be smaller than that of the randomized extended Kaczmarz method for several typical choices of the numbers of inner iteration steps. Numerical experiments also show that the multi-step randomized extended Kaczmarz method can perform better than the randomized extended Kaczmarz method if we choose the numbers of inner iteration steps appropriately.

RegCloser: a robust regression approach to closing genome gaps

BackgroundClosing gaps in draft genomes leads to more complete and continuous genome assemblies. The ubiquitous genomic repeats are challenges to the existing gap-closing methods, based on either the k-mer representation by the de Bruijn graph or the overlap-layout-consensus paradigm. Besides, chimeric reads will cause erroneous k-mers in the former and false overlaps of reads in the latter.ResultsWe propose a novel local assembly approach to gap closing, called RegCloser. It represents read coordinates and their overlaps respectively by parameters and observations in a linear regression model. The optimal overlap is searched only in the restricted range consistent with insert sizes. Under this linear regression framework, the local DNA assembly becomes a robust parameter estimation problem. We solved the problem by a customized robust regression procedure that resists the influence of false overlaps by optimizing a convex global Huber loss function. The global optimum is obtained by iteratively solving the sparse system of linear equations. On both simulated and real datasets, RegCloser outperformed other popular methods in accurately resolving the copy number of tandem repeats, and achieved superior completeness and contiguity. Applying RegCloser to a plateau zokor draft genome that had been improved by long reads further increased contig N50 to 3-fold long. We also tested the robust regression approach on layout generation of long reads.ConclusionsRegCloser is a competitive gap-closing tool. The software is available at https://github.com/csh3/RegCloser. The robust regression approach has a prospect to be incorporated into the layout module of long read assemblers.

A two‐step matrix splitting iteration paradigm based on one single splitting for solving systems of linear equations

AbstractFor solving large sparse systems of linear equations, we construct a paradigm of two‐step matrix splitting iteration methods and analyze its convergence property for the nonsingular and the positive‐definite matrix class. This two‐step matrix splitting iteration paradigm adopts only one single splitting of the coefficient matrix, together with several arbitrary iteration parameters. Hence, it can be constructed easily in actual applications, and can also recover a number of representatives of the existing two‐step matrix splitting iteration methods. This result provides systematic treatment for the two‐step matrix splitting iteration methods, establishes rigorous theory for their asymptotic convergence, and enriches algorithmic family of the linear iteration solvers, for the iterative solutions of large sparse linear systems.

An efficient hybrid direct-iterative solver for three-dimensional higher-order edge-based finite element simulation for magnetotelluric data in anisotropic media

This work investigates the use of the high-order finite element method (FEM) for 3D magnetotelluric (MT) modeling in anisotropic media. The high-order FE is more accurate but also more difficult than the linear elements. Numerically, the resulting linear system is larger and denser than the linear case with the same grid. In 3D cases, finite element discretization of Maxwell's equations often leads to large but sparse linear systems of equations, requiring effective methods to solve them. However, for the domain-frequency Maxwell's equations, the discrete rotor operator with a large kernel hampers the convergence of common iterative solvers. Robust and effective solvers are needed to solve 3D electromagnetic (EM) problems. In this study, we presented a robust direct-iterative solver based on a block-diagonal preconditioner. The proposed solver is efficient for high-order FE modeling with arbitrary anisotropy. We implemented the algorithm in a high-level Julia language, which is easy to read, maintain, and extend. Numerical examples show that the high-order elements produce a more accurate solution on a coarse grid than the linear case with the same grid. Compared with the direct solver, the hybrid direct-iterative solver reduces the computation time and memory usage. In addition, the solver only needs a few iterations to converge compared with the traditional BCISGTAB solver with an SSOR preconditioner. Our modeling results suggested the FGMRES solver with the block-diagonal preconditioner is quite effective for 3D MT high-order FE modeling.

SIMD vectorization for simultaneous solution of locally varying linear systems with multiple right-hand sides

Developments in numerical simulation of flows and high-performance computing influence one another. More detailed simulation methods create a permanent need for more computational power, while new hardware developments often require changes to the software to exploit new hardware features. This dependency is very pronounced in the case of vector-units which are featured by all modern processors to increase their numerical throughput but require vectorization of the software to be used efficiently. We study the vectorization of a simulation method that exhibits an inherent level of vector-parallelism. This is of particular interest as SIMD operations will hopefully be available with std::simd in a future C++ standard. The simulation method considered here results in the simultaneous solution of multiple sparse linear systems of equations which only differ by their main diagonal and right-hand sides. Such structure arises in the simulation of unsteady flow in turbomachinery by means of a frequency domain approach called harmonic balance.

Performance Analysis and Optimal Node-aware Communication for Enlarged Conjugate Gradient Methods

Krylov methods are a key way of solving large sparse linear systems of equations but suffer from poor strong scalability on distributed memory machines. This is due to high synchronization costs from large numbers of collective communication calls alongside a low computational workload. Enlarged Krylov methods address this issue by decreasing the total iterations to convergence, an artifact of splitting the initial residual and resulting in operations on block vectors. In this article, we present a performance study of an enlarged Krylov method, Enlarged Conjugate Gradients (ECG), noting the impact of block vectors on parallel performance at scale. Most notably, we observe the increased overhead of point-to-point communication as a result of denser messages in the sparse matrix-block vector multiplication kernel. Additionally, we present models to analyze expected performance of ECG, as well as motivate design decisions. Most importantly, we introduce a new point-to-point communication approach based on node-aware communication techniques that increases efficiency of the method at scale.

Generalized Gearhart-Koshy acceleration for the Kaczmarz method

The Kaczmarz method is an iterative numerical method for solving large and sparse rectangular systems of linear equations. Gearhart, Koshy and Tam have developed an acceleration technique for the Kaczmarz method that minimizes the distance to the desired solution in the direction of a full Kaczmarz step. The present paper generalizes this technique to an acceleration scheme that minimizes the Euclidean norm error over an affine subspace spanned by a number of previous iterates and one additional cycle of the Kaczmarz method. The key challenge is to find a formulation in which all parameters of the least-squares problem defining the unique minimizer are known, and to solve this problem efficiently. When only a single Kaczmarz cycle is considered, the proposed affine search is more effective than the Gearhart-Koshy/Tam line-search, which in turn is more effective than the underlying Kaczmarz method. A numerical experiment from the context of computerized tomography suggests that the proposed affine search has the potential to outperform the the Gearhart-Koshy/Tam line-search and the underlying Kaczmarz method in terms of the computational cost that is needed to achieve a given error tolerance.

Deflation techniques applied on mixed model equations

In this paper, we consider a block Jacobi preconditioner and various deflation techniques applied in the Deflated Preconditioned Conjugate Gradient (DPCG) method for solving a sparse system of linear equations derived from a statistical linear mixed model that analyses simultaneously phenotypic and pedigree information of genotyped and ungenotyped animals with Single Polymorphism Nucleotide genotypes of genotyped animals. In livestock production systems, evaluating the genetic merit of the animals through such a model is a key process to ensure an improvement of animals for some characteristics of interest at each generation. First, we propose to define the deflation vectors using a subdomain deflation approach that considers some biological properties of the genotypes. Using simulated data, this approach reduces the number of iterations by up to 87% in comparison to a Preconditioned Conjugate Gradient method with a Jacobi preconditioner. Furthermore, compared to a DPCG method with same number of subdomains but defined randomly, this approach reduces the number of iterations by up to 20% for the same computational costs of one DPCG iteration. The properties of the resulting systems show that this approach annihilates the largest eigenvalues of the preconditioned coefficient matrix. Second, we propose the use of solution vectors of 12 systems of equations that include between 0.25% and 3% less data, as deflation vectors. For reducing the computational costs, we also consider a Proper Orthogonal Decomposition-reduced set of these 12 vectors. The properties of the resulting systems show that this recycling information approach annihilates the smallest eigenvalues of the preconditioned coefficient matrix, and results in a reduction of up to 39% in comparison to the PCG method. Finally, based on our experiment, the combination of the subdomain deflation approach relying on biological properties and of the POD-based approach to recycle previous solution vectors, for defining the deflation vectors, results in annihilating both the smallest and largest eigenvalues, and in a reduction of up to 88 % of the number of iterations in comparison to the PCG method.

Resolvent analysis of a finite wing in transonic flow

Shock waves interacting with turbulent boundary layers on wings can result first in self-sustained flow unsteadiness and eventually in structural vibration. Due to its importance to modern wing design and aircraft certification, the transonic flow physics continue to be investigated intensively. Herein we focus the discussion on three main aspects. First, we assess a practical implementation of an iterative resolvent algorithm in the linear harmonic incarnation of an industrial computational fluid dynamics code for computing optimal forcing and response modes. This heavily relies on the efficient solution of large sparse linear systems of equations. Second, we showcase its application as a predictive tool to detect transonic buffet flow unsteadiness, well before a global stability analysis can first identify its dynamics through weakly damped eigenmodes, using the NASA common research model at wind-tunnel conditions. Third, we discuss its ability to uncover modal physics, not identifiable through global stability analysis, revealing higher-frequency wake and wingtip vortex modes while shedding some light on the elusive finite wing equivalent of the aerofoil buffet mode. We demonstrate that earlier computational limitations of resolvent analysis, when solving the truncated singular value decomposition using matrix-forming methods with direct matrix factorisation, have been overcome ready for industrial use.

Использование метода декомпозиции области для распараллеливания моделирования течения вязкой несжимаемой среды методом LS-STAG и дополнительного предобуславливания

As a rule, the main part of the computational costs in the numerical solution of problems in continuum mechanics consists in the solving of large sparse systems of linear algebraic equations. For this reason, efficient parallelization of this particular procedure can significantly speed up the simulation. To solve this problem, two main approaches can be used. The simplest approach consists in parallelizing of matrix-vector operations in a usual iterative solver. It requires several synchronization points and exchanges of coefficients at each iteration of the solver, which does not significantly speed up the simulation process as a whole. Therefore, domain decomposition methods are preferable. These methods involve dividing the computational domain into subdomains, constructing and solving separate problems in them, as well as some procedure to coordinate the solution between subdomains to ensure global convergence. Subdomains can overlap, as in the Schwartz method used in OpenFOAM, or they can be separated by interface sections, on which their own interface task is built, as in the Schur complement method. The latter method is used in this research to construct a parallel algorithm for viscous incompressible flow simulation by using the immersed boundary method LS-STAG with cut-cells and level-set functions. The resulting matrix of the interface system has a block tridiagonal structure. To speed up prototyping, OpenMP parallel programming technology is used in the software implementation of the developed algorithm, so computational experiments are carried out only on systems with shared memory, in particular on individual nodes of the educational and experimental cluster of the Applied Mathematics Department, Bauman Moscow State Technical University. To verify and evaluate the effectiveness of the developed parallel algorithm, a well-studied test problem about simulation of two-dimensional flow around a stationary circular airfoil is considered. Computations on a sequence of meshes with different numbers of subdomains show that the parallel algorithm allows one to obtain the same numerical solution as the original algorithm of the LS-STAG method, and the computed values of the Strouhal number and drag coefficient are in good agreement with the experimental and computational data known in the literature. Experiments demonstrate that the developed algorithm with domain decomposition allows to accelerate simulation even in sequential mode by reducing the number of solver iterations, i.e. the domain decomposition method acts as an additional preconditioner. Due to this property, the acceleration is superlinear when simulating in parallel mode with developed algorithm. This effect persists up to a certain number of subdomains, which depends on the size of the problem.

Solution of Large Sparse System of Linear Equations over GF(2) on a Multi Node Multi GPU Platform

We provide an efficient multi-node, multi-GPU implementation of the Block Wiedemann Algorithm (BWA)to find the solution of a large sparse system of linear equations over GF(2). One of the important applications ofsolving such systems arises in most integer factorization algorithms like Number Field Sieve. In this paper, wedescribe how hybrid parallelization can be adapted to speed up the most time-consuming sequence generation stage of BWA. This stage involves generating a sequence of matrix-matrix products and matrix transpose-matrix products where the matrices are very large, highly sparse, and have entries over GF(2). We describe a GPU-accelerated parallel method for the computation of these matrix-matrix products using techniques like row-wise parallel distribution of the first matrix over multi-node multi-GPU platform using MPI and CUDA and word-wise XORing of rows of the second matrix. We also describe the hybrid parallelization of matrix transpose-matrix product computation, where we divide both the matrices row-wise into equal-sized blocks using MPI. Then after a GPU-accelerated matrix transpose-matrix product generation, we combine all those blocks using MPI_BXOR operation in MPI_Reduce to obtain the result. The performance of hybrid parallelization of the sequence generation step on a hybrid cluster using multiple GPUs has been compared with parallelization on only multiple MPI processors. We have used this hybrid parallel sequence generation tool for the benchmarking of an HPC cluster. Detailed timings of the complete solution of number field sieve matrices of RSA-130, RSA-140, and RSA-170 are also compared in this paper using up to 4 NVidia V100 GPUs of a DGX station. We got a speedup of 2.8 after parallelization on 4 V100 GPUs compared to that over 1 GPU.

Iterative solver with folded preconditioner for finite element simulation of magnetotelluric fields

We propose an iterative solver with a folded preconditioner for large sparse systems of linear algebraic equations (SLAEs), which are obtained as a result of a vector finite element approximation of electromagnetic fields in magnetotelluric problems. The solver is based on the computational scheme of the Conjugate A-Orthogonal Conjugate Residual method (COCR), diagonal preconditioner, and residual smoothing. To build a finite element approximation, the primary-secondary field approach is used. To construct a preconditioner, we use a matrix, the elements of which are the weights of the representation of gradients of nodal basis functions as a linear combination of edge basis functions. For Dublin Test Model 1 in a wide frequency range from 0.0001 Hz to 10 Hz, we compare the accuracy and computational efficiency of our code with the results of other authors, as well as compare the computational efficiency of the proposed diagonal folded preconditioned (DFP) COCR solver with other solvers. It is shown that the apparent resistivities obtained on the coarsest mesh with our code differ by no more than a few percent from the solution on a finer mesh and from some solutions obtained by other authors using the finite element, finite difference, and integral equation methods. The computational time required for our code is 8 times less than the minimum time presented by other authors (with a corresponding recalculation of the performance of the computers used). It is also shown that the DFP COCR solver practically does not increase the number of iterations at low frequencies and its computational time is about 20 times less than the time required for the COCR solver with a standard diagonal preconditioner. A comparison of the DFP COCR solver with the Parallel Direct Sparse solver (PARDISO) shows that even at low and medium frequencies, the computational time of the DFP COCR solver is about 20 times less than the computational time of the PARDISO, and at high frequencies this difference reaches already 50 times.

Parallel implementation of the SHYFEM (System of HydrodYnamic Finite Element Modules) model

Abstract. This paper presents the message passing interface (MPI)-based parallelization of the three-dimensional hydrodynamic model SHYFEM (System of HydrodYnamic Finite Element Modules). The original sequential version of the code was parallelized in order to reduce the execution time of high-resolution configurations using state-of-the-art high-performance computing (HPC) systems. A distributed memory approach was used, based on the MPI. Optimized numerical libraries were used to partition the unstructured grid (with a focus on load balancing) and to solve the sparse linear system of equations in parallel in the case of semi-to-fully implicit time stepping. The parallel implementation of the model was validated by comparing the outputs with those obtained from the sequential version. The performance assessment demonstrates a good level of scalability with a realistic configuration used as benchmark.

Smoothness on rank-order path graphs and its use in compressive spectral imaging with side information

This paper proposes a novel reconstruction approach to compressive spectral imaging (CSI) with panchromatic side information, which is based on the notion of approximate rank-order statistics. To that end, we assume that the signal of interest is sufficiently smooth on an unknown graph. When restricted to the family of path graphs, we show that the best path is indeed the rank-order path graph induced by the signal. That is, the path graph whose edge structure is given by the permutation that sorts the entries of the signal in ascending order. Our goal is to show that smoothness on rank-order path graphs inferred from the rank-order statistics of a co-registered panchromatic signal can be used to find accurate spectral image estimates from a compressive snapshot of the scene. We derive theoretical properties of rank-order path graphs and give illustrative examples of their use in signal recovery from undersampled measurements. Our approach leads to solutions with a closed-form, found efficiently by iterative inversion of highly sparse systems of linear equations. We evaluate our method through an experimental demonstration and extensive simulations. Our method performs notably better against a bilateral-filter graph model, adapted to the task, and some traditional and state-of-the-art algorithms.

ParGeMSLR: A parallel multilevel Schur complement low-rank preconditioning and solution package for general sparse matrices

This paper discusses parGeMSLR, a C++/MPI software library for the solution of sparse systems of linear algebraic equations via preconditioned Krylov subspace methods in distributed-memory computing environments. The preconditioner implemented in parGeMSLR is based on algebraic domain decomposition and partitions the symmetrized adjacency graph recursively into several non-overlapping partitions via a p-way vertex separator, where p is an integer multiple of the total number of MPI processes. From a numerical perspective, parGeMSLR builds a Schur complement approximate inverse preconditioner as the sum between the matrix inverse of the interface coupling matrix and a low-rank correction term. To reduce the cost associated with the computation of the approximate inverse matrices, parGeMSLR exploits a multilevel partitioning of the algebraic domain. The parGeMSLR library is implemented on top of the Message Passing Interface and can solve both real and complex linear systems. Furthermore, parGeMSLR can take advantage of hybrid computing environments with in-node access to one or more Graphics Processing Units. Finally, the parallel efficiency (weak and strong scaling) of parGeMSLR is demonstrated on a few model problems arising from discretizations of 3D Partial Differential Equations.