Lower-upper Factorization Research Articles

NUMA-aware parallel sparse LU factorization for SPICE-based circuit simulators on ARM multi-core processors

In circuit simulators that resemble the Simulation Program with Integrated Circuit Emphasis (SPICE), one of the most crucial steps is the solution of numerous sparse linear equations generated by frequency domain analysis or time domain analysis. The sparse direct solvers based on lower-upper (LU) factorization are extremely time-consuming, so their performance has become a significant bottleneck. Despite the existence of some parallel sparse direct solvers for circuit simulation problems, they remain challenging to adapt in terms of performance and scalability in the face of rapidly evolving parallel computers with multiple NUMA hardware based on ARM architecture. In this paper, we introduce a parallel sparse direct solver named HLU, which re-examines the performance of the parallel algorithm from the viewpoint of parallelism in pipeline mode and the computing efficiency of each task. To maximize task-level parallelism and further minimize the thread waiting time, HLU devises a fine-grained scheduling method based on an elimination tree in pipeline mode, which employs depth-first search (DFS-like) to iteratively search for parent tasks and then place dependent tasks in the same task queue. HLU also suggests two NUMA node affinity strategies: thread affinity optimization based on NUMA nodes topology to guarantee computational load balancing and data affinity optimization to enable effective memory placement when threads access data. The rationality and effectiveness of the sparse solver HLU are validated by the SuiteSparse Matrix Collection. In comparison with KLU and NICSLU, the experimental results and analysis show that HLU attains a speedup of up to 9.14× and 1.26x (geometric mean) on a Huawei Kunpeng 920 Server, respectively.

Unique Symbolic Factorization for Fast Contingency Analysis Using Full Newton–Raphson Method

Contingency analysis plays an important role in assessing the static security of a network. Its purpose is to check whether a system can operate safely when some elements are out of service. In a real-time application, the computational time required to perform the calculation is paramount for operators to take immediate actions to prevent cascading outages. Therefore, the numerical performance of the contingency analysis is the main focus of this current research. In power flow calculation, when solving the network equations with a sparse matrix solver, most of the time is spent factorizing the Jacobian matrix. In terms of computation time, the symbolic factorization is the costliest operation in the LU (Lower-upper) factorization process. This paper proposes a novel method to perform the calculation with only one symbolic factorization using a full Newton–Raphson-based generic formulation and modular approach (GFMA). The symbolic factorization retained can be used during the iterations of any power flow contingency scenario. A computer study demonstrates that reusing the same symbolic factorization greatly reduces computation time and improves numerical performance. Power system security assessment under N-1 and N-2 contingency conditions is performed for the IEEE standard 54-bus and 108-bus to evaluate the numerical performance of the proposed method. A comparison with the conventional power flow method shows that the time required for the analysis is shortened considerably, with a minimum gain of 228%. The comparative analysis demonstrates that the proposed solution has better numerical performance for large-scale networks.

Orthogonal Time Frequency Space Detection via Low-Complexity Expectation Propagation

Orthogonal time frequency space (OTFS) can provide better error performance than orthogonal frequency division modulation in the high-speed scenario. However, the two-dimensional convolution between the information-bearing symbols and the channel response in the delay-Doppler domain induces high-complexity detection, which hinders the implementation of the OTFS. In this paper, expectation propagation (EP) is introduced in the OTFS system to improve the reliability compared to the message passing approaches. The low-complexity EP detection with a log-linear order is proposed by observing the sparsity and the block quasi-banded matrix structure of the time-domain matrix. The exploitation of the lower-upper factorization for the banded matrix, the forward or the back substitution algorithm for the lower or the upper triangular matrices reduces the complexity of the large-scale matrix inversion involved in the EP detection from a cubic order to a linear order. The error performance of the proposed scheme is further analyzed by utilizing the state evolution and discussed under the different frame sizes. Moreover, the proposed low-complexity detector is then explored together with the likelihood decoder in an iterative manner to further improve the reliability of the proposed receiver. Finally, the better error performance and the lower computational complexity of the proposed receiver are validated by the numerical results.

A fully coupled, fully implicit simulation method for unsteady flames using Jacobian approximation and clustering

