Sparse Direct Solver Research Articles

Parallel Block Sequential Closed-Form Matting With Fan-Shaped Partitions.

Applying alpha matting to large images is a challenging task because of its computational complexity. This paper provides a divide and conquer strategy for performing closed-form matting. The matting problem, defined for an entire image, is broken down into systems of linear equations defined for very small blocks of the image. The sizes of the small systems are small enough for us to find solutions efficiently using a direct sparse linear equation system solver. The small systems are solved following a sequential order such that the alpha matte grows from a user scribble. With the block sequential application, matting is performed on fan-shaped partitions in parallel on multiple processing cores. Experiments on large test images as well as on standard benchmark test images show that the proposed parallel block sequential matting provides high quality alpha mattes with good scalability.

A Schur complement approach to preconditioning sparse linear least-squares problems with some dense rows

The effectiveness of sparse matrix techniques for directly solving large-scale linear least-squares problems is severely limited if the system matrix A has one or more nearly dense rows. In this paper, we partition the rows of A into sparse rows and dense rows (As and Ad) and apply the Schur complement approach. A potential difficulty is that the reduced normal matrix AsTAs is often rank-deficient, even if A is of full rank. To overcome this, we propose explicitly removing null columns of As and then employing a regularization parameter and using the resulting Cholesky factors as a preconditioner for an iterative solver applied to the symmetric indefinite reduced augmented system. We consider complete factorizations as well as incomplete Cholesky factorizations of the shifted reduced normal matrix. Numerical experiments are performed on a range of large least-squares problems arising from practical applications. These demonstrate the effectiveness of the proposed approach when combined with either a sparse parallel direct solver or a robust incomplete Cholesky factorization algorithm.

Modeling kV X-ray-Induced Coloration in Radiochromic Films

Commercially available radiochromic films are primarily designed for clinical X-ray dosimetry. These films change color upon exposures to radiation as a result of solid-state polymerization (SSP). Built on a previous model developed for SSP upon exposures to ultraviolet (UV) radiation, a new model was developed in the present work for X-ray-induced coloration in Gafchromic EBT3 films. Monte Carlo simulations using the Monte Carlo N-Particle (MCNP) code were employed to model the transport and interaction of photons and the generated secondary electrons within the film active layer. The films were exposed to continuous-energy photon beams. The dose DE in the external radiation detector (i.e., ionization chamber) was determined and the realistic dose DA in the film active layer was then obtained using the calibration coefficient R (=DA/DE). The finite element method (FEM) was used to solve the classical steady-state Helmholtz equation using the multifrontal massively parallel sparse direct solver (MUMPS). An extensive grid independence test was carried out and the numerical stability of the present model was ensured. The reflected light intensity from the film surface was used to theoretically obtain the net reflective optical density of the film exposed to X-ray. Good agreement was obtained between the experimental and theoretical results of the net reflective optical density of the film. For X-ray doses >~600 cGy, due to the already formed densely cross-linked structure in the active layer of the EBT3 film so further bond formation was less likely, the reflected light intensity from the film surface increased at a relatively lower rate when compared to those for dose values <~600 cGy.

Fourth-Order Time-Stepping For Stiff PDEs On The Sphere

We present in this paper algorithms for solving stiff PDEs on the unit sphere with spectral accuracy in space and fourth-order accuracy in time. These are based on a variant of the double Fourier sphere method in coefficient space with multiplication matrices that differ from the usual ones, and implicit-explicit time-stepping schemes. Operating in coefficient space with these new matrices allows one to use a sparse direct solver, avoids the coordinate singularity and maintains smoothness at the poles, while implicit-explicit schemes circumvent severe restrictions on the time-steps due to stiffness. A comparison is made against exponential integrators and it is found that implicit-explicit schemes perform best. Implementations in MATLAB and Chebfun make it possible to compute the solution of many PDEs to high accuracy in a very convenient fashion.

Selected inversion as key to a stable Langevin evolution across the QCD phase boundary

We present new results of full QCD at nonzero chemical potential. In PRD 92, 094516 (2015) the complex Langevin method was shown to break down when the inverse coupling decreases and enters the transition region from the deconfined to the confined phase. We found that the stochastic technique used to estimate the drift term can be very unstable for indefinite matrices. This may be avoided by using the full inverse of the Dirac operator, which is, however, too costly for four-dimensional lattices. The major breakthrough in this work was achieved by realizing that the inverse elements necessary for the drift term can be computed efficiently using the selected inversion technique provided by the parallel sparse direct solver package PARDISO. In our new study we show that no breakdown of the complex Langevin method is encountered and that simulations can be performed across the phase boundary.

Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths

In this paper new asymptotic-preserving schemes are introduced for the iterative resolution of anisotropic elliptic equations arising in magnetized plasma simulations. These methods overcome the resolution of the saddle point problems systematically implemented in precedent realizations. This allows us to easily derive iterative solvers for these asymptotic-preserving schemes. It brings a leap forward in the computational efficiency of the method for three-dimensional problems. This is conclusively outlined thanks to three-dimensional serial computations carrying out tens of millions of unknowns. The gains are substantial in terms of memory as well as computational requirements compared to sparse direct solvers, which represent the only alternative successfully operated so far.

A distributed-memory hierarchical solver for general sparse linear systems

Abstract We present a parallel hierarchical solver for general sparse linear systems on distributed-memory machines. For large-scale problems, this fully algebraic algorithm is faster and more memory-efficient than sparse direct solvers because it exploits the low-rank structure of fill-in blocks. Depending on the accuracy of low-rank approximations, the hierarchical solver can be used either as a direct solver or as a preconditioner. The parallel algorithm is based on data decomposition and requires only local communication for updating boundary data on every processor. Moreover, the computation-to-communication ratio of the parallel algorithm is approximately the volume-to-surface-area ratio of the subdomain owned by every processor. We present various numerical results to demonstrate the versatility and scalability of the parallel algorithm.

An assessment of some solvers for saddle point problems emerging from the incompressible Navier–Stokes equations

Efficient incompressible flow simulations, using inf–sup stable pairs of finite element spaces, require the application of efficient solvers for the arising linear saddle point problems. This paper presents an assessment of different solvers: the sparse direct solver UMFPACK, the flexible GMRES (FGMRES) method with different coupled multigrid preconditioners, and FGMRES with Least Squares Commutator (LSC) preconditioners. The assessment is performed for steady-state and time-dependent flows around cylinders in 2d and 3d. Several pairs of inf–sup stable finite element spaces with second order velocity and first order pressure are used. It turns out that for the steady-state problems often FGMRES with an appropriate multigrid preconditioner was the most efficient method on finer grids. For the time-dependent problems, FGMRES with LSC preconditioners that use an inexact iterative solution of the velocity subproblem worked best for smaller time steps.

2.5D marine CSEM modeling in the frequency-domain based on an improved interpolation scheme at receiver positions

In the marine controlled-source electromagnetic (CSEM) survey, the receivers are usually placed at the seafloor. The resistivity contrast between the seawater and seafloor sediments is large, which can cause difficulties in numerical modeling of CSEM fields at receiver locations. In this paper, we present an improved interpolating method for calculating electric and magnetic fields at the seafloor with a resistivity contrast. This method is applied to the 2. 5 dimensional (2. 5D) frequency-domain CSEM modeling with towed transmitters and receivers located at the seafloor. Considering the discontinuity of the normal electric fields, we use the normal current electric density for interpolation. We simulate the 2. 5D marine CSEM responses by the staggered finite-difference (SFD) method with Fourier transform to the strike direction. The final SFD equations are solved by the direct solver MUMPS (MUltifrontal Massively Parallel Sparse direct Solver). To avoid the source singularities, the secondary-field approach is used and the primary fields excited by the electric dipole source can be calculated quasi-analytically for the one-dimensional (1D) layered background model. We focus on interpolating of electric and magnetic fields in the wavenumber domain to the receiver locations at the seafloor interface between the conductive seawater and resistive seafloor formation. The secondary electric and magnetic fields are used for interpolation instead of the total fields for high numerical accuracy. After performing the inverse Fourier transform to the wavenumbers, the electric and magnetic fields in the space domain are obtained. To check the accuracy of our 2. 5D marine CSEM SFD modeling algorithm with the improved receiver interpolating technique, we compare our results with both the 1D analytical results and the adaptive finite element results. The SFD numerical results are approved to be accurate. We also compare the numerical accuracy between our improved interpolation scheme and others, i.e., the conventional linear interpolation and the rigorous interpolation. The proposed interpolation only utilizes the nodes below/above the seafloor interface, and is proved to be much more accurate than the other two interpolating methods used.

Shared-Memory parallelization of consistent particle method for violent wave impact problems

A shared-memory parallelization is implemented to the recently developed Consistent Particle. Method (CPM) for violent wave impact problems. The advantages of this relatively new particle method lie in four key aspects: (1) accurate computation of Laplacian and gradient operators based on Taylor series expansion, alleviating spurious pressure fluctuation and being able to model two-phase flows characterized by large density difference, (2) a thermodynamics-based compressible solver for modelling compressible air that eliminates the need of determining artificial sound speed, (3) seamless coupling of the compressible air solver and incompressible water solver, and (4) parallelization of the numerical model based on Open Multi-Processing (OpenMP) and a parallel direct sparse solver (Pardiso) to significantly improve computational efficiency. Strong and weak scaling analyses of the parallelized CPM are conducted, showing an efficiency speedup of 100 times or more depending on the size of simulated problem. To demonstrate the accuracy of the developed numerical model, three numerical examples are studied including the benchmark study of wave impact on seawall, and our experimental studies of violent water sloshing under rotational excitations and sloshing impact with entrapped air pocket. CPM is shown to accurately capture highly deformed breaking waves and violent wave impact pressure including pressure oscillation induced by air cushion effect.

