Large Systems Of Linear Equations Research Articles

A novel HPL-AI approach for FP16-only accelerator and its instantiation on Kunpeng+Ascend AI-specific platform

HPL-AI, also known as HPL-MxP, is a new benchmark program used to evaluate the upper-bound performance of AI-related tasks on a specific computing cluster. It solves a large linear equation system in FP64, preconditioned by complete LU factorization in lower precision. In this paper, we propose a new HPL-AI approach that relies on the factorization of the coefficient matrix in mixed precision: FP32 diagonals and FP16 off-diagonals. Without compromising the quality of the resultant LU preconditioner, the proposed approach only utilizes the primitive of dense matrix multiplication in FP16 on the accelerator, maximizing the FP16 throughput. Numerical analysis and experiments validate our approach, ensuring avoidance of numerical underflow or overflow during factorization. We implement the proposed approach on Kunpeng+Ascend clusters, a novel AI-specific platform with exceedingly high FP16 peak performance. By applying various optimization techniques, including 2D lookahead, HCCL-based communication pipeline, and SYCL-based tasks overlapping, we achieve 975 TFlops on a single node and nearly 100 PFlops on a cluster of 128 nodes, with a weak scalability of 79.8%.

A study of concurrent multi-frontal solvers for modern massively parallel architectures

Leveraging Trace Theory, we investigate the efficient parallelization of direct solvers for large linear equation systems. Our focus lies on a multi-frontal algorithm, and we present a methodology for achieving near-optimal scheduling on modern massively parallel machines. By employing trace theory with Diekert Graphs and Foata Normal Form, we rigorously validate the correctness of our proposed solution. To establish a strong link between the mesh and elimination tree of the multi-frontal solver, we conduct extensive testing on matrices derived from the Finite Element Method (FEM). Furthermore, we assess the performance of computations on both GPU and CPU platforms, employing practical implementation strategies.

A GPU-based computational framework that bridges neuron simulation and artificial intelligence

Biophysically detailed multi-compartment models are powerful tools to explore computational principles of the brain and also serve as a theoretical framework to generate algorithms for artificial intelligence (AI) systems. However, the expensive computational cost severely limits the applications in both the neuroscience and AI fields. The major bottleneck during simulating detailed compartment models is the ability of a simulator to solve large systems of linear equations. Here, we present a novel Dendritic Hierarchical Scheduling (DHS) method to markedly accelerate such a process. We theoretically prove that the DHS implementation is computationally optimal and accurate. This GPU-based method performs with 2-3 orders of magnitude higher speed than that of the classic serial Hines method in the conventional CPU platform. We build a DeepDendrite framework, which integrates the DHS method and the GPU computing engine of the NEURON simulator and demonstrate applications of DeepDendrite in neuroscience tasks. We investigate how spatial patterns of spine inputs affect neuronal excitability in a detailed human pyramidal neuron model with 25,000 spines. Furthermore, we provide a brief discussion on the potential of DeepDendrite for AI, specifically highlighting its ability to enable the efficient training of biophysically detailed models in typical image classification tasks.

Privacy-Preserving Outsourcing Algorithm for Solving Large Systems of Linear Equations

Privacy-Preserving Outsourcing Algorithm for Solving Large Systems of Linear Equations

Node's residual descent method for linear elastic boundary value problems

This study proposes a node's residual descent method (NRDM) for the linear elasticity boundary value problem. In this method, the boundary information and deformation states are attached to the nodes, and the partial differential equations (PDEs) are then transformed into a generalized difference form with the first-order generalized finite difference (GFD) scheme. Furthermore, this study thoroughly describes the NRDM, confirms the accuracy and convergence by numerical examples, performs an error analysis, and discusses the influence of computational parameters. Results demonstrate that NRDM has high computational accuracy and superliner convergence rate. The double derivation is implemented to solve the second-order PDE boundary value problem with the first-order GFD scheme, thus simplifying the requirement for star connectivity and reducing the computational complexity. NRDM is a matrix-free iterative method that eliminates the task of assembling large global sparse matrices, thus avoiding the numerous matrix–matrix product operations. Rather than solving large systems of linear equations, many extremely small matrix multiplications are present during iterations, which are extremely friendly to parallel computing and can be adapted to the computation of ultra-multi-degree-of-freedom problems. The computational accuracy and convergence rate are considerably improved in NRDM compared to the classical GFD method.