Jacobian clustering is proposed to reduce the computational cost associated with matrix operations encountered in the Newton iteration in fully coupled, fully implicit schemes for unsteady reactive flow simulations with detailed chemistry. The iterative solver is based on the Lower-Upper Symmetric Gauss-Seidel (LUSGS) algorithm and sparse matrix technique. The evaluation and sparse Lower-Upper (LU) factorization of the diagonal block of the system Jacobian are performed within clusters rather than individual cells for Computational Fluid Dynamics (CFD) simulations. Cells that are close in the state space are clustered to provide the averaged states for calculating the Jacobians of chemical source terms and transport fluxes. The cells retrieve the factorized sparse matrices from the belonging clusters to perform the necessary iterations. For the purpose of clustering, the spatial dependency of transport Jacobian in the diagonal block is eliminated. To further reduce the computational cost, the sparsity of chemical Jacobian is augmented by removing the insignificant matrix elements. The method is tested in various one-dimensional hydrocarbon flames with both the second order Crank-Nicolson scheme and a third order implicit Runge-Kutta scheme. Various chemical mechanisms with 9 to 111 species are used to test the performance of the iterative solver. Fast convergence of Newton iteration is achieved, and the formal order of accuracy is demonstrated with Jacobian clustering. The overall costs of evaluation and factorization of the block diagonal Jacobian are negligible compared to the cost of calculating transport fluxes and chemical source terms. The averaged costs of Jacobian evaluation, LU factorization and Newton iteration, all increase only linearly with the number of chemical species. The fully coupled, fully implicit Crank-Nicolson scheme with Jacobian clustering shows 4 to 42 times speedup in computational time compared to the decoupled implicit scheme with Strang operator splitting. Jacobian clustering is promising to increase the computational efficiency of high order fully coupled, fully implicit schemes for unsteady reactive flow simulations with detailed chemistry.

Approximation of P-, S1-, and S2-wave reflection coefficients for orthorhombic media

Orthorhombic (ORT) media are commonly used in the descriptions of subsurface strata during the exploration and exploitation of unconventional hydrocarbon reservoirs. In anisotropic media, seismic waves propagating in ORT media usually occur as quasi-P, quasi-S1, and quasi-S2 waves. Approximating the phase velocities and reflection coefficients of these waves is of great significance for predicting the properties of unconventional reservoirs, such as their elastic parameters, fracture azimuths, and fracture densities. At the incidence angles of P waves, we have developed formulas to approximate the reflection coefficients of quasi-P, quasi-S1, and quasi-S2 waves based on their first-order phase velocities and polarization vectors. Under weak anisotropy and weak impedance contrast assumptions, in the [ X, Z] and [ Y, Z] symmetry planes, our approximations basically match Rüger’s approximations at small-to-moderate incidence angles. Numerical analyses verify that our approximations are accurate at all azimuth angles under high impedance contrast and strong anisotropy. Two important innovations are provided in our manuscript. First, the errors in the proposed reflection coefficients only lie in the first-order phase velocities and polarization vectors. Second, we use lower-upper factorization to decompose matrix-dependent reflection coefficients into multiple simple nested formulas, which can be used to calculate the reflection coefficients without matrix operations.

A novel and efficient engine for P-/S-wave-mode vector decomposition for vertical transverse isotropic elastic reverse time migration

Wave-mode decomposition plays a very important role in elastic reverse time migration (ERTM). Improved imaging quality can be achieved due to reduced wave-mode crosstalk artifacts. The current state-of-the-art methods for anisotropic wavefield separation are based on either splitting model strategy, low-rank approximation, or lower-upper (LU) factorization. Most of these involve expensive matrix computation and Fourier transforms with strong model assumptions. Based on the anisotropic-Helmholtz (ani-Helmholtz) decomposition operator and decoupled formulations, we develop a novel and efficient P-/S-wave-mode vector decomposition method in vertical transverse isotropic (VTI) media with application in ERTM. We first review the basics of classical Helmholtz decomposition and isotropic decoupled formulations. In addition, the ani-Helmholtz decomposition operator is built from the P- and S-wave polarizations of the Christoffel equation in VTI media. We then derive novel decoupled formulations of anisotropic P-/S-waves based on the obtained ani-Helmholtz operator. Moreover, we use the first-order Taylor expansion to approximate the normalization term from the derived decoupled formulations and obtain an efficient ani-Helmholtz decomposition approach, which produces vector P- and S-wavefields with correct units, phases, and amplitudes. Compared with the previous studies, our approach mitigates model assumptions, avoids intricate calculations, such as LU factorization or low-rank approximation, and only needs three fast Fourier transforms at each time step. In addition, the graphic processing unit technique is used to dramatically accelerate various functions of ERTM, such as wavefields extrapolation, decomposition, and imaging. Three synthetic examples demonstrate the effectiveness and feasibility of our proposed approach.

