Fast Direct Solver Research Articles

Second kind integral equation formulation for the modified biharmonic equation and its applications

A system of Fredholm second kind integral equations (SKIEs) is constructed for the modified biharmonic equation in two dimensions with gradient boundary conditions. Such boundary value problem arises naturally when the incompressible Navier–Stokes equations are solved via a semi-implicit discretization scheme and the resulting boundary value problem at each time step is then solved using the pure stream-function formulation. The advantages of such an approach (Greengard and Kropinski, 1998) [14] are two fold: first, the velocity is automatically divergence free, and second, complicated (nonlocal) boundary conditions for the vorticity are avoided. Our construction, though similar to that of Farkas (1989) [10] for the biharmonic equation, is modified such that the SKIE formulation has low condition numbers for large values of the parameter. We illustrate the performance of our numerical scheme with several numerical examples. Finally, the scheme can be easily coupled with standard fast algorithms such as FFT, fast multipole methods (Greengard and Rokhlin, 1987) [15], or fast direct solvers (Ho and Greengard, 2012; Martinsson and Rokhlin, 2005) [17,25] to achieve optimal complexity, bringing accurate large-scale long-time fluid simulations within practical reach.

A fast direct solver for quasi-periodic scattering problems

We consider the numerical solution of the scattering of time-harmonic plane waves from an infinite periodic array of reflection or transmission obstacles in a homogeneous background medium, in two dimensions. Boundary integral formulations are ideal since they reduce the problem to N unknowns on the obstacle boundary. However, for complex geometries and/or higher frequencies the resulting dense linear system becomes large, ruling out dense direct methods, and often ill-conditioned (despite being 2nd-kind), rendering fast multipole-based iterative schemes also inefficient. We present an integral equation based solver with O(N) complexity, which handles such ill-conditioning, using recent advances in “fast” direct linear algebra to invert hierarchically the isolated obstacle matrix. This is combined with a recent periodizing scheme that is robust for all incident angles, including Wood’s anomalies, based upon the free space Green’s function kernel. The resulting solver is extremely efficient when multiple incident angles are needed, as occurs in many applications. Our numerical tests include a complicated obstacle several wavelengths in size, with N=105 and solution error of 10-10, where the solver is 66 times faster per incident angle than a fast multipole based iterative solution, and 600 times faster when incident angles are chosen to share Bloch phases.

An $$\mathcal O (N \log N)$$ O ( N log N ) Fast Direct Solver for Partial Hierarchically Semi-Separable Matrices

This article describes a fast direct solver (i.e., not iterative) for partial hierarchically semi-separable systems. This solver requires a storage of $$\mathcal O (N \log N)$$ O ( N log N ) and has a computational complexity of $$\mathcal O (N \log N)$$ O ( N log N ) arithmetic operations. The numerical benchmarks presented illustrate the method in the context of interpolation using radial basis functions. The key ingredients behind this fast solver are recursion, efficient low rank factorization using Chebyshev interpolation, and the Sherman---Morrison---Woodbury formula. The algorithm and the analysis are worked out in detail. The performance of the algorithm is illustrated for a variety of radial basis functions and target accuracies.

Fast ${\cal H}$-Matrix-Based Direct Integral Equation Solver With Reduced Computational Cost for Large-Scale Interconnect Extraction

In this paper, we propose a fast H-matrix-based direct solution with a significantly reduced computational cost for an integral-equation-based capacitance extraction of large-scale 3-D interconnects in multiple dielectrics. We reduce the computational cost of an H-matrix-based computation by simultaneously optimizing the H-matrix partition to minimize the number of matrix blocks and minimizing the rank of each matrix block based on a prescribed accuracy. With the proposed cost-reduction method, we develop a fast LU-based direct solver. This solver possesses a complexity of kCspO (NlogN) in storage, a complexity of k2Csp2O(Nlog2N) in LU factorization, and a complexity of kCspO(NlogN) in LU solution, where k is the maximal rank, Csp is a constant dependent on matrix partition, and the constant kCsp is minimized based on accuracy by the proposed cost-reduction method. The proposed solver successfully factorizes dense matrices that involve millions of unknowns in fast CPU time and modest memory consumption, and with the prescribed accuracy satisfied. As an algebraic method, the underlying fast technique is kernel independent.

2008 2次元Laplace方程式のDirichlet-Neumann混合境界値問題におけるQR分解を用いた境界要素法の高速直接解法について

2008 2次元Laplace方程式のDirichlet-Neumann混合境界値問題におけるQR分解を用いた境界要素法の高速直接解法について

Parallelizing the Conjugate Gradient Algorithm for Multilevel Toeplitz Systems