Stein-based preconditioners for weak-constraint 4D-var

Algorithms for data assimilation try to predict the most likely state of a dynamical system by combining information from observations and prior models. Variational approaches, such as the weak-constraint four-dimensional variational data assimilation formulation considered in this problem, can ultimately be interpreted as a minimization problem. One of the main challenges of such a formulation is the solution of large linear systems of equations which arise within the inner linear step of the adopted nonlinear solver. Depending on the adopted approach, these linear algebraic problems amount to either a saddle point linear system or a symmetric positive definite (SPD) one. Both formulations can be solved by means of a Krylov method, like GMRES or CG, that needs to be preconditioned to ensure fast convergence in terms of the number of iterations. In this paper we illustrate novel, efficient preconditioning operators which involve the solution of certain Stein matrix equations. In addition to achieving better computational performance, the latter machinery allows us to derive tighter bounds for the eigenvalue distribution of the preconditioned linear system for certain problem settings. A panel of diverse numerical results displays the effectiveness of the proposed methodology compared to current state-of-the-art approaches.

Quantum radial basis function method for the Poisson equation

The radial basis function (RBF) method is widely used for the numerical solution of the Poisson problem in high dimension, where the approximate solution can be found by solving a large system of linear equations. We demonstrate that the RBF method can be accelerated on a quantum computer by using an efficient quantum algorithm for linear equations. We compare the theoretical performance of our quantum algorithm with that of a standard classical algorithm, and find that the quantum algorithm can achieve a polynomial speedup.

Quadratic programming based partitioning for Block Cimmino with correct value representation

The block Cimmino method is successfully used for the parallel solution of large linear systems of equations due to its amenability to parallel processing. Since the convergence rate of block Cimmino depends on the orthogonality between the row blocks, advanced partitioning methods are used for faster convergence. In this work, we propose a new partitioning method that is superior to the state-of-the-art partitioning method, GRIP, in several ways. Firstly, our proposed method exploits the Mongoose partitioning library which can outperform the state-of-the-art methods by combining the advantages of classical combinatoric methods and continuous quadratic programming formulations. Secondly, the proposed method works on the numerical values in a floating-point format directly without converting them to integer format as in GRIP. This brings an additional advantage of obtaining higher quality partitionings via better representation of numerical values. Furthermore, the preprocessing time is also improved since there is no overhead in converting numerical values to integer format. Finally, we extend the Mongoose library, which originally partitions graphs into only two parts, by using the recursive bisection paradigm to partition graphs into more than two parts. Extensive experiments conducted on both shared and distributed memory architectures demonstrate the effectiveness of the proposed method for solving different types of real-world problems.

Parallel implementations of randomized vector algorithm for solving large systems of linear equations

The results of a parallel implementation of a randomized vector algorithm for solving systems of linear equations are presented in the paper. The solution is represented in the form of a Neumann series. The stochastic method computes this series by sampling only random columns, avoiding multiplication of matrix by matrix and matrix by vector. We consider the case when the matrix is too large to fit in random-access memory (RAM). We use two approaches to solve this problem. In the first approach, the matrix is divided into parts that are distributed among MPI processes and stored in the available RAM of the cluster nodes. In the second approach, the entire matrix is stored on each node’s hard drive, loaded into RAM, and processed in parts. Independent Monte Carlo experiments for random column indices are distributed among MPI processes or OpenMP threads for both approaches to matrix storage. The efficiency of parallel implementations is analyzed. Results are given for a system governed by dense matrices of size $$10^4$$ and $$10^5$$ .

An Aggregation-Based Algebraic Multigrid Method with Deflation Techniques and Modified Generic Factored Approximate Sparse Inverses