Accurate 3D frequency-domain seismic wave modeling with the wavelength-adaptive 27-point finite-difference stencil: A tool for full-waveform inversion

Efficient frequency-domain full-waveform inversion (FWI) of long-offset node data can be performed with a few frequencies. The seismic response of these frequencies can be computed with compact finite-difference stencils on regular Cartesian grid with direct or hybrid direct/iterative methods. Compactness, which is necessary to mitigate the fill-in induced by the lower-upper factorization, is implemented with second-order stencils, whereas accuracy is achieved by building a consistent mass matrix and a compound stiffness matrix with optimal weights so that the stencil covers the eight cells surrounding the central point, leading to a 27-point stencil. Classical approaches estimate constant weights by jointly minimizing numerical dispersion in homogeneous media for several numbers of grid points per wavelength ([Formula: see text]). Then, the impedance matrix is built by using these constant weights at each central point covering the heterogeneous subsurface model, leading to nonuniform wavefield accuracy. Instead, we estimate [Formula: see text]-dependent weights once and for all by minimizing dispersion separately for each value of [Formula: see text] in a range found in FWI applications. Then, we build the impedance matrix without computational overhead by selecting at each central point the weights corresponding to the local [Formula: see text], hence leading to a wavelength-adaptive stencil. This separate approach with adaptive weights makes wavefield accuracy uniform. We benchmark wavefields, which are computed in several 3D large-scale subsurface models with a sparse multifrontal direct solver and the nonadaptive/adaptive stencils, against analytical solutions when available and the highly accurate discretization-free convergent Born series method. Each benchmark reveals the higher accuracy of the adaptive stencil relative to the nonadaptive one. In the presence of sharp contrasts, the adaptive stencil also is more accurate than a classical [Formula: see text] finite-different time-domain staggered-grid stencil. The relevance of the adaptive stencil is finally illustrated with a 3D FWI case study in the 3.5–13 Hz frequency band.

Numerical Simulation of Multiphase Multicomponent Flow in Porous Media: Efficiency Analysis of Newton-Based Method

Newton’s method has been widely used in simulation multiphase, multicomponent flow in porous media. In addition, to solve systems of linear equations in such problems, the generalized minimal residual method (GMRES) is often used. This paper analyzed the one-dimensional problem of multicomponent fluid flow in a porous medium and solved the system of the algebraic equation with the Newton-GMRES method. We calculated the linear equations with the GMRES, the GMRES with restarts after every m steps—GMRES (m) and preconditioned with Incomplete Lower-Upper factorization, where the factors L and U have the same sparsity pattern as the original matrix—the ILU(0)-GMRES algorithms, respectively, and compared the computation time and convergence. In the course of the research, the influence of the preconditioner and restarts of the GMRES (m) algorithm on the computation time was revealed; in particular, they were able to speed up the program.

Performance Analysis of the χMD Matrix Solver Package for MODFLOW-USG.

The χMD matrix solver package is incorporated into USGS groundwater modeling software, such as MODFLOW-NWT, MODFLOW-USG, and MT3D. The solver is used to solve matrices assembled through numerical discretization of the groundwater flow equation, and solute transport equations. χMD has demonstrated its higher robustness, faster execution speed, and more efficient memory usage compared to the existing solvers for many types of groundwater flow problems. χMD uses preconditioned iterative Krylov-subspace methods and consists of preconditioning and acceleration modules. Because the solver package uses a variety of preconditioning features including level-based incomplete lower-upper (ILU) factorization method with a drop tolerance scheme, users must choose optimal preconditioning parameters to improve execution speed and robustness. In order to examine how the preconditioning parameters, ILU factorization level, and drop tolerance values affect the overall performance of the matrix solver, we evaluated five different groundwater model applications using MODFLOW-USG that include different numerical complexities. For those five cases, the number of discretization nodes varied from 10,000 cells to 730,300 cells. From the analysis, we found that the preconditioning parameters greatly affect execution times and memory usage of the preconditioning and acceleration procedures. In addition, a combination of the ILU level between five to seven and the drop tolerance value between 10-2 and 10-3 usually resulted in shorter overall execution time. Our study suggests that the users can elicit higher performance and robustness of the χMD matrix solver using this combination of the parameters and enhance computational efficiency of solving groundwater and solute transport problems.

A spectral element method to compute approximations of the anisotropic diffusion operator with bidimensional tensor coefficient

We derived a continuous Galerkin spectral element method to compute approximations of the anisotropic diffusion operator with the tensorial coefficient of order two. After writing the operator in the weak form, we calculated spatial integration using the Legendre base. Then, we applied three different methods for the solution of the associated linear system: The Lower-Upper factorization, the biconjugate gradient method and the biconjugate gradient stabilized method with the Richardson preconditioner. To validate the algorithm, we present a convergence study for the Poisson equation when varying the functions of the tensor coefficient. Results show how the algorithm can be used to construct solvers for more complex problems like the anisotropic diffusion-reaction equation.

Explicit nonlinear predictive control algorithms for Laguerre filter and sparse least square support vector machine-based Wiener model

In this paper, three computationally proficient model predictive control (MPC) algorithms for least square support vector machine (LSSVM)-based Wiener model are described. A Wiener model with Laguerre filter as dynamic linear part and LSSVM approximator as nonlinear static part is considered. Even though having excellent approximation abilities, LSSVM suffers from lack of sparseness. A pruning algorithm for LSSVM model is proposed and its comparison is made with classical pruning algorithm. The proposed pruning algorithm is able to remove 99% of support vectors with no remarkable drop in modelling accuracy. Using pruned Wiener model, three computationally efficient MPC algorithms are described. In the first algorithm, linearization of Wiener model is performed at every sampling interval and therefore control vector is determined by carrying out a quadratic optimization task. In the second algorithm, control signal is determined by an explicit control law and parameters of this control law are computed by performing lower-upper (LU) factorization of a matrix and solving linear equations without any online optimization. In the third algorithm, the parameters of explicit control law are calculated directly by another LSSVM approximator, which is trained offline. The advantages and effectiveness of proposed methods are demonstrated on the benchmark pH neutralization reactor. The control performance and computational efficiency of proposed algorithms are compared with computationally complex nonlinear MPC, which repeats a nonlinear optimization task at every sampling interval. The impact of pruning on model accuracy, computational efficiency and control accuracy is also discussed.

GLU3.0: Fast GPU-based Parallel Sparse LU Factorization for Circuit Simulation

Editor's note: Many scientific computing problems, including circuit simulations, rely on efficient lower-upper (LU) decomposition of sparse matrices. Prior studies took advantage of GPUs to parallelize LU decomposition, but they suffer from nontrivial data dependencies. This article presents a new method, called GLU3.0, to accelerate GPU-based sparse LU factorization. -Umit Ogras, Arizona State University.

Online Thevenin Equivalent Parameter Identification Method of Large Power Grids Using LU Factorization

Real-time tracking and identification of Thevenin equivalent parameters is crucial for voltage stability online monitoring of large power grids. An online parameter identification method based on wide-area measurements and nodal analysis is proposed in this paper. The node-voltage equation is utilized directly, and the open-circuit voltage is calculated accurately as the Thevenin equivalent voltage of load buses. Moreover, the equivalent parameters can be identified precisely under the generator reactive power over-limit condition. To avoid repetitive Gaussian elimination of the high-dimensional linear system, lower upper factorization and optimized substitution are used to accelerate the open-circuit voltage solution process. The simulation results of multiple test systems demonstrate the accuracy and rapidity of the proposed method.

Three-Dimensional Wide-Band Electromagnetic Forward Modelling Using Potential Technique

The efficacy of Krylov subspace solvers is strongly dependent on the preconditioner applied to solve the large sparse linear systems of equation for electromagnetic problems. In this study, we present a three-dimensional (3-D) plane wave electromagnetic forward simulation over a broadband frequency range. The Maxwell’s equation is solved in a secondary formulation of the Lorentz gauge coupled-potential technique. A finite-volume scheme is employed for discretizing the system of equations on a structured rectilinear mesh. We employed a block incomplete lower-upper factorization (ILU) preconditioner that is suitable for our potential formulation to enhance the computing time and convergence of the systems of equation by comparing with other preconditioners. Furthermore, we observe their effect on the iterative solvers such as the quasi-minimum residual and bi-conjugate gradient stabilizer. Several applications were used to validate and test the effectiveness of our method. Our scheme shows good agreement with the analytical solution. Notably, from the marine hydrocarbon and the crustal model, the utilisation of the bi-conjugate gradient stabilizer with block ILU preconditioner is the most appropriate. Thus, our approach can be incorporated to optimize the inversion process.

High-Resolution Three-Dimensional Displacement Retrieval of Mining Areas From a Single SAR Amplitude Pair Using the SPIKE Algorithm

High-resolution three-dimensional (3-D) displacements of mining areas are crucial to assess mining-related geohazards and understand the mining deformation mechanism. In 2018, we proposed a cost-effective and robust method for retrieving mining-induced 3-D displacements from a single SAR amplitude pair (SAP) using offset tracking (OT) procedures. Hereafter, we refer to this method as the “alternative OT-SAP” (AOT-SAP) method. A key step in the AOT-SAP method is solving the 3-D surface displacements from the AOT-SAP-constructed linear system using routine lower–upper (LU) factorization. However, if the AOT-SAP method is used to retrieve high-resolution 3-D displacements, the dimension of the linear system becomes very large (in the order millions), and a high-end supercomputer is often needed to perform the LU-based solving procedure. This significantly narrows the practical application of the AOT-SAP method, considering the limited availability of supercomputers. In this paper, owing to the banded nature of the AOT-SAP-constructed linear system, we introduce the SPIKE algorithm as an alternative to LU factorization to solve high-resolution mining-induced 3-D displacements. The SPIKE algorithm is a divide-and-conquer direct solver of a large banded system, which can parallelly or sequentially solve a large banded linear system, with a much smaller memory requirement and a shorter time cost than LU factorization. This allows us to retrieve the high-resolution 3-D mining-induced displacements with the AOT-SAP method on either a supercomputer or a standard personal computer. Finally, the accuracy of the retrieved 3-D displacements and the efficiency improvement of the SPIKE algorithm were tested using both simulation analysis and a real dataset.

Stochastic LU factorizations, Darboux transformations and urn models

Abstract We consider upper‒lower (UL) (and lower‒upper (LU)) factorizations of the one-step transition probability matrix of a random walk with the state space of nonnegative integers, with the condition that both upper and lower triangular matrices in the factorization are also stochastic matrices. We provide conditions on the free parameter of the UL factorization in terms of certain continued fractions such that this stochastic factorization is possible. By inverting the order of the factors (also known as a Darboux transformation) we obtain a new family of random walks where it is possible to state the spectral measures in terms of a Geronimus transformation. We repeat this for the LU factorization but without a free parameter. Finally, we apply our results in two examples; the random walk with constant transition probabilities, and the random walk generated by the Jacobi orthogonal polynomials. In both situations we obtain urn models associated with all the random walks in question.

An Efficient Sixth-Order Newton-Type Method for Solving Nonlinear Systems

In this paper, we present a new sixth-order iterative method for solving nonlinear systems and prove a local convergence result. The new method requires solving five linear systems per iteration. An important feature of the new method is that the LU (lower upper, also called LU factorization) decomposition of the Jacobian matrix is computed only once in each iteration. The computational efficiency index of the new method is compared to that of some known methods. Numerical results are given to show that the convergence behavior of the new method is similar to the existing methods. The new method can be applied to small- and medium-sized nonlinear systems.

An optimized high payload audio watermarking algorithm based on LU-factorization

In the field of audio watermarking, reliably embedding the large number of watermarking bits per second in an audio signal without affecting the audible quality of the host audio with good robustness against signal processing attacks is still one of the most challenging issues. In this paper, a high payload, perceptually transparent and robust audio watermarking solution for such a problem by optimizing the existing problem using genetic algorithm is presented. The genetic algorithm in this paper is used to find the optimal number of audio samples required for hiding each watermarking bit. The embedding is done using the imperceptible properties of LU (lower upper) factorization in wavelet domain. This paper addresses the robustness within perceptual constraints at high payload rate in both mathematical analysis and experimental testing by representing behavior of various attacks using attack characterization. Experimental results show that the proposed audio watermarking algorithm can achieve 1280 bps capacity at an average Signal-to-noise ratio (SNR) of 31.02 dB with good robustness to various signal processing attacks such as noise addition, filtering, and compression. In addition, the proposed watermarking algorithm is blind as it does not require the original signal or watermark during extraction. The comparison of the proposed algorithm with the existing techniques also shows that the proposed algorithm is able to achieve high payload with good robustness under perceptual constraints.

LU factorization on heterogeneous systems: an energy-efficient approach towards high performance

Dense lower–upper (LU) factorization (hereafter referred to as LU) is a critical kernel that is widely used to solve dense linear algebra problems. Hybrid LU algorithms have been well designed to exploit the full capacity of heterogeneous systems. However, existing heterogeneous implementations are typically CPU-centric, which rely highly on CPU cores and suffer from a large amount of data transfers via the PCIe bus, and thus reduce the overall energy efficiency of the entire computer system. In this paper, we provide a coprocessor-resident implementation of LU for a heterogeneous platform to improve energy efficiency by relieving the CPUs from performing heavy load computations and avoiding excessive data transfers via PCIe. To maintain the performance, we conduct optimizations to pipeline the CPU computation, coprocessor computation, MPI communication, and PCIe transfer between the CPUs and coprocessors. The experiments on the Tianhe-2 supercomputer show that our LU implementation can compete with the highly optimized Intel MKL implementation in performance and overcome the limitations of energy efficiency.

GPU-Accelerated Parallel Sparse LU Factorization Method for Fast Circuit Analysis

Lower upper (LU) factorization for sparse matrices is the most important computing step for circuit simulation problems. However, parallelizing LU factorization on the graphic processing units (GPUs) turns out to be a difficult problem due to intrinsic data dependence and irregular memory access, which diminish GPU computing power. In this paper, we propose a new sparse LU solver on GPUs for circuit simulation and more general scientific computing. The new method, which is called GPU accelerated LU factorization (GLU) solver (for GPU LU), is based on a hybrid right-looking LU factorization algorithm for sparse matrices. We show that more concurrency can be exploited in the right-looking method than the left-looking method, which is more popular for circuit analysis, on GPU platforms. At the same time, the GLU also preserves the benefit of column-based left-looking LU method, such as symbolic analysis and column-level concurrency. We show that the resulting new parallel GPU LU solver allows the parallelization of all three loops in the LU factorization on GPUs. While in contrast, the existing GPU-based left-looking LU factorization approach can only allow parallelization of two loops. Experimental results show that the proposed GLU solver can deliver $5.71\times $ and $1.46\times $ speedup over the single-threaded and the 16-threaded PARDISO solvers, respectively, $19.56\times $ speedup over the KLU solver, $47.13\times $ over the UMFPACK solver, and $1.47\times $ speedup over a recently proposed GPU-based left-looking LU solver on the set of typical circuit matrices from the University of Florida (UFL) sparse matrix collection. Furthermore, we also compare the proposed GLU solver on a set of general matrices from the UFL, GLU achieves $6.38\times $ and $1.12\times $ speedup over the single-threaded and the 16-threaded PARDISO solvers, respectively, $39.39\times $ speedup over the KLU solver, $24.04\times $ over the UMFPACK solver, and $2.35\times $ speedup over the same GPU-based left-looking LU solver. In addition, comparison on self-generated $RLC$ mesh networks shows a similar trend, which further validates the advantage of the proposed method over the existing sparse LU solvers.