Number Of Right-hand Sides Research Articles

Portable CPU implementation of Wilson, Brillouin and Susskind fermions in lattice QCD

A modern Fortran implementation of three Dirac operators (Wilson, Brillouin, Susskind) in lattice QCD is presented, based on OpenMP shared-memory parallelization and SIMD pragmas. The main idea is to apply a Dirac operator to Nv vectors simultaneously, to ease the memory bandwidth bottleneck. All index computations are left to the compiler and maximum weight is given to portability and flexibility. The lattice volume, NxNyNzNt, the number of colors, Nc, and the number of right-hand sides, Nv, are parameters defined at compile time. Several memory layout options are compared. The code performs well on modern many-core architectures (480 Gflop/s, 880 Gflop/s, and 780 Gflop/s with Nv=12 for the three operators in single precision on a 72-core KNL processor, a 2×24-core Skylake node yields similar results). Explicit run-time tests with CG/BiCGstab inverters confirm that the memory layout is relevant for the KNL, but less so for the Skylake architecture. The ancillary code distribution contains all routines, including the single, double, and mixed precision Krylov space solvers, to render it self-contained and ready-to-use. Program summaryProgram title: testknl_main.f90CPC Library link to program files:https://doi.org/10.17632/7km9wcspn3.1Licensing provisions: Creative Commons by 4.0Programming language: fortran 2008Supplementary material: Detailed description of each module of the code distribution.Nature of problem: Solving the Dirac equation Du=b for three choices of the Dirac operator D relevant to Lattice QCD. The choices are known as Wilson fermions, Brillouin fermions and Susskind or staggered fermions, respectively.Solution method: Standard solution through sparse iterative solvers CG and BiCGstab, albeit with several right-hand vectors packed into a right-hand matrix. The number of colors, Nc, the number of right-hand sides, Nv, and the box lengths, Nx, Ny, Nz and Nt, can be chosen freely.Additional comments including restrictions and unusual features: The code allows the user to compare several memory layout options for the block of right-hand sides, that is several orderings of the slots pertinent to Nc, Nv and possibly the Dirac spinor, in the array representation of u and b. No machine dependent optimizations are used, only standard OpenMP shared-memory parallelization and SIMD pragmas are used.

Flexible subspace iteration with moments for an effective contour integration‐based eigensolver

AbstractContour integration schemes are a valuable tool for the solution of difficult interior eigenvalue problems. However, the solution of many large linear systems with multiple right hand sides may prove a prohibitive computational expense. The number of right hand sides, and thus, computational cost may be reduced if the projected subspace is created using multiple moments. In this work, we explore heuristics for the choice and application of moments with respect to various other important parameters in a contour integration scheme. We provide evidence for the expected performance, accuracy, and robustness of various schemes, showing that good heuristic choices can provide a scheme featuring good properties in all three of these measures.

The method of polarized traces for the 3D Helmholtz equation

We have developed a fast solver for the 3D Helmholtz equation, in heterogeneous, constant density, acoustic media, in the high-frequency regime. The solver is based on the method of polarized traces, a layered domain-decomposition method, where the subdomains are connected via transmission conditions prescribed by the discrete Green’s representation formula and artificial reflections are avoided by enforcing nonreflecting boundary conditions between layers. The method of polarized traces allows us to consider only unknowns at the layer interfaces, reducing the overall cost and memory footprint of the solver. We determine that polarizing the wavefields in this manner yields an efficient preconditioner for the reduced system, whose rate of convergence is independent of the problem frequency. The resulting preconditioned system is solved iteratively using generalized minimum residual, where we never assemble the reduced system or preconditioner; rather, we implement them via solving the Helmholtz equation locally within the subdomains. The method is parallelized using Message Passing Interface and coupled with a distributed linear algebra library and pipelining to obtain an empirical on-line runtime [Formula: see text], where [Formula: see text] is the total number of degrees of freedom, [Formula: see text] is the number of subdomains, and [Formula: see text] is the number of right-hand sides (RHS). This scaling is favorable for regimes in which the number of sources (distinct RHS) is large, for example, enabling large-scale implementations of frequency-domain full-waveform inversion.

Development of the Block BiCGGR2 method for linear systems with multiple right-hand sides