In this paper, we examine deflation-based algebraic multigrid methods for solving large systems of linear equations. Aggregation of the unknown terms is applied for coarsening, while deflation techniques are proposed for improving the rate of convergence. More specifically, the V-cycle strategy is adopted, in which, at each iteration, the solution is computed by initially decomposing it utilizing two complementary subspaces. The approximate solution is formed by combining the solution obtained using multigrids and deflation. In order to improve performance and convergence behavior, the proposed scheme was coupled with the Modified Generic Factored Approximate Sparse Inverse preconditioner. Furthermore, a parallel version of the multigrid scheme is proposed for multicore parallel systems, improving the performance of the techniques. Finally, characteristic model problems are solved to demonstrate the applicability of the proposed schemes, while numerical results are given.

Simulation toolkit for digital material characterization of large image-based microstructures

In this paper, an efficient image-based simulation toolkit for material characterization is presented, which is scalable to work from personal computers to workstations. The effective thermal conductivity, elasticity, and permeability are evaluated employing a computational homogenization framework based on the Finite Element Method (FEM). Two complementary open-source packages are presented: one developed in Python, which can convert digital images into voxel meshes (pyTomoviewer); the other developed in Julia, that can run numerical simulations to compute effective material properties (chpack). Also, a CUDA C version of chpack is provided (chfem_gpu). They were designed to deal with large multi-phase models, so strategies were devised to minimize their memory footprint, while avoiding a high toll on execution time. The voxel-based approach significantly simplifies the FEM meshes and allows efficient matrix-free implementations. In that sense, to handle large linear systems of equations, the element-by-element (EBE) technique is adopted, in conjunction with a low-memory implementation of the Preconditioned Conjugate Gradient (PCG) method. The code was thoroughly tested on an artificial geometry made of a square array of cylinders, for which analytical solutions exist, as well as on a real micro-tomographic reconstruction of FiberFormTM, a carbon preform commonly used in thermal protection systems.

Fast computational E-field dosimetry for transcranial magnetic stimulation using adaptive cross approximation and auxiliary dipole method (ACA-ADM)

Transcranial Magnetic Stimulation (TMS) is a non-invasive brain stimulation technique that uses a coil to induce an electric field (E-field) in the brain and modulate its activity. Many applications of TMS call for the repeated execution of E-field solvers to determine the E-field induced in the brain for different coil placements. However, the usage of solvers for these applications remains impractical because each coil placement requires the solution of a large linear system of equations. We develop a fast E-field solver that enables the rapid evaluation of the E-field distribution for a brain region of interest (ROI) for a large number of coil placements, which is achieved in two stages. First, during the pre-processing stage, the mapping between coil placement and brain ROI E-field distribution is approximated from E-field results for a few coil placements. Specifically, we discretize the mapping into a matrix with each column having the ROI E-field samples for a fixed coil placement. This matrix is approximated from a few of its rows and columns using adaptive cross approximation (ACA). The accuracy, efficiency, and applicability of the new ACA approach are determined by comparing its E-field predictions with analytical and standard solvers in spherical and MRI-derived head models. During the second stage, the E-field distribution in the brain ROI from a specific coil placement is determined by the obtained rows and columns in milliseconds. For many applications, only the E-field distribution for a comparatively small ROI is required. For example, the solver can complete the pre-processing stage in approximately 4 hours and determine the ROI E-field in approximately 40 ms for a 100 mm diameter ROI with less than 2% error enabling its use for neuro-navigation and other applications. Highlight: We developed a fast solver for TMS computational E-field dosimetry, which can determine the ROI E-field in approximately 40 ms for a 100 mm diameter ROI with less than 2% error.

Effective Preconditioners for Mixed‐Dimensional Scalar Elliptic Problems

AbstractDiscretization of flow in fractured porous media commonly lead to large systems of linear equations that require dedicated solvers. In this work, we develop an efficient linear solver and its practical implementation for mixed‐dimensional scalar elliptic problems. We design an effective preconditioner based on approximate block factorization and algebraic multigrid techniques. Numerical results on benchmarks with complex fracture structures demonstrate the effectiveness of the proposed linear solver and its robustness with respect to different physical and discretization parameters.

Effective 3D-transient electromagnetic inversion using finite-element method with a parallel direct solver