Modeling Coloration of a Radiochromic Film with Molecular Dynamics-Coupled Finite Element Method

Radiochromic films change color upon exposures to radiation doses as a result of solid-state polymerization (SSP). Commercially available radiochromic films are primarily designed for, and have become widely used in, clinical X-ray dosimetry. However, many intriguing properties of radiochromic films are not yet fully understood. The present work aimed at developing a theoretical model at both atomic and macroscopic scales to provide a platform for future works to understand these intriguing properties. Despite the fact that radiochromic films were primarily designed for clinical X-ray dosimetry, dose-response curves for the Gafchromic EBT3 film obtained for ultraviolet (UV) radiation were employed to develop our model in order to avoid complications of ionization, non-uniform energy deposition, as well as dispersed doses caused by secondary electrons set in motion by the indirectly ionizing X-ray photons, which might introduce added uncertainties to the model and overshadow the basic SSP processes. The active layer in the EBT3 film consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers, which were modelled using molecular dynamics (MD). The degrees of SSP in the atomic scale upon different UV exposures were obtained to determine the absorption coefficients of the active layer, which were then input into the finite element method (FEM). The classical steady-state Helmholtz equation was engaged to model the reflection from the active layer using the FEM technique. The multifrontal massively parallel sparse direct solver (MUMPS) was employed to solve the present numerical problem. Very good agreement between experimentally and theoretically obtained coloration in terms of net reflective optical density was achieved for different UV exposures. In particular, for UV exposures larger than ~40 J/cm2, the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was explained by the densely cross-linked structure under near-complete polymerization, and thus the lower efficiency for further bond formation between DA monomer strands.

APPLICATION OF RIGID LINKS IN STRUCTURAL DESIGN MODELS

A special finite element modelling rigid links is proposed for the linear static and buckling analysis. Unlike the classical approach based on the theorems of rigid body kinematics, the proposed approach preserves the similarity between the adjacency graph for a sparse matrix and the adjacency graph for nodes of the finite element model, which allows applying sparse direct solvers more effectively. Besides, the proposed approach allows significantly reducing the number of nonzero entries in the factored stiffness matrix in comparison with the classical one, which greatly reduces the duration of the solution. For buckling problems of structures containing rigid bodies, this approach gives correct results. Several examples demonstrate its efficiency.

Robust full‐wave Maxwell solver in time‐domain using magnetic vector potential with edge elements

Magnetic vector potential (MVP) formulations are widely used in low-frequency eddy-current computation. For full-wave Maxwell problems with coupled resistive, inductive and capacitive effects, existing MVP formulations are either ungauged or with too complicated stabilisation procedure or with non-symmetric system matrix. Since sophisticated preconditioners are necessary to ensure good convergence rates of iterative solvers when used to solve the resultant linear system, it is also important to investigate symmetric gauged formulations which can be solved by more robust state-of-the-art sparse direct solvers. Traditional Coulomb gauged formulation using penalty technique is not feasible for edge elements since the divergence of the edge element basis function is zero within each element. In this study, a novel MVP formulation with Coulomb Gauge for full-wave Maxwell problems using edge elements is proposed. The proposed formulation is symmetric and stable for all frequencies, easy to implement and uniquely solvable. Numerical results obtained using the proposed MVP formulation are demonstrated to showcase its stability and accuracy. The method is very promising in further engineering applications.

Block row projection method based on M-matrix splitting

Abstract We propose a hybrid sparse linear system solver based on M-matrix splitting and block-row projection (BRP). We split the sparse coefficient matrix A into two (nonsingular) M-matrices, and construct an augmented larger linear system which we solve using a BRP method. The robustness of BRP is compared with those of ILUT-preconditioned GMRES, and the sparse direct solver Pardiso. We also demonstrate the parallel scalability of BRP on a cluster of multicore nodes.

Numerically Aware Orderings for Sparse Symmetric Indefinite Linear Systems

Sparse symmetric indefinite problems arise in a large number of important application areas; they are often solved through the use of an LDL T factorization via a sparse direct solver. While for many problems prescaling the system matrix A is sufficient to maintain stability of the factorization, for a small but important fraction of problems numerical pivoting is required. Pivoting often incurs a significant overhead, and consequently, a number of techniques have been proposed to try and limit the need for pivoting. In particular, numerically aware ordering algorithms may be used, that is, orderings that depend not only on the sparsity pattern of A but also on the values of its (scaled) entries. Current approaches identify large entries of A and symmetrically permute them onto the subdiagonal, where they can be used as part of a 2 × 2 pivot. This is numerically effective, but the fill in the factor L and hence the runtime of the factorization and subsequent triangular solves may be significantly increased over a standard ordering if no pivoting is required. We present a new algorithm that combines a matching-based approach with a numerically aware nested dissection ordering. Numerical comparisons with current approaches for some tough symmetric indefinite problems are given.

Dimensionless modeling for convective drying of tuberous crop (Solanum tuberosum) by considering shrinkage

The shrinkage of biological materials (agri-food), is a common phenomenon developed during a drying operation that directly affects the product quality. In this work a dimensionless mathematical model including the reduction at the interface of a high moisture content food (tuber) in 1D geometry was developed. The mass and heat conservation equations coupled to the Leibniz-Reynolds Transport Theorem (LRTT) at the interface of a high moisture content food were solved. The system of nonlinear unsteady state equations was solved by using a multi-frontal massively parallel sparse direct Solver. In order to validate this model, tuberous crop slices were dried at three temperatures: 40, 50 and 60 °C, relative humidity of 25% and a constant airflow of 1.5 m/s. Drying kinetics and the evolution of temperature within the product were logged. The model is able to simulate the moisture loss and predict the highest thickness change. The simulations show a good agreement with the experimental drying kinetics, as well as for the evolution of temperature. Simulations show 1-D shrinkage is not linear and the model improves its prediction when the drying temperature is higher (intensive mass flux).

Intel Xeon Phi acceleration of Hybrid Total FETI solver

This paper describes an approach for acceleration of the Hybrid Total FETI (HTFETI) domain decomposition method using the Intel Xeon Phi coprocessors. The HTFETI method is a memory bound algorithm which uses sparse linear BLAS operations with irregular memory access pattern. The presented local Schur complement (LSC) method has regular memory access pattern, that allows the solver to fully utilize the Intel Xeon Phi fast memory bandwidth.This translates to speedup over 10.9 of the HTFETI iterative solver when solving 3 billion unknown heat transfer problem (3D Laplace equation) on almost 400 compute nodes. The comparison is between the CPU computation using sparse data structures (PARDISO sparse direct solver) and the LSC computation on Xeon Phi. In the case of the structural mechanics problem (3D linear elasticity) of size 1 billion DOFs the respective speedup is 3.4.The presented speedups are asymptotic and they are reached for problems requiring high number of iterations (e.g., ill-conditioned problems, transient problems, contact problems). For problems which can be solved with under hundred iterations the local Schur complement method is not optimal. For these cases we have implemented sparse matrix processing using PARDISO also for the Xeon Phi accelerators.

Extending fields in a level set method by solving a biharmonic equation

We present an approach for computing extensions of velocities or other fields in level set methods by solving a biharmonic equation. The approach differs from other commonly used approaches to velocity extension because it deals with the interface fully implicitly through the level set function. No explicit properties of the interface, such as its location or the velocity on the interface, are required in computing the extension. These features lead to a particularly simple implementation using either a sparse direct solver or a matrix-free conjugate gradient solver. Furthermore, we propose a fast Poisson preconditioner that can be used to accelerate the convergence of the latter.We demonstrate the biharmonic extension on a number of test problems that serve to illustrate its effectiveness at producing smooth and accurate extensions near interfaces. A further feature of the method is the natural way in which it deals with symmetry and periodicity, ensuring through its construction that the extension field also respects these symmetries.

Direct‐Iterative Hybrid Solution in Nonlinear Dynamic Structural Analysis

AbstractAlthough many advanced sparse direct solvers are widely used in structural analysis, these often require longer computing times than iterative solvers for well‐conditioned structural systems. However, iterative solvers cannot efficiently solve an ill‐conditioned system when a structure becomes highly nonlinear. This work proposes a hybrid solution integrating a direct solver and an iterative solver to reduce overall computing time in solving a series of linear equations arising from nonlinear dynamic structural analysis. The hybrid solution selects the iterative solver, which reuses the factorized matrices for preconditioning, and switches to the direct solver in the initial stage or when factorized matrices need updating. The performance of the hybrid solution is tested on the OpenSees platform. The results show that the hybrid solution outperforms the direct solver, even when the structure becomes highly nonlinear during analysis.

Sparse direct solver for large finite element problems based on the minimum degree algorithm

A sparse direct solver for large problems from solid continuum mechanics based on the minimum degree algorithm is proposed and tested. The solver is designed to take advantage of the properties of the finite element method, particularly the structure of the finite element mesh. For the minimization of the fill-in in the matrix factors a modification of the approximate minimum degree ordering algorithm of Amestoy, Davis and Duff is utilized. The employed sparse matrix storage format and the algorithms for each of the solver phases are also described. The results of numerical tests of the solver on large real-world finite element problems are presented and its performance is compared to a frontal solver and the PARDISO sparse direct solver.