Iterative Krylov Research Articles

AFMPB: An adaptive fast multipole Poisson–Boltzmann solver for calculating electrostatics in biomolecular systems

A Fortran program package is introduced for rapid evaluation of the electrostatic potentials and forces in biomolecular systems modeled by the linearized Poisson–Boltzmann equation. The numerical solver utilizes a well-conditioned boundary integral equation (BIE) formulation, a node-patch discretization scheme, a Krylov subspace iterative solver package with reverse communication protocols, and an adaptive new version of the fast multipole method in which the exponential expansions are used to diagonalize the multipole-to-local translations. The program and its full description, as well as several closely related libraries and utility tools are available at http://lsec.cc.ac.cn/~lubz/afmpb.html and a mirror site at http://mccammon.ucsd.edu/. This paper is a brief summary of the program: the algorithms, the implementation and the usage. New version program summaryProgram title: AFMPBCatalogue identifier: AEGB_v1_1Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEGB_v1_1.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU General Public License, version 2No. of lines in distributed program, including test data, etc.: 440784No. of bytes in distributed program, including test data, etc.: 8187139Distribution format: tar.gzProgramming language: FortranComputer: AnyOperating system: AnyRAM: Depends on the size of the discretized biomolecular systemClassification: 3External routines: Pre- and post-processing tools are required for generating the boundary elements and for visualization. Users can use MSMS (http://www.scripps.edu/sanner/html/msmshome.html) for pre-processing, and VMD (http://www.ks.uiuc.edu/Research/vmd/) for visualization. Sub-programs included: An iterative Krylov subspace solvers package from SPARSKIT by Yousef Saad (http://www-users.cs.umn.edu/saad/software/SPARSKIT/sparskit.html), and the fast multipole methods subroutines from FMMSuite (http://www.fastmultipole.org/).Catalogue identifier of previous version: AEGB_v1_0Journal reference of previous version: Comput. Phys. Comm. 181 (2010) 1150Does the new version supersede the previous version?: YesNature of problem: Numerical solution of the linearized Poisson–Boltzmann equation that describes electrostatic interactions of molecular systems in ionic solutions.Solution method: A novel node-patch scheme is used to discretize the well-conditioned boundary integral equation formulation of the linearized Poisson–Boltzmann equation. Various Krylov subspace solvers can be subsequently applied to solve the resulting linear system, with a bounded number of iterations independent of the number of discretized unknowns. The matrix–vector multiplication at each iteration is accelerated by the adaptive new versions of fast multipole methods. The AFMPB solver requires other stand-alone pre-processing tools for boundary mesh generation, post-processing tools for data analysis and visualization, and can be conveniently coupled with different time stepping methods for dynamics simulation.Reasons for new version: Some bugs are fixed in the new version.Summary of revisions:•The type definition of ippt1 in line 88 of FBEM/bempb.f and line 32 of FBEM/closecoef.f is changed from real *8 to integer*4, and a similar change is made for ippt in line 105 of FBEM/solvpb.f and in line 32 of FBEM/closecoef.f.•In FBEM/elmgeom.f, line 239 “ELSEIF (meshfmt.EQ.1.OR. meshfmt.EQ. 4.OR. meshfmt.EQ.5) THEN” is changed to “ELSEIF (meshfmt.EQ.1.OR. meshfmt.EQ. 4) THEN”, line 478 “KJ=IDFCL(K)+J-1” is changed to “i=IDFCL(K)+J-1”, line 479 “KJ=NE(KJ)” is changed to “KJ=NE(i)”, line 480 “KJ1=NE(KJ+1)” is changed to “KJ1=NE(i+1)”, and line 647 “ STOP” is changed to “c STOP”.•Five subroutines in FMM part (syukadap.f, syukdn.f, slapadap.f, slapdn.f, and treeadap.f) are substituted with the new ones in the new version.Restrictions: Only three or six significant digits options are provided in this version.Unusual features: Most of the codes are in Fortran77 style. Memory allocation functions from Fortran90 and above are used in a few subroutines.Additional comments: The current version of the codes is designed and written for single core/processor desktop machines. Check http://lsec.cc.ac.cn/lubz/afmpb.html for updates and changes.Running time: The running time varies with the number of discretized elements (N) in the system and their distributions. In most cases, it scales linearly as a function of N.

Project APhiD: A Lorenz-gauged A-[formula omitted] decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth

An essential element for computational hypothesis testing, data inversion and experiment design for electromagnetic geophysics is a robust forward solver, capable of easily and quickly evaluating the electromagnetic response of arbitrary geologic structure. The usefulness of such a solver hinges on the balance among competing desires like ease of use, speed of forward calculation, scalability to large problems or compute clusters, parsimonious use of memory access, accuracy and by necessity, the ability to faithfully accommodate a broad range of geologic scenarios over extremes in length scale and frequency content. This is indeed a tall order. The present study addresses recent progress toward the development of a forward solver with these properties. Based on the Lorenz-gauged Helmholtz decomposition, a new finite volume solution over Cartesian model domains endowed with complex-valued electrical properties is shown to be stable over the frequency range 10−2–1010Hz and range 10−3–105m in length scale. Benchmark examples are drawn from magnetotellurics, exploration geophysics, geotechnical mapping and laboratory-scale analysis, showing excellent agreement with reference analytic solutions. Computational efficiency is achieved through use of a matrix-free implementation of the quasi-minimum-residual (QMR) iterative solver, which eliminates explicit storage of finite volume matrix elements in favor of “on the fly” computation as needed by the iterative Krylov sequence. Further efficiency is achieved through sparse coupling matrices between the vector and scalar potentials whose non-zero elements arise only in those parts of the model domain where the conductivity gradient is non-zero. Multi-thread parallelization in the QMR solver through OpenMP pragmas is used to reduce the computational cost of its most expensive step: the single matrix–vector product at each iteration. High-level MPI communicators farm independent processes to available compute nodes for simultaneous computation of multi-frequency or multi-transmitter responses.

Local Fourier analysis of the complex shifted Laplacian preconditioner for Helmholtz problems

SUMMARYIn this paper, we solve the Helmholtz equation with multigrid preconditioned Krylov subspace methods. The class of shifted Laplacian preconditioners is known to significantly speed up Krylov convergence. However, these preconditioners have a parameter , a measure of the complex shift. Because of contradictory requirements for the multigrid and Krylov convergence, the choice of this shift parameter can be a bottleneck in applying the method. In this paper, we propose a wavenumber‐dependent minimal complex shift parameter, which is predicted by a rigorous k‐grid local Fourier analysis (LFA) of the multigrid scheme. We claim that, given any (regionally constant) wavenumber, this minimal complex shift parameter provides the reader with a parameter choice that leads to efficient Krylov convergence. Numerical experiments in one and two spatial dimensions validate the theoretical results. It appears that the proposed complex shift is both the minimal requirement for a multigrid V‐cycle to converge and being near optimal in terms of Krylov iteration count. Copyright © 2013 John Wiley & Sons, Ltd.

Communication Avoiding ILU0 Preconditioner

In this paper we present a communication avoiding ILU0 preconditioner for solving large linear systems of equations by using iterative Krylov subspace methods. Recent research has focused on communication avoiding Krylov subspace methods based on so-called s-step methods. However, there are not many communication avoiding preconditioners yet, and this represents a serious limitation of these methods. Our preconditioner allows us to perform $s$ iterations of the iterative method with no communication, through ghosting some of the input data and performing redundant computation. To avoid communication, an alternating reordering algorithm is introduced for structured and well partitioned unstructured matrices, which requires the input matrix to be ordered by using a graph partitioning technique such as k-way or nested dissection. We show that the reordering does not affect the convergence rate of the ILU0 preconditioned system as compared to k-way or nested dissection ordering, while it reduces data movement and is expected to reduce the time needed to solve a linear system. In addition to communication avoiding Krylov subspace methods, our preconditioner can be used with classical methods such as GMRES to reduce communication.

Multigrid in energy preconditioner for Krylov solvers

We have added a new multigrid in energy (MGE) preconditioner to the Denovo discrete-ordinates radiation transport code. This preconditioner takes advantage of a new multilevel parallel decomposition. A multigroup Krylov subspace iterative solver that is decomposed in energy as well as space-angle forms the backbone of the transport solves in Denovo. The space-angle-energy decomposition facilitates scaling to hundreds of thousands of cores. The multigrid in energy preconditioner scales well in the energy dimension and significantly reduces the number of Krylov iterations required for convergence. This preconditioner is well-suited for use with advanced eigenvalue solvers such as Rayleigh Quotient Iteration and Arnoldi.

Hierarchical Schur complement preconditioner for the stochastic Galerkin finite element methods

Use of the stochastic Galerkin finite element methods leads to large systems of linear equations obtained by the discretization of tensor product solution spaces along their spatial and stochastic dimensions. These systems are typically solved iteratively by a Krylov subspace method. We propose a preconditioner which takes an advantage of the recursive hierarchy in the structure of the global matrices. In particular, the matrices posses a recursive hierarchical two-by-two structure, with one of the submatrices block diagonal. Each one of the diagonal blocks in this submatrix is closely related to the deterministic mean-value problem, and the action of its inverse is in the implementation approximated by inner loops of Krylov iterations. Thus our hierarchical Schur complement preconditioner combines, on each level in the approximation of the hierarchical structure of the global matrix, the idea of Schur complement with loops for a number of mutually independent inner Krylov iterations, and several matrix-vector multiplications for the off-diagonal blocks. Neither the global matrix, nor the matrix of the preconditioner need to be formed explicitly. The ingredients include only the number of stiffness matrices from the truncated Karhunen-Lo\`{e}ve expansion and a good preconditioned for the mean-value deterministic problem. We provide a condition number bound for a model elliptic problem and the performance of the method is illustrated by numerical experiments.

A Parallel Adaptive Newton-Krylov-Schwarz Method for 3D Compressible Inviscid Flow Simulations

A parallel adaptive pseudo transient Newton-Krylov-Schwarz (αΨNKS) method for the solution of compressible flows is presented. Multidimensional upwind residual distribution schemes are used for space discretisation, while an implicit time-marching scheme is employed for the discretisation of the (pseudo)time derivative. The linear system arising from the Newton method applied to the resulting nonlinear system is solved by the means of Krylov iterations with Schwarz-type preconditioners. A scalable and efficient data structure for theαΨNKS procedure is presented. The main computational kernels are considered, and an extensive analysis is reported to compare the Krylov accelerators, the preconditioning techniques. Results, obtained on a distributed memory computer, are presented for 2D and 3D problems of aeronautical interest on unstructured grids.

Hybrid Deterministic/Monte Carlo Neutronics

In this paper we describe a hybrid deterministic/Monte Carlo algorithm for neutron transport simulation. The algorithm is based on nonlinear accelerators for source iteration, using Monte Carlo methods for the purely absorbing high-order problem and a Jacobian-free Newton--Krylov iteration for the low-order problem. We couple the Monte Carlo solution with the low-order problem using filtering to smooth the flux and current from the Monte Carlo solver and an analytic Jacobian-vector product to avoid numerical differentiation of the Monte Carlo results. We use a continuous energy deposition tally for the Monte Carlo simulation. We conclude the paper with numerical results which illustrate the effectiveness of the new algorithm.

Hierarchical Preconditioning for PTAP‐Systems arising from particulate flows

AbstractParticulate flows, i.e. flow of an incompressible, Newtonian carrier fluid loaded with (many) rigid bodies, play an important role in diverse technical applications. In [1] a finite element method was introduced to simulate such flows. The method relies on a splitting method and a subspace projection method to incorporate the rigid body motion of the particles in a so called one domain approach [2]. The resulting systems arising after discretization are ill conditioned in general. Thus preconditioning is mandatory for instance when using an iterative Krylov subspace methods. In this paper we present a Schur complement based preconditioner well suited for this type of application. (© 2012 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

An optimal linear solver for the Jacobian system of the extreme type-II Ginzburg–Landau problem

This paper considers the extreme type-II Ginzburg–Landau equations, a nonlinear PDE model for describing the states of a wide range of superconductors. Based on properties of the Jacobian operator and an AMG strategy, a preconditioned Newton–Krylov method is constructed. After a finite-volume-type discretization, numerical experiments are done for representative two- and three-dimensional domains. Strong numerical evidence is provided that the number of Krylov iterations is independent of the dimension n of the solution space, yielding an overall solver complexity of O(n).

Comparative performance of exponential, implicit, and explicit integrators for stiff systems of ODEs

Exponential integrators have enjoyed a resurgence of interest in recent years, but there is still limited understanding of how their performance compares with that of state-of-the-art integrators, most notably the commonly used Newton–Krylov implicit methods. In this paper we present comparative performance analysis of Krylov-based exponential, implicit and explicit integrators on a suite of stiff test problems and demonstrate that exponential integrators have computational advantages compared to the other methods, particularly as problems become larger and more stiff. We argue that the faster convergence of the Krylov iteration within exponential integrators accounts for the main proportion of the computational savings that they provide and illustrate how the structure of these methods ensures such efficiency. In addition, we demonstrate the computational advantages of the newly introduced Tokman and Loffeld (2010) [17] exponential propagation Runge–Kutta (EpiRK) fifth-order methods. The detailed analysis of the performance of the methods that is presented provides guidelines for the construction and implementation of efficient exponential methods and the quantitative comparisons inform the selection of appropriate schemes for other problems.

Comparison between iteration schemes for three-dimensional coordinate-transformed saturated–unsaturated flow model

Three different iteration methods for a three-dimensional coordinate-transformed saturated–unsaturated flow model are compared in this study. The Picard and Newton iteration methods are the common approaches for solving Richards’ equation. The Picard method is simple to implement and cost-efficient (on an individual iteration basis). However it converges slower than the Newton method. On the other hand, although the Newton method converges faster, it is more complex to implement and consumes more CPU resources per iteration than the Picard method. The comparison of the two methods in finite-element model (FEM) for saturated–unsaturated flow has been well evaluated in previous studies. However, two iteration methods might exhibit different behavior in the coordinate-transformed finite-difference model (FDM). In addition, the Newton–Krylov method could be a suitable alternative for the coordinate-transformed FDM because it requires the evaluation of a 19-point stencil matrix. The formation of a 19-point stencil is quite a complex and laborious procedure. Instead, the Newton–Krylov method calculates the matrix–vector product, which can be easily approximated by calculating the differences of the original nonlinear function. In this respect, the Newton–Krylov method might be the most appropriate iteration method for coordinate-transformed FDM. However, this method involves the additional cost of taking an approximation at each Krylov iteration in the Newton–Krylov method. In this paper, we evaluated the efficiency and robustness of three iteration methods—the Picard, Newton, and Newton–Krylov methods—for simulating saturated–unsaturated flow through porous media using a three-dimensional coordinate-transformed FDM.

Алгоритмы решения СЛАУ на системах с распределенной памятью в применении к задачам электромагнетизма

Paper presents various aspects of harmonic electromagnetic fields simulation on clusters. The major computational complexity comes from the solution of the systems of linear algebraic equations (SLAEs) arising from the approximations of corresponding electromagnetic boundary value problems by Nedelec elements of various orders. Effective and efficient approaches to the decomposition of the computational domain and the matrix of the system are considered. Distributed SLAEs are solved using iterative Krylov subspace methods preconditioned by additive Schwarz method. In order to increase the effectiveness of the algorithms iterations are performed in the trace space. Implementation of the solvers is based on MPI for data transfers. The solution of the systems in subdomains is performed by PARDISO direct solver from IntelR MKL library. Numerical experiments results on a series of model and real-life problems show the effectiveness of the presented algorithms.

Old and new parameter choice rules for discrete ill-posed problems

Linear discrete ill-posed problems are difficult to solve numerically because their solution is very sensitive to perturbations, which may stem from errors in the data and from round-off errors introduced during the solution process. The computation of a meaningful approximate solution requires that the given problem be replaced by a nearby problem that is less sensitive to disturbances. This replacement is known as regularization. A regularization parameter determines how much the regularized problem differs from the original one. The proper choice of this parameter is important for the quality of the computed solution. This paper studies the performance of known and new approaches to choosing a suitable value of the regularization parameter for the truncated singular value decomposition method and for the LSQR iterative Krylov subspace method in the situation when no accurate estimate of the norm of the error in the data is available. The regularization parameter choice rules considered include several L-curve methods, Regińska’s method and a modification thereof, extrapolation methods, the quasi-optimality criterion, rules designed for use with LSQR, as well as hybrid methods.

3D finite-difference modeling of elastic wave propagation in the Laplace-Fourier domain

With the recent interest in the Laplace-Fourier domain full waveform inversion, we have developed new heterogeneous 3D fourth- and second-order staggered-grid finite-difference schemes for modeling seismic wave propagation in the Laplace-Fourier domain. Our approach is based on the integro-interpolation technique for the velocity-stress formulation in the Cartesian coordinate system. Five averaging elastic coefficients and three averaging densities are necessary to describe the heterogeneous medium, with harmonic averaging of the bulk and shear moduli, and arithmetic averaging of density. In the fourth-order approximation, we improved the accuracy of the scheme using a combination of integral identities for two elementary volumes — “small” and “large” around spatial gridpoints where the wave variables are defined. Two solution approaches are provided, both of which are solved with iterative Krylov methods. In the first approach, the stress variables are eliminated and a linear system for the velocity components is solved. In the second approach, we worked directly with the first-order system of velocity and stress variables. This reduced the computer memory required to store the complex matrix, along with reducing (by 30%) the number of arithmetic operations needed for the solution compared to the fourth-order scheme for velocity only. Numerical examples show that our finite-difference formulations for elastic wavefield simulations can achieve more accurate solutions with fewer grid points than those from previously published second and fourth-order frequency-domain schemes. We applied our simulator to the investigation of wavefields from the SEG/EAGE model in the Laplace-Fourier domain. The calculation is sensitive to the heterogeneity of the medium and capable of describing the structures of complex objects. Our technique can also be extended to 3D elastic modeling within the time domain.

Convergence of the Iterative Rational Krylov Algorithm

The iterative rational Krylov algorithm (IRKA) of Gugercin et al. (2008) [8] is an interpolatory model reduction approach to the optimal H2 approximation problem. Even though the method has been illustrated to show rapid convergence in various examples, a proof of convergence has not been provided yet. In this note, we show that in the case of state-space-symmetric systems, IRKA is a locally convergent fixed-point iteration to a local minimum of the underlying H2 approximation problem.

A Refined Arnoldi Algorithm Based Krylov Subspace Technique for MEMS Model Order Reduction

A refined approach producing MEMS numerical macromodels is proposed in this paper by generating the iterative Krylov subspace using a refined Arnoldi algorithm, which can reduce the degrees of freedom of the original system equations described by the state space method. Projection of the original system matrix onto the Krylov subspace which is spanned by a refined Arnoldi algorithm is still based on the transfer function moment matching principle. The idea of the iterative version is to expect that a new initial vector will contain more and more information on the required eigenvectors that is called refined vector. The refined approach improves approximation accuracy of the system matrix eigenvalues equivalent to a more accurate approximation to the poles of the system transfer function, obtaining a more accurate reduced-order model. The clamped beam model and the FOM model are reduced order by classical Arnoldi and refined Arnoldi algorithm in numerical experiments. From the computing result it is concluded that the refined Arnoldi algorithm based Krylov subspace technique for MEMS model order reduction has more accuracy and reaches lower order number of reduced order model than the classical Arnoldi process.

Optimizing large parameter sets in variational quantum Monte Carlo

We present a technique for optimizing hundreds of thousands of variational parameters in variational quantum Monte Carlo. By introducing iterative Krylov subspace solvers and by multiplying by the Hamiltonian and overlap matrices as they are sampled, we remove the need to construct and store these matrices and thus bypass the most expensive steps of the stochastic reconfiguration and linear method optimization techniques. We demonstrate the effectiveness of this approach by using stochastic reconfiguration to optimize a correlator product state wave function with a Pfaffian reference for four example systems. In two examples on the two dimensional Fermionic Hubbard model, we study 16 and 64 site lattices, recovering energies accurate to 1$%$ in the smaller lattice and predicting particle-hole phase separation in the larger. In two examples involving an ab initio Hamiltonian, we investigate the potential energy curve of a symmetrically dissociated $4\ifmmode\times\else\texttimes\fi{}4$ hydrogen lattice as well as the singlet-triplet gap in free base porphin. In the hydrogen system we recover 98$%$ or more of the correlation energy at all geometries, while for porphin we compute the gap in a 24 orbital active space to within 0.02 eV of the exact result. The number of variational parameters in these examples ranges from $4\ifmmode\times\else\texttimes\fi{}{10}^{3}$ to $5\ifmmode\times\else\texttimes\fi{}{10}^{5}$.

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.

Recycling BiCG with an Application to Model Reduction

Science and engineering problems frequently require solving a sequence of dual linear systems. Besides having to store only a few Lanczos vectors, using the biconjugate gradient method (BiCG) to solve dual linear systems has advantages for specific applications. For example, using BiCG to solve the dual linear systems arising in interpolatory model reduction provides a backward error formulation in the model reduction framework. Using BiCG to evaluate bilinear forms---for example, in quantum Monte Carlo (QMC) methods for electronic structure calculations---leads to a quadratic error bound. Since our focus is on sequences of dual linear systems, we introduce recycling BiCG, a BiCG method that recycles two Krylov subspaces from one pair of dual linear systems to the next pair. The derivation of recycling BiCG also builds the foundation for developing recycling variants of other bi-Lanczos based methods, such as CGS, BiCGSTAB, QMR, and TFQMR. We develop an augmented bi-Lanczos algorithm and a modified two-term recurrence to include recycling in the iteration. The recycle spaces are approximate left and right invariant subspaces corresponding to the eigenvalues closest to the origin. These recycle spaces are found by solving a small generalized eigenvalue problem alongside the dual linear systems being solved in the sequence. We test our algorithm in two application areas. First, we solve a discretized partial differential equation (PDE) of convection-diffusion type. Such a problem provides well-known test cases that are easy to test and analyze further. Second, we use recycling BiCG in the iterative rational Krylov algorithm (IRKA) for interpolatory model reduction. IRKA requires solving sequences of slowly changing dual linear systems. We analyze the generated recycle spaces and show up to 70% savings in iterations. For our model reduction test problem, we show that solving the problem without recycling leads to (about) a 50% increase in runtime.