The correct interpretation of time-domain or transient electromagnetic (TEM) data relies on effective inversion techniques. Herein, we have developed an effective 3D-TEM inversion method based on the vector finite-element method. The domain was discretized with unstructured tetrahedra to simulate complex geologic structures. We develop a second-order backward-Euler scheme, rigorously with nonuniform time-step sizes, to approximate the time derivative of the electric field. The waveform was directly simulated and considered in the forward modeling and inversion. In a 3D inversion, the most time-consuming element of the workflow is solving a large system of linear equations. To alleviate this problem, we adopt a parallel direct solver for solving a large system of equations involved in forward modeling and sensitivity calculations. We use the Gauss-Newton approach to solve the inverse problem and formulate the least-squares problem to obtain model updates. We avoided the explicit sensitivity calculation by solving the matrix-vector multiplication between the sensitivity matrix (or its transpose) and vector. The developed inversion algorithm is validated using several realistic synthetic models. Finally, we apply the method to practical TEM data with complex topography for coal mine goaf detection. We thus obtain a detailed 3D resistivity model that revealed the layering structure and some localized anomalies related to the coal mine goaf.

A method for computing a few eigenpairs of large generalized eigenvalue problems

Many methods for the computation of selected eigenpairs of generalized eigenproblems for matrix pairs use a shift-and-invert technique. When applied to large-scale problems, this requires the solution of large linear systems of equations. This paper proposes an application of an Arnoldi method described in Voss (2004) [21] to the computation of a few extreme eigenpairs of a matrix pair. An advantage of this approach, when compared to methods that use the shift-and-invert technique, is that no large systems of equations have to be solved. We compare this approach to using a technique for simultaneously reducing a pair of large matrices to a pair of small matrices by a generalized Arnoldi process described in Li and Ye (2003) [12] and Hoffnung et al. (2006) [6]. The latter technique does not require the solution of large linear systems of equations either. Computed examples show the proposed method to yield approximations of the desired eigenpairs of higher accuracy when using about the same amount of computer storage space.

Randomized Monte Carlo algorithms for matrix iterations and solving large systems of linear equations

Abstract Randomized scalable vector algorithms for calculation of matrix iterations and solving extremely large linear algebraic equations are developed. Among applications presented in this paper are randomized iterative methods for large linear systems of algebraic equations governed by M-matrices. The crucial idea of the randomized method is that the iterations are performed by sampling random columns only, thus avoiding not only matrix-matrix but also matrix-vector multiplications. The suggested vector randomized methods are highly efficient for solving linear equations of high dimension, the computational cost depends only linearly on the dimension.

MAPPRAISER: A massively parallel map-making framework for multi-kilo pixel CMB experiments

Forthcoming cosmic microwave background (CMB) polarized anisotropy experiments have the potential to revolutionize our understanding of the Universe and fundamental physics. The sought-after, tale-telling signatures will be however distributed over voluminous data sets which these experiments will collect. These data sets will need to be efficiently processed and unwanted contributions due to astrophysical, environmental, and instrumental effects characterized and efficiently mitigated in order to uncover the signatures. This poses a significant challenge to data analysis methods, techniques, and software tools which will not only have to be able to cope with huge volumes of data but to do so with unprecedented precision driven by the demanding science goals posed for the new experiments.A keystone of efficient CMB data analysis is solvers of very large linear systems of equations. Such systems appear in very diverse contexts throughout CMB data analysis pipelines, however they typically display similar algebraic structures and can therefore be solved using similar numerical techniques. Linear systems arising in the so-called map-making problem are one of the most prominent and common ones.In this work we present a massively parallel, flexible and extensible framework, comprised of a numerical library, MIDAPACK, and a high level code, MAPPRAISER, which provide tools for solving efficiently such systems. The framework implements iterative solvers based on conjugate gradient techniques: enlarged and preconditioned using different preconditioners. We demonstrate the framework on simulated examples reflecting basic characteristics of the forthcoming data sets issued by ground-based and satellite-borne instruments, executing it on as many as 16,384 compute cores. The software is developed as an open source project freely available to the community at: https://github.com/B3Dcmb/midapack.