Block Krylov subspace methods are efficient iterative methods for solving linearsystems with multiple right-hand sides. The Block BiCGGR method is one of the Block Krylov subspace methods, and can generate highly accurate approximatesolutions. However, when the number of right-hand sides is large, this method becomes numerically unstable. In this paper, the Block BiCGGR2 method is proposed. This method is a numerically stable and highly accurate method than the Block BiCGGR method and the Block BiCGSTAB method. Numerical results illustrate the robustness and the accuracy of the approximate solutions of the Block BiCGGR2 method.

Accuracy improvement of the Block BiCGSTAB method for linear systems with multiple right-hand sides by group-wise updating technique

The Block BiCGSTAB method is an efficient method for solving linear systems with multiple right-hand sides. However, when the number of right-hand sides is large, this method may not generate high accuracyapproximate solution due to an error in the computation of matrix-matrix multiplication. In this paper, in order to improve the accuracy of the approximate solution, the recursions of the Block BiCGSTAB method are reconstructed by using the group-wise updating technique. Moreover, the convergence property of the proposed method is also improved when the number of right-hand sides is large.

Математические проблемы организации метаболических путей из биохимических реакций

The problem of metabolic pathway building from biochemical reactions is considered. For the given pair “main substrate – main product” it is necessary to select the reactions performing this conversion and find all the attendant substrates and products. The approach is based on a system of stoichiometric equations for reaction rates. The mathematical problems are as follows: 1) a number of right-hand sides are initially unknown, 2) a part of the system with the known right-hand sides has less number of equations than variables, and 3) presence of irreversible reactions imposes restrictions on signs of the corresponding rates in the form of inequalities. The method is described enabling to solve these problems. It is found that many restrictions are represented as series of parallel hyperplanes in the space of variables. This fact makes it possible to eliminate many restrictions and to find out a big number of fixed variables the values of which, thus, are found before the complete solving of the whole problem. The search of unnecessary restrictions and fixed variables is of the form of specific iterations due to which the problem is substantially simplified. Examples are given which display effectiveness of the approach described.

A review of block Krylov subspace methods for multisource electromagnetic modelling

Practical applications of controlled-source electromagnetic (EM) modelling require solutions for multiple sources at several frequencies, thus leading to a dramatic increase of the computational cost. In this paper, we present an approach using block Krylov subspace solvers that are iterative methods especially designed for problems with multiple right-hand sides (RHS). Their main advantage is the shared subspace for approximate solutions, hence, these methods are expected to converge in less iterations than the corresponding standard solver applied to each linear system. Block solvers also share the same preconditioner, which is constructed only once. Simultaneously computed block operations have better utilization of cache due to the less frequent access to the system matrix. In this paper, we implement two different block solvers for sparse matrices resulting from the finite-difference and the finite-element discretizations, discuss the computational cost of the algorithms and study their dependence on the number of RHS given at once. The effectiveness of the proposed methods is demonstrated on two EM survey scenarios, including a large marine model. As the results of the simulations show, when a powerful preconditioning is employed, block methods are faster than standard iterative techniques in terms of both iterations and time.

3D Frequency-Domain Seismic Inversion with Controlled Sloppiness

Seismic waveform inversion aims at obtaining detailed estimates of subsurface medium parameters, such as the spatial distribution of soundspeed, from multiexperiment seismic data. A formulation of this inverse problem in the frequency domain leads to an optimization problem constrained by a Helmholtz equation with many right-hand sides. Application of this technique to industry-scale problems faces several challenges: First, we need to solve the Helmholtz equation for high wave numbers over large computational domains. Second, the data consist of many independent experiments, leading to a large number of PDE solves. This results in high computational complexity both in terms of memory and CPU time as well as input/output costs. Finally, the inverse problem is highly nonlinear and a lot of art goes into preprocessing and regularization. Ideally, an inversion needs to be run several times with different initial guesses and/or tuning parameters. In this paper, we discuss the requirements of the various components (PDE solver, optimization method, \dots) when applied to large-scale three-dimensional seismic waveform inversion and combine several existing approaches into a flexible inversion scheme for seismic waveform inversion. The scheme is based on the idea that in the early stages of the inversion we do not need all the data or very accurate PDE solves. We base our method on an existing preconditioned Krylov solver (CARP-CG) and use ideas from stochastic optimization to formulate a gradient-based (quasi-Newton) optimization algorithm that works with small subsets of the right-hand sides and uses inexact PDE solves for the gradient calculations. We propose novel heuristics to adaptively control both the accuracy and the number of right-hand sides. We illustrate the algorithms on synthetic benchmark models for which significant computational gains can be made without being sensitive to noise and without losing the accuracy of the inverted model.

An Effective Method for Parameter Estimation with PDE Constraints with Multiple Right-Hand Sides

Often, parameter estimation problems of parameter-dependent PDEs involve multiple right-hand sides. The computational cost and memory requirements of such problems increase linearly with the number of right-hand sides. For many applications this is the main bottleneck of the computation. In this paper we show that problems with multiple right-hand sides can be reformulated as stochastic programming problems by combining the right-hand sides into a few „simultaneous” sources. This effectively reduces the cost of the forward problem and results in problems that are much cheaper to solve. We discuss two solution methodologies: namely sample average approximation and stochastic approximation. To illustrate the effectiveness of our approach we present two model problems, direct current resistivity and seismic tomography.

A block IDR( [formula omitted]) method for nonsymmetric linear systems with multiple right-hand sides

The IDR( s ) based on the induced dimension reduction (IDR) theorem, is a new class of efficient algorithms for large nonsymmetric linear systems. IDR(1) is mathematically equivalent to BiCGStab at the even IDR(1) residuals, and IDR( s ) with s > 1 is competitive with most Bi-CG based methods. For these reasons, we extend the IDR( s ) to solve large nonsymmetric linear systems with multiple right-hand sides. In this paper, a variant of the IDR theorem is given at first, then the block IDR( s ), an extension of IDR( s ) based on the variant IDR( s ) theorem, is proposed. By analysis, the upper bound on the number of matrix-vector products of block IDR( s ) is the same as that of the IDR( s ) for a single right-hand side in generic case, i.e., the total number of matrix-vector products of IDR( s ) may be m times that of of block IDR( s ), where m is the number of right-hand sides. Numerical experiments are presented to show the effectiveness of our proposed method.

Deflation for inversion with multiple right-hand sides in QCD

Most calculations in lattice Quantum Chromodynamics (QCD) involve the solution of a series of linear systems of equations with exceedingly large matrices and a large number of right hand sides. Iterative methods for these problems can be sped up significantly if we deflate approximations of appropriate invariant spaces from the initial guesses. Recently we have developed eigCG, a modification of the Conjugate Gradient (CG) method, which while solving a linear system can reuse a window of the CG vectors to compute eigenvectors almost as accurately as the Lanczos method. The number of approximate eigenvectors can increase as more systems are solved. In this paper we review some of the characteristics of eigCG and show how it helps remove the critical slowdown in QCD calculations. Moreover, we study scaling with lattice volume and an extension of the technique to nonsymmetric problems.

Acceleration of free-space discrete body of revolution codes by exploiting circulant submatrices

Discrete body of revolution (DBOR) codes can accelerate moment method (MoM) solutions for scattering geometries possessing discrete circular symmetry. Traditional DBOR formulations require modifications to the incident field, such as decomposition into circular modes. While appropriate for methods employing direct solvers; approaches employing iterative solvers are disadvantaged by this dramatic increase in the number of right-hand sides that must be solved. In this correspondence, we show how the block circulant property of the impedance matrix can be exploited to avoid incident field modification while compressing matrix storage and accelerating computation of the iterative solver's associated matrix-vector product above that of the typical iterative MoM approach.

Estsol—a subroutine library for establishment and solution of linear equation systems

One of the crucial steps in the solution of problems with Finite Element Methods is establishment and solution of systems of equations. A numerous number of papers have been presented concerning different techniques of solution of equation systems such as direct methods, Mondkar and Powell (1974), Wilson (1978), Kamel and McCabe (1978), Argyris and Bronlund (1975), and iterative methods such as conjugate gradient, Axelsson and Gustafsson (1980), Gustafsson and Lindskog (1984), and viscous relaxations, Hughes et al. (1983), Zienkiewicz (1984). However, the user of these methods normally has to write a number of subroutines for the establishment and solution which are different for each method. The presented subroutine library will establish and solve a linear system of equations for a number of right hand sides using almost the same data for all methods. It provides the user of an application program with the possibility to select the solution method which is best suited concerning efficiency and accuracy for the problem at hand.