LU Factorization Research Articles

Balanced incomplete factorization preconditioner with pivoting

In this work we study pivoting strategies for the preconditioner presented in Bru (SIAM J Sci Comput 30(5):2302–2318, 2008) which computes the LU factorization of a matrix A. This preconditioner is based on the Inverse Sherman Morrison (ISM) decomposition [Preconditioning sparse nonsymmetric linear systems with the Sherman–Morrison formula. Bru (SIAM J Sci Comput 25(2):701–715, 2003), that using recursion formulas derived from the Sherman-Morrison formula, obtains the direct and inverse LU factors of a matrix. We present a modification of the ISM decomposition that allows for pivoting, and so the computation of preconditioners for any nonsingular matrix. While the ISM algorithm at a given step computes only a new pair of vectors, the new pivoting algorithm in the k-th step also modifies all the remaining vectors from k+1 to n. Thus, it can be seen as a right looking version of the ISM decomposition. The results of numerical experiments with ill-conditioned and highly indefinite matrices arising from different applications show the robustness of the new algorithm, since it is able to solve problems that are not possible to solve otherwise.

Research on application of LU decomposition and CG method in nonlinear iteration of transmission line finite element method

PurposeThe purpose of this study is to develop a reliable finite element algorithm based on the transmission line method (TLM) to solve the nonlinear problem in electromagnetic field calculation.Design/methodology/approachIn this paper, the TLM has been researched and applied to solve nonlinear iteration in FEM. LU decomposition method and the Jacobi preconditioned conjugate gradient method have been investigated to solve the equations in transmission line finite element method (FEM-TLM). The algorithms have been developed in C++ language. The algorithm is applied to analyze the magnetic field of a long straight current-carrying wire and a single-phase transformer.FindingsFEM-TLM is more effective than traditional FEM in nonlinear iteration. The results of FEM-TLM have been compared and analyzed under different calculation scales. It is found that the LU decomposition method is more suitable for FEM-TLM because there is no need to repeatedly assemble the global coefficient matrix in the iterative solution process and it is not affected by the disturbance of the right-hand vector.Originality/valueAn effective algorithm is provided for solving nonlinear problems in the electromagnetic field, which can save a lot of computing costs. The efficiency of LU decomposition and CG method in FEM-TLM nonlinear iteration is investigated, which also makes a preliminary exploration for the research of FEM-TLM parallel algorithms.

Robust level-3 BLAS Inverse Iteration from the Hessenberg Matrix

Inverse iteration is known to be an effective method for computing eigenvectors corresponding to simple and well-separated eigenvalues. In the non-symmetric case, the solution of shifted Hessenberg systems is a central step. Existing inverse iteration solvers approach the solution of the shifted Hessenberg systems with either RQ or LU factorizations and, once factored, solve the corresponding systems. This approach has limited level-3 BLAS potential since distinct shifts have distinct factorizations. This paper rearranges the RQ approach such that data shared between distinct shifts can be exploited. Thereby the backward substitution with the triangular R factor can be expressed mostly with matrix–matrix multiplications (level-3 BLAS). The resulting algorithm computes eigenvectors in a tiled, overflow-free, and task-parallel fashion. The numerical experiments show that the new algorithm outperforms existing inverse iteration solvers for the computation of both real and complex eigenvectors.

HIFIR: Hybrid Incomplete Factorization with Iterative Refinement for Preconditioning Ill-Conditioned and Singular Systems

We introduce a software package called Hybrid Incomplete Factorization with Iterative Refinement (HIFIR) for preconditioning sparse, unsymmetric, ill-conditioned, and potentially singular systems. HIFIR computes a hybrid incomplete factorization (HIF) , which combines multilevel incomplete LU factorization with a truncated, rank-revealing QR (RRQR) factorization on the final Schur complement. This novel hybridization is based on the new theory of ϵ-accurate approximate generalized inverse (AGI) . It enables near-optimal preconditioners for consistent systems and enables flexible GMRES to solve inconsistent systems when coupled with iterative refinement. In this article, we focus on some practical algorithmic and software issues of HIFIR. In particular, we introduce a new inverse-based rook pivoting (IBRP) into ILU, which improves the robustness and the overall efficiency for some ill-conditioned systems by significantly reducing the size of the final Schur complement for some systems. We also describe the software design of HIFIR in terms of its efficient data structures for supporting rook pivoting in a multilevel setting, its template-based generic programming interfaces for mixed-precision real and complex values in C++, and its user-friendly high-level interfaces in MATLAB and Python. We demonstrate the effectiveness of HIFIR for ill-conditioned or singular systems arising from several applications, including the Helmholtz equation, linear elasticity, stationary incompressible Navier–Stokes (INS) equations, and time-dependent advection-diffusion equation.

Inverting a k-heptadiagonal matrix based on Doolitle LU factorization

Inverting a k-heptadiagonal matrix based on Doolitle LU factorization

Mixed precision low-rank approximations and their application to block low-rank LU factorization

Abstract We introduce a novel approach to exploit mixed precision arithmetic for low-rank approximations. Our approach is based on the observation that singular vectors associated with small singular values can be stored in lower precisions while preserving high accuracy overall. We provide an explicit criterion to determine which level of precision is needed for each singular vector. We apply this approach to block low-rank (BLR) matrices, most of whose off-diagonal blocks have low rank. We propose a new BLR LU factorization algorithm that exploits the mixed precision representation of the blocks. We carry out the rounding error analysis of this algorithm and prove that the use of mixed precision arithmetic does not compromise the numerical stability of the BLR LU factorization. Moreover, our analysis determines which level of precision is needed for each floating-point operation (flop), and therefore guides us toward an implementation that is both robust and efficient. We evaluate the potential of this new algorithm on a range of matrices coming from real-life problems in industrial and academic applications. We show that a large fraction of the entries in the LU factors and flops to perform the BLR LU factorization can be safely switched to lower precisions, leading to significant reductions of the storage and expected time costs, of up to a factor three using fp64, fp32, and bfloat16 arithmetics.

Monte Carlo-transformed field expansion method for simulating electromagnetic wave scattering by multilayered random media.

We present an efficient numerical method for simulating the scattering of electromagnetic fields by a multilayered medium with random interfaces. The elements of this algorithm, the Monte Carlo-transformed field expansion method, are (i) an interfacial problem formulation in terms of impedance-impedance operators, (ii) simulation by a high-order perturbation of surfaces approach (the transformed field expansions method), and (iii) efficient computation of the wave field for each random sample by forward and backward substitutions. Our perturbative formulation permits us to solve a sequence of linear problems featuring an operator that is deterministic, and its LU decomposition matrices can be reused, leading to significant savings in computational effort. With an extensive set of numerical examples, we demonstrate not only the robust and high-order accuracy of our scheme for small to moderate interface deformations, but also how Padé summation can be used to address large deviations.

Iterative Method of Thomas Algorithm on The Case Study of Energy Equation

Implicit method is one of the finite difference method and is widely used for discretization some of ordinary or partial differential equations, such like: advection equation, heat transfer equation, burger equation, and many others. Implicit method is unconditionally stable and has been proved with the approximation of Von-Neumann stability criterion. Actually, implicit method is always identical to block matrices (tri-diagonal matrices or penta-diagonal matrices). These matrices can be solved numerically by Thomas algorithm including Gauss elimination using pivot or not, backward or forward substitution. Furthermore, it can be also solved using LU decomposition method with the elimination of lower triangle matrices first and then the elimination of upper triangle matrices. In this research, Thomas algorithm is used to solve numerically for the problem of convective flow on boundary layer, especially for energy equation with the variation of Prandtl number ( ).

Refined isogeometric analysis of quadratic eigenvalue problems

Certain applications that analyze damping effects require the solution of quadratic eigenvalue problems (QEPs). We use refined isogeometric analysis (rIGA) to solve quadratic eigenproblems. rIGA discretization, while conserving desirable properties of maximum-continuity isogeometric analysis (IGA), reduces the interconnection between degrees of freedom by adding low-continuity basis functions. This connectivity reduction in rIGA’s algebraic system results in faster matrix LU factorizations when using multifrontal direct solvers. We compare computational costs of rIGA versus those of IGA when employing Krylov eigensolvers to solve quadratic eigenproblems arising in 2D vector-valued multifield problems. For large problem sizes, the eigencomputation cost is governed by the cost of LU factorization, followed by costs of several matrix–vector and vector–vector multiplications, which correspond to Krylov projections. We minimize the computational cost by introducing C0 and C1 separators at specific element interfaces for our rIGA generalizations of the curl-conforming Nédélec and divergence-conforming Raviart–Thomas finite elements. Let p be the polynomial degree of basis functions; the LU factorization is up to O((p−1)2) times faster when using rIGA compared to IGA in the asymptotic regime. Thus, rIGA theoretically improves the total eigencomputation cost by O((p−1)2) for sufficiently large problem sizes. Yet, in practical cases of moderate-size eigenproblems, the improvement rate deteriorates as the number of computed eigenvalues increases because of multiple matrix–vector and vector–vector operations. Our numerical tests show that rIGA accelerates the solution of quadratic eigensystems by O(p−1) for moderately sized problems when we seek to compute a reasonable number of eigenvalues.

Two Interpolation Matrix Triangularization Methods for Parametric Level Set-Based Structural Topology Optimization

As an implementation form of basis function, interpolation matrices (IMs) have a crucial impact on parametric level set method (PLSM)-based structural topology optimization (STO). However, there are few studies on compressing IM into triangular matrix (TM) with less storage and computation. Algorithm using LU decomposition and Algorithm using innovative asymmetric basis functions that transform the IMs of compactly supported radial basis functions (CSRBFs) into highly sparse TMs are proposed. Theoretical derivation and numerical experiments show that they effectively improve computational efficiency.

A Krylov accelerated Newton–Raphson scheme for efficient pseudo-arclength pathfollowing

The computational efficiency of an enhanced version of a pseudo-arclength pathfollowing scheme tailored for general multi-degree-of-freedom (multi-dof) nonlinear dynamical systems is discussed. The pathfollowing approach is based on the numerical computation of the Poincaré map and its Jacobian in order to tackle nonautonomous systems with discontinuous vector fields. The scheme is applied to obtain frequency response curves of multi-dof hysteretic systems with a state vector size up to 120, as well as various reduced-order models of single and multiple cantilever beams on a shuttle mass. The proposed approach is shown to drastically increase the speed of convergence in the modified Newton–Raphson scheme thanks to a Krylov sub-space iteration which makes use of the LU decomposition of a frozen Jacobian matrix, which, upon convergence, becomes the monodromy matrix. Several numerical tests performed on mechanical systems with material or geometric nonlinearities corroborate the efficiency of the numerical strategy.The C++ code implementing the proposed methodology is freely available at https://doi.org/10.5281/zenodo.6616482.

An MOR Algorithm Based on the Immittance Zero and Pole Eigenvectors for Fast FEM Simulations of Two-Port Microwave Structures

The aim of this article is to present a novel model-order reduction (MOR) algorithm for fast finite-element frequency-domain simulations of microwave two-port structures. The projection basis used to construct the reduced-order model (ROM) comprises two sets: singular vectors and regular vectors. The first set is composed of the eigenvectors associated with the poles of the finite-element method (FEM) state-space system, while the second one is made up from the eigenvectors corresponding to the zeros of the diagonal elements of the matrix-valued immittance transfer function. Importantly, just one LU factorization of the FEM system is required to construct the projection basis during the reduction process, due to the application of a new formulation based on the Schur complement. The sets of eigenvectors that are used in the basis are independent of one another, which makes the new technique better suited for parallel computing compared with previously developed methods, which are sequential in nature. The reliability and accuracy of the proposed scheme are compared with that of the standard MOR technique, namely, the reduced-basis method (RBM), and verified through the analysis of three microwave structures: an eighth-order dual-mode waveguide filter, a dielectric resonator filter, and a folded waveguide filter.

3D anisotropic TEM modeling with loop source using model reduction method

Abstract For model reduction techniques, there have been relatively few studies performed regarding the forward modeling of anisotropic media in comparison to transient electromagnetic (TEM) forward modeling of isotropic media. The transient electromagnetic method (TEM) responses after the current has been turned off can be represented as a homogeneous ordinary differential equation (ODE) with an initial value, and the ODE can be solved using a matrix exponential function. However, the order of the matrix exponential function is large and solving it directly is challenging, thus this study employs the Shift-and-Invert (SAI-Krylov) subspace algorithm. The SAI-Krylov subspace technique is classified as a single-pole approach compared to the multi-pole rational Krylov subspace approach. It only takes one LU factorization of the coefficient matrix, along with hundreds of backward substitutions. The research in this paper shows that the anisotropic medium has little effect on the optimal shift ${\gamma _{opt}}$ and subspace order m. Furthermore, as compared to the mimetic finite volume method (SAI-MFV) of the SAI-Krylov subspace technique, the method proposed in this paper (SAI-FEM) can further improve the computing efficiency by roughly 13%. In contrast to the standard implicit time step iterative technique, the SAI-FEM method does not require discretization in time, and the TEM response at any moment within the off-time period can be easily computed. Next, the accuracy of the SAI-FEM algorithm was verified by 1D solutions for an anisotropic layer model and a 3D anisotropic model. Finally, the electromagnetic characteristics of the anisotropic anomalous body of the center loop device and separated device of the airborne transient electromagnetic method were analyzed, and it was found that horizontal conductivity has a considerable influence on the TEM response of the anisotropic medium.

Proton-exchange membrane fuel cell ionomer hydration model using finite volume method

In this paper a dynamic membrane electrode assembly water transport model, based on the Finite Volume Method, is presented. The purpose of this paper is to provide an accessible and reproductible model capable of real time simulation. To this aim, a detailed explanation is provided regarding the equations and methods used to compute the physical-based fuel cell model. Additionally, the model is purposely developed using basic code (on Matlab™), to not be limited to a single programming language. Two phase water transport through multi-gaseous porous media (electrodes), interfacial transport, as well as diffusion, convection, and electro-osmosis within the polymer are considered. The main novelty relies in the restructuring of all equations into a single implicit system, which can iteratively be resolved through LU decomposition. This computationally efficient method allows the model to be capable of real-time simulation, by displaying the membrane water content profile evolution on a 3D figure. For nominal PEMFC operating conditions, a dry membrane reaches 35% of its final water concentration value after 2 s, and fully converges after 20 s. The final water content profile displays an 18% gradient (9 and 11 molecules per sulfonic acid sites on the anode and cathode sides, respectively). To calibrate and validate this model, mass transfer (flowmeter) and electrical (ohmmeter) methods have been applied.

Algorithm 1021: SPEX Left LU, Exactly Solving Sparse Linear Systems via a Sparse Left-looking Integer-preserving LU Factorization

SPEX Left LU is a software package for exactly solving unsymmetric sparse linear systems. As a component of the sparse exact (SPEX) software package, SPEX Left LU can be applied to any input matrix, A , whose entries are integral, rational, or decimal, and provides a solution to the system \( Ax = b \) , which is either exact or accurate to user-specified precision. SPEX Left LU preorders the matrix A with a user-specified fill-reducing ordering and computes a left-looking LU factorization with the special property that each operation used to compute the L and U matrices is integral. Notable additional applications of this package include benchmarking the stability and accuracy of state-of-the-art linear solvers and determining whether singular-to-double-precision matrices are indeed singular. Computationally, this article evaluates the impact of several novel pivoting schemes in exact arithmetic, benchmarks the exact iterative solvers within Linbox, and benchmarks the accuracy of MATLAB sparse backslash. Most importantly, it is shown that SPEX Left LU outperforms the exact iterative solvers in run time on easy instances and in stability as the iterative solver fails on a sizeable subset of the tested (both easy and hard) instances. The SPEX Left LU package is written in ANSI C, comes with a MATLAB interface, and is distributed via GitHub, as a component of the SPEX software package, and as a component of SuiteSparse.

Modeling of unsteady flows of multiphase viscous fluid in a pore space

The authors have developed and implemented a numerical algorithm to model unsteady flows of a viscous multiphase isothermal fluid by finite difference method using the projection method for the numerical solution of the Navier-Stokes equation. The projection method implies splitting the initial system of equations by physical processes, in which convective transport and the effect of the pressure gradient are separately taken into account. As a result, at each step, it is necessary to solve the Poisson equation to find the pressure field. The solution of SLAE is performed by a parallel direct solver based on LU decomposition. An explicit scheme is used to solve the Cahn-Hilliard equation to update the phase field, the parameter of which is taken into account when adding surface forces to the Navier-Stokes equation. Computational experiments showing qualitative and quantitative agreement with experimental and numerical data from the literature are presented.

An efficient cascadic multigrid solver for 3-D magnetotelluric forward modelling problems using potentials

SUMMARY The fast and accurate 3-D magnetotelluric (MT) forward modelling is core engine of the interpretation and inversion of MT data. In this study, we develop an improved extrapolation cascadic multigrid method (EXCMG) to solve the large sparse complex linear system arising from the finite-element (FE) discretization on non-uniform orthogonal grids of the Maxwell’s equations using potentials. First, the vector Helmholtz equation and the scalar auxiliary equation are derived from the Maxwell’s equations using Coulomb-gauged potentials. The weighted residual method is adopted to discretize the weak formulation and assemble the FE equation. Secondly, carefully choosing the preconditioned complex stable bi-conjugate gradient method (BiCGStab) as multigrid smoother, we develop an improved EXCMG method on non-uniform grids to solve the resulting large sparse complex non-Hermitian linear systems. Finally, several examples including three standard testing models (COMMEMI3D-1, COMMEMI3D-2 and DTM1.0) and a topographic model are used to validate the accuracy and efficiency of the proposed multigrid solver. Numerical results show that the proposed EXCMG algorithm greatly improves the efficiency of 3-D MT forward modelling, is more efficient than some existing solvers, such as Pardiso, incomplete LU factorization preconditioned biconjugate gradients stabilized method (ILU-BiCGStab) and flexible generalized minimum residual method with auxiliary space Maxwell preconditioner (FGMRES-AMS), and capable to simulate large-scale problems with more than 100 million unknowns.

3D marine controlled-source electromagnetic modeling using an edge-based finite element method with a block Krylov iterative solver

Abstract An edge-based finite element method for the numerical modeling of electromagnetic fields in complex media is presented. We used the analytical solution on an electric field in a homogeneous half space to develop a source correct factor to reduce the influence of source singularity and boundary conditions on the numerical accuracy, so that we can minimize the time required to construct the field source term in the scattered field formula. The modeling domain was discretized using an unstructured tetrahedral mesh. We transformed the complex equations of the electrical field into two real-valued equations by decomposing the field into real and imaginary components. Thereafter, we adopted a block conjugate orthogonal conjugate gradient (BL_COCR) iterative solver with an incomplete LU decomposition preconditioner, which was robust for ill-conditioned systems and efficient for multiple source electromagnetic modeling to solve the real-valued equation systems. Using the analytical solution on an electric field in a homogeneous layer model, we evaluated the accuracy of our numerical forward solution and the results showed that the source correct factor can reduce forward modeling errors associated with boundary effects and source singularities. We also applied the developed algorithm to compute the CSEM responses for typical 3D marine geo-electric models with different number of sources and compared with different iterative solvers, and the results showed that the BL_COCR solver has high computational efficiency when solving multiple right-hand term problems.

GSoFa: Scalable Sparse Symbolic LU Factorization on GPUs

Decomposing a matrix <inline-formula><tex-math notation="LaTeX">$\mathbf {A}$</tex-math></inline-formula> into a lower matrix <inline-formula><tex-math notation="LaTeX">$\mathbf {L}$</tex-math></inline-formula> and an upper matrix <inline-formula><tex-math notation="LaTeX">$\mathbf {U}$</tex-math></inline-formula> , which is also known as LU decomposition, is an essential operation in numerical linear algebra. For a sparse matrix, LU decomposition often introduces more nonzero entries in the <inline-formula><tex-math notation="LaTeX">$\mathbf {L}$</tex-math></inline-formula> and <inline-formula><tex-math notation="LaTeX">$\mathbf {U}$</tex-math></inline-formula> factors than in the original matrix. A <i>symbolic factorization</i> step is needed to identify the nonzero structures of <inline-formula><tex-math notation="LaTeX">$\mathbf {L}$</tex-math></inline-formula> and <inline-formula><tex-math notation="LaTeX">$\mathbf {U}$</tex-math></inline-formula> matrices. Attracted by the enormous potentials of the Graphics Processing Units (GPUs), an array of efforts have surged to deploy various LU factorization steps except for the symbolic factorization, to the best of our knowledge, on GPUs. This article introduces <small>gSoFa</small> , the first <u>G</u> PU-based <u>s</u> ymb <u>o</u> lic <u>fa</u> ctorization design with the following three optimizations to enable scalable LU symbolic factorization for <i>nonsymmetric pattern</i> sparse matrices on GPUs. First, we introduce a novel fine-grained parallel symbolic factorization algorithm that is well suited for the <i>Single Instruction Multiple Thread</i> (SIMT) architecture of GPUs. Second, we tailor supernode detection into a SIMT friendly process and strive to balance the workload, minimize the communication and saturate the GPU computing resources during supernode detection. Third, we introduce a three-pronged optimization to reduce the excessive space consumption problem faced by multi-source concurrent symbolic factorization. Taken together, <small>gSoFa</small> achieves up to 31× speedup from 1 to 44 Summit nodes (6 to 264 GPUs) and outperforms the state-of-the-art CPU project, on average, by 5×. Notably, <small>gSoFa</small> also achieves up to 47 percent of the peak memory throughput of a V100 GPU in the Summit Supercomputer.

Accelerating Restarted GMRES With Mixed Precision Arithmetic

The generalized minimum residual method (GMRES) is a commonly used iterative Krylov solver for sparse, non-symmetric systems of linear equations. Like other iterative solvers, data movement dominates its run time. To improve this performance, we propose running GMRES in reduced precision with key operations remaining in full precision. Additionally, we provide theoretical results linking the convergence of finite precision GMRES with classical Gram-Schmidt with reorthogonalization (CGSR) and its infinite precision counterpart which helps justify the convergence of this method to double-precision accuracy. We tested the mixed-precision approach with a variety of matrices and preconditioners on a GPU-accelerated node. Excluding the incomplete LU factorization without fill in (ILU(0)) preconditioner, we achieved average speedups ranging from 8 to 61 percent relative to comparable double-precision implementations, with the simpler preconditioners achieving the higher speedups.