Decision Scheduling for Cloud Computing Tasks Relying on Solving Large Linear Systems of Equations.

With the continuous reform and innovation of Internet technology and the continuous development and progress of social economy, Big Data cloud computing technology is more and more widely used in people's work and life. Many parallel algorithms play a very important role in solving large linear equations in various applications. To this end, this article aims to propose and summarize a cloud computing task scheduling model that relies on the solution of large linear equations. The method of this paper is to study the technology of solving large-scale linear equations and propose an M-QoS-OCCSM scheduling model. The function of the experimental method is to solve the problem of efficiently executing N mutually dependent parallel tasks within limited resources, while fully satisfying users' expectations of task completion time, bandwidth rate, reliability, and cost. In this paper, the application experiment of large-scale linear equations in task scheduling is used to study task scheduling algorithms. The results show that when the task load is 10 and 20, the convergence speed of the MPQGA algorithm is 32 seconds and 95 seconds faster than that of the BGA algorithm, respectively.

An efficient multigrid solver based on a four-color cell-block Gauss-Seidel smoother for 3D magnetotelluric forward modeling

Practical application of 3D magnetotelluric inversion requires efficient forward modeling of electromagnetic (EM) fields in the earth. To resolve realistic 3D structures, large computational domains and extremely large linear systems of equations are required. The iterative solvers, which are almost exclusively used to solve these systems, can be inefficient due to the abundant null space of the curl-curl operator. Multigrid (MG) solvers are considered a potentially efficient technique for solving such problems. However, due to the abundant null solution space and existence of the air layer, MG solvers can still converge slowly or even diverge. We have developed an efficient MG solver for finite-difference frequency-domain EM solution. In this algorithm, the excellent smoothing property of an efficient four-color cell-block Gauss-Seidel (GS) is exploited to remove the short-range errors effectively, and the interpolation and prolongation operators are used to handle the long-range errors. They work as a whole to speed the convergence of our algorithm remarkably. Because all of the nodes for the four-color cell-block GS are grouped into four colors and the edge components attached to different nodes in each color are completely decoupled, this can be used to develop a highly vectorized or parallelized algorithm. Another important property is that our algorithm is locally current divergence free, effectively eliminating spurious solutions in the null space of the curl-curl operator. The accuracy and efficiency of the algorithm are verified by comparing the numerical solutions obtained with our MG solver to those from the biconjugate gradient stabilized solver with different preconditioners based on synthetic models and a model from 3D inversion. Comparisons, in terms of iteration number and computational time, indicate that our algorithm is extremely stable and efficient relative to the other solvers. Our MG algorithm will be suitable for massively parallel computing as well.

Numerical aspects of a Godunov-type stabilization scheme for the Boltzmann transport equation

We discuss the numerical aspects of the Boltzmann transport equation (BE) for electrons in semiconductor devices, which is stabilized by Godunov’s scheme. The k-space is discretized with a grid based on the total energy to suppress spurious diffusion in the stationary case. Band structures of arbitrary shape can be handled. In the stationary case, the discrete BE yields always nonnegative distribution functions and the corresponding system matrix has only eigenvalues with positive real parts (diagonally dominant matrix) resulting in an excellent numerical stability. In the transient case, this property yields an upper limit for the time step ensuring the stability of the CPU-efficient forward Euler scheme and a positive distribution function. Similar to the Monte-Carlo (MC) method, the discrete BE can be solved in time together with the Poisson equation (PE), where the time steps for the PE are split into shorter time steps for the BE, which can be performed at minor additional computational cost. Thus, similar to the MC method, the transient approach is matrix-free and the solution of memory and CPU intensive large systems of linear equations is avoided. The numerical properties of the approach are demonstrated for a silicon nanowire NMOSFET, for which the stationary I–V characteristics, small-signal admittance parameters and the switching behavior are simulated with and without strong scattering. The spurious damping introduced by Godunov’s (upwind) scheme is found to be negligible in the technically relevant frequency range. The inherent asymmetry of the upwind scheme results in an error for very strong scattering that can be alleviated by a finer grid in transport direction.