Multilevel Toeplitz linear systems appear in a wide range of scientific and engineering applications. While several fast direct solvers exist for the basic 1-level Toeplitz matrices, in the multilevel case an iterative solver provides the most general and practical solution. This paper proposes several parallelization techniques that enable an efficient implementation of the conjugate gradient algorithm for solving multilevel Toeplitz systems on distributed-memory machines. The two major differences between this implementation and that for a general sparse linear solver are (1) a communication-efficient approach to handle data expansion and truncation and data transpose simultaneously; (2) the interleaving of matrix-vector multiplications and vector inner product calculations to reduce synchronization cost and latency. Similar ideas can be applied to the implementation of other iterative methods for Toeplitz systems that are not necessarily symmetric positive definite. Scaling results are shown to demonstrate the usefulness of the proposed techniques.

A Fast Direct Solver for a Class of 3-D Elliptic Partial Differential Equation with Variable Coefficient

AbstractWe propose a direct solver for the three-dimensional Poisson equation with a variable coefficient, and an algorithm to directly solve the associated sparse linear systems that exploits the sparsity pattern of the coefficient matrix. Introducing some appropriate finite difference operators, we derive a second-order scheme for the solver, and then two suitable high-order compact schemes are also discussed. For a cube containing N nodes, the solver requires arithmetic operations and memory to store the necessary information. Its efficiency is illustrated with examples, and the numerical results are analysed.

A fast direct matrix solver for surface integral equation methods for electromagnetic wave scattering from non‐penetrable targets

The implementation details of a fast direct solver is described herein for solving dense matrix equations from the application of surface integral equation methods for electromagnetic field scatterings from non‐penetrable targets. The proposed algorithm exploits the smoothness of the far field and computes a low rank decomposition of the off‐diagonal coupling blocks of the matrices through a set of skeletonization processes. Moreover, an artificial surface (the Huygens' surface) is introduced for each clustering group to efficiently account for the couplings between well‐separated groups. Furthermore, a recursive multilevel version of the algorithm is presented. Although asymptotically the algorithm would not alter the bleak outlook of the complexity of the worst case scenario,O(N3) for required CPU time where N denotes the number of unknowns, for electrically large electromagnetic (EM) problems; through numerical examples, we found that the proposed multilevel direct solver can scale as good as O(N1.3) in memory consumption and O(N1.8) in CPU time for moderate‐sized EM problems. Note that our conclusions are drawn based on a few sample examples that we have conducted and should not be taken as a true complexity analysis for general electrodynamic applications. However, for the fixed frequency (h‐refinement) scenario, where the discretization size decreases, the computational complexities observed agree well with the theoretical predictions. Namely, the algorithm exhibits O(N) and O(N1.5) complexities for memory consumption and CPU time, respectively.

Sparse Matrix Factorization Using Overlapped Localizing LOGOS Modes on a Shifted Grid

Different types of local global solution (LOGOS) modes have been proposed to develop fast direct solvers for electromagnetic wave problems. While radiating LOGOS modes are suitable for high frequency problems, this paper focuses on the localizing LOGOS modes, which are suitable for low-to-mid frequency problems. Factorizations relying on non-overlapped, localizing LOGOS (NL-LOGOS) modes are provably error controllable, but fail to provide acceptable asymptotic computational efficiencies. Factorizations using overlapped, localizing LOGOS (OL-LOGOS) modes provide improved asymptotic computational performance, but have been observed to suffer from significant overhead and, in some cases, poorer error control. This paper describes a new procedure for obtaining the OL-LOGOS modes. It is observed that the asymptotic efficiency is similar to previously reported OL-LOGOS schemes while providing stable error control for standard problem formulations. The overhead associated with computing the OL-LOGOS modes is also reduced significantly.

Hierarchical Matrices Method and Its Application in Electromagnetic Integral Equations

Hierarchical (H-) matrices method is a general mathematical framework providing a highly compact representation and efficient numerical arithmetic. When applied in integral-equation- (IE-) based computational electromagnetics,H-matrices can be regarded as a fast algorithm; therefore, both the CPU time and memory requirement are reduced significantly. Its kernel independent feature also makes it suitable for any kind of integral equation. To solveH-matrices system, Krylov iteration methods can be employed with appropriate preconditioners, and direct solvers based on the hierarchical structure ofH-matrices are also available along with high efficiency and accuracy, which is a unique advantage compared to other fast algorithms. In this paper, a novel sparse approximate inverse (SAI) preconditioner in multilevel fashion is proposed to accelerate the convergence rate of Krylov iterations for solvingH-matrices system in electromagnetic applications, and a group of parallel fast direct solvers are developed for dealing with multiple right-hand-side cases. Finally, numerical experiments are given to demonstrate the advantages of the proposed multilevel preconditioner compared to conventional “single level” preconditioners and the practicability of the fast direct solvers for arbitrary complex structures.

Generalized Equivalence Integral Equations

A generalized equivalence integral equation (GEIE) approach to formulating scattering from essentially convex closed surfaces is proposed. The GEIE approach invokes the generalized surface field equivalence to partially fill the volume originally occupied by the scatterer with judiciously selected materials, as opposed to the conventional replacement of the scatterer by the free space. The type and shape of the material inclusions can be selected to allow for a numerically efficient construction of the modified Green's function. Introduction of impenetrable and lossy materials confines the field interaction along the scatterer surface and reduces the coupling between the distant parts of the scatterer, which essentially makes the impedance matrix banded. The presence of lossy materials also resolves the nonuniqueness problem of the electric and magnetic field integral equations by eliminating the internal resonances. The formulation provides a pathway for developing fast iterative and direct electromagnetic integral equation solvers.

A Fast Direct Solver for Structured Linear Systems by Recursive Skeletonization

We present a fast direct solver for structured linear systems based on multilevel matrix compression. Using the recently developed interpolative decomposition of a low-rank matrix in a recursive manner, we embed an approximation of the original matrix into a larger, but highly structured sparse one that allows fast factorization and application of the inverse. The algorithm extends the Martinsson/Rokhlin method developed for 2D boundary integral equations and proceeds in two phases: a precomputation phase, consisting of matrix compression and factorization, followed by a solution phase to apply the matrix inverse. For boundary integral equations which are not too oscillatory, e.g., based on the Green's functions for the Laplace or low-frequency Helmholtz equations, both phases typically have complexity O(N) in two dimensions, where $N$ is the number of discretization points. In our current implementation, the corresponding costs in three dimensions are $O(N^{3/2})$ and $O(N \log N)$ for precomputation and solution, respectively. Extensive numerical experiments show a speedup of $\sim 100$ for the solution phase over modern fast multipole methods; however, the cost of precomputation remains high. Thus, the solver is particularly suited to problems where large numbers of iterations would be required. Such is the case with ill-conditioned linear systems or when the same system is to be solved with multiple right-hand sides. Our algorithm is implemented in Fortran and freely available.

A fast direct solver for the integral equations of scattering theory on planar curves with corners

We describe an approach to the numerical solution of the integral equations of scattering theory on planar curves with corners. It is rather comprehensive in that it applies to a wide variety of boundary value problems; here, we treat the Neumann and Dirichlet problems as well as the boundary value problem arising from acoustic scattering at the interface of two fluids. It achieves high accuracy, is applicable to large-scale problems and, perhaps most importantly, does not require asymptotic estimates for solutions. Instead, the singularities of solutions are resolved numerically. The approach is efficient, however, only in the low- and mid-frequency regimes. Once the scatterer becomes more than several hundred wavelengths in size, the performance of the algorithm of this paper deteriorates significantly. We illustrate our method with several numerical experiments, including the solution of a Neumann problem for the Helmholtz equation given on a domain with nearly 10000 corner points.

A fast direct solver for elliptic problems on general meshes in 2D

We present a fast direct algorithm for solutions to linear systems arising from 2D elliptic equations. We follow the approach in Xia et al. (2009) on combining the multifrontal method with hierarchical matrices. We present a variant of that approach with additional hierarchical structure, extend it to quasi-uniform meshes, and detail an adaptive decomposition procedure for general meshes. Linear time complexity is shown for a quasi-regular grid and demonstrated via numerical results for the adaptive algorithm.

An adaptive fast direct solver for boundary integral equations in two dimensions

We describe an algorithm for the rapid direct solution of linear algebraic systems arising from the discretization of boundary integral equations of potential theory in two dimensions. The algorithm is combined with a scheme that adaptively rearranges the parameterization of the boundary in order to minimize the ranks of the off-diagonal blocks in the discretized operator, thus obviating the need for the user to supply a parameterization r of the boundary for which the distance ‖ r ( s ) − r ( t ) ‖ between two points on the boundary is related to their corresponding distance | s − t | in the parameter space. The algorithm has an asymptotic complexity of O ( N log 2 N ) , where N is the number of nodes in the discretization. The performance of the algorithm is illustrated with several numerical examples.

Tuning a Hybrid GPU-CPU V-Cycle Multilevel Preconditioner for Solving Large Real and Complex Systems of FEM Equations

This letter presents techniques for tuning an accelerated preconditioned conjugate gradient solver with a multilevel preconditioner. The solver is optimized for a fast solution of sparse systems of equations arising in computational electromagnetics in a finite element method using higher-order elements. The goal of the tuning is to increase the throughput while at the same time reducing the memory requirements in order to allow one to process very large complex or real systems in single and double precision using commodity graphic processing units (GPUs). A threefold memory footprint reduction is achieved by means of a new format of storing sparse matrices. The acceleration is achieved by optimizing a sparse matrix-vector product on a GPU by applying new features of the Fermi architecture. Further improvements are obtained by introducing more levels into the preconditioner and the application of a fast sparse direct solver for the operations executed on a CPU. Numerical results for a setup consisting of a Fermi GPU (GTX 480) and a Xeon six-core CPU showed that the proposed approach allows one to handle systems involving millions of unknowns and reach the speedup factor of almost 4 compared to the CPU-only implementation.

Existence of ${\cal H}$-Matrix Representations of the Inverse Finite-Element Matrix of Electrodynamic Problems and ${\cal H}$-Based Fast Direct Finite-Element Solvers

In this work, we prove that the sparse matrix resulting from a finite-element-based analysis of electrodynamic problems can be represented by an <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">H</i> matrix without any approximation, and the inverse of this sparse matrix has a data-sparse <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">H</i> -matrix approximation with error well controlled. Two proofs are developed. One is based on the general eigenvalue-based solution to the ordinary differential equations, and the other is based on the relationship between a partial differential operator and an integral operator. Both proofs have reached the same conclusion. Based on the proof, we develop an <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">H</i> -matrix-based direct finite-element solver of <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">O</i> ( <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">kN</i> log <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> ) memory complexity and <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">O</i> ( <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</i> <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sup> <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> log <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sup> <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> ) time complexity for solving electromagnetic problems, where <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">k</i> is a small parameter that is adaptively determined based on accuracy requirements, and <i xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> is the number of unknowns. Both inverse-based and LU-based direct solutions are developed. The LU-based solution is further accelerated by nested dissection. A comparison with the state-of-the-art direct finite element solver that employs the most advanced sparse matrix solution has shown clear advantages of the proposed direct solver. In addition, the proposed solver is applicable to arbitrarily-shaped three-dimensional structures and arbitrary inhomogeneity.

Numerical methods for on-line power system load flow analysis

Newton-Raphson method is the most widely accepted load flow solution algorithm. However LU factorization remains a computationally challenging task to meet the real-time needs of the power system. This paper proposes the application of very fast multifrontal direct linear solvers for solving the linear system sub-problem of power system real-time load flow analysis by utilizing the state-of-the-art algorithms for ordering and preprocessing. Additionally the unsymmetric multifrontal method for LU factorization and highly optimized Intel® Math Kernel Library BLAS has been used. Two state-of-the-art multifrontal algorithms for unsymmetric matrices namely UMFPACK V5.2.0 and sequential MUMPS 4.8.3 (“Multifrontal Massively Parallel Solver”) are customized for the AC power system Newton-Raphson based load flow analysis. The multifrontal solvers are compared against the state-of-the-art sparse Gaussian Elimination based HSL sparse solver MA48. This study evaluates the performance of above multifrontal solvers in terms of number of factors, computational time, number of floating-point operations and memory, in the context of load flow solution on nine systems including very large real power systems. The results of the performance evaluation are reported. The proposed method achieves significant reduction in computational time.

Approximate Inverse Circulant-plus-Diagonal Preconditioners for Toeplitz-plus-Diagonal Matrices

We consider the solutions of Hermitian positive definite Toeplitz-plus-diagonal systems $(T+D)x=b$, where $T$ is a Toeplitz matrix and $D$ is diagonal and positive. However, unlike the case of Toeplitz systems, no fast direct solvers have been developed for solving them. In this paper, we employ the preconditioned conjugate gradient method with approximate inverse circulant-plus-diagonal preconditioners to solving such systems. The proposed preconditioner can be constructed and implemented efficiently using fast Fourier transforms. We show that if the entries of $T$ decay away exponentially from the main diagonals, the preconditioned conjugate gradient method applied to the preconditioned system converges very quickly. Numerical examples including spatial regularization for image deconvolution application are given to illustrate the effectiveness of the proposed preconditioner.

An H-LU Based Direct Finite Element Solver Accelerated by Nested Dissection for Large-scale Modeling of ICs and Packages

In this work, we prove that for the sparse matrix resulting from a flnite-element- based analysis of electrodynamic problems, its inverse has a data-sparse H-matrix approximation with error well controlled. Based on this proof, we develop a fast direct flnite element solver. In this direct solver, the H-matrix-based LU factorization is developed, which is further accelerated by nested dissection. We show that the proposed direct solver has an O(kNlogN) memory complexity and O(k 2 Nlog 2 N) time complexity, where k is a small number that is adaptively determined based on accuracy requirements, and N is the number of unknowns. A comparison with the state-of-the-art direct flnite element solver that employs the most advanced sparse matrix solution has shown a clear advantage of the proposed solver. Applications to large-scale package modeling involving millions of unknowns have demonstrated the accuracy and almost linear complexity of the proposed direct solver. In addition, the proposed method is applicable to arbitrarily-shaped three-dimensional structures and arbitrary inhomogeneity.