Iterative Solvers

A discrete Fourier transform (DFT)-based iterative method of moments (IMoM) algorithm is developed to provide an O(N tot ) computational complexity and memory storages for the efficient analysis of electromagnetic radiation/scattering from large phased arrays. Here, N tot is the total number of unknowns. Numerical results for both printed and free-standing dipole arrays are presented to validate the algorithm's efficiency and accuracy.

In this paper we present some computational strategies tailored for the Finite Element solution of large-scale flow problems in high performance computers. Reduced Integration techniques and edge-based data structures are studied. We also introduce new preconditioners for implicit schemes based on superedges clustering. The performance of these new computational strategies is evaluated on the Cray J90 and SGI/Origin 2000 parallel computers.

Writing efficient iterative solvers for irregular sparse matrices in High Performance Fortran (HPF) is hard. The locality in the computations is unclear, and for efficiency we use storage schemes that obscure any structure in the matrix. Moreover, the limited capabilities of HPF to distribute and align data structures make it hard to implement the desired distributions, or to indicate these such that the compiler recognizes the efficient implementation. We propose techniques to handle these problems. We combine strategies that have become popular in message-passing parallel programming, like mesh partitioning and splitting the matrix in local submatrices, with the functionality of HPF and HPF compilers, like the implicit handling of communication and distribution. The implementation of these techniques in HPF is not trivial, and we describe in detail how we propose to solve the problems. Our results demonstrate that very efficient implementations are possible. We indicate how some of the 'approved extensions' of HPF-2.0 can be used, but they do not solve all problems. For comparison we show the results for regular sparse matrices.

Writing efficient iterative solvers for irregular sparse matrices in High Performance Fortran (HPF) is hard. The locality in the computations is unclear, and for efficiency we use storage schemes that obscure any structure in the matrix. Moreover, the limited capabilities of HPF to distribute and align data structures make it hard to implement the desired distributions, or to indicate these such that the compiler recognizes the efficient implementation. We propose techniques to handle these problems. We combine strategies that have become popular in message-passing parallel programming, like mesh partitioning and splitting the matrix in local submatrices, with the functionality of HPF and HPF compilers, like the implicit handling of communication and distribution. The implementation of these techniques in HPF is not trivial, and we describe in detail how we propose to solve the problems. Our results demonstrate that very efficient implementations are possible. We indicate how some of the 'approved extensions' of HPF-2.0 can be used, but they do not solve all problems. For comparison we show the results for regular sparse matrices.

The Conjugate Gradient Squared (CGS) is an iterative method for solving nonsymmetric linear systems of equations. However, during the iteration large residual norms may appear, which may lead to inaccurate approximate solutions or may even deteriorate the convergence rate. Instead of squaring the Bi-CG polynomial as in CGS, we propose to consider products of two nearby Bi-CG polynomials which leads to generalized CGS methods, of which CGS is just a particular case. This approach allows the construction of methods that converge less irregularly than CGS and that improve on other convergence properties as well. Here, we are interested in a property that got less attention in literature: we concentrate on retaining the excellent approximation qualities of CGS with respect to components of the solution in the direction of eigenvectors associated with extreme eigenvalues. This property seems to be important in connection with Newton's scheme for nonlinear equations: our numerical experiments show that the number of Newton steps may decrease significantly when using a generalized CGS method as linear solver for the Newton correction equations.

New compact approximation schemes for the Laplace operator of 4th-and 6thorder are proposed. The schemes are based on a Padé approximation of the Taylor expansion for the discretized Laplace operator. The new schemes are compared with other finite difference approximations in several test cases. It is found that the new schemes are superior in performance and accuracy with respect to other methods.

Boundary Element Method frequency sweep analyses in acoustics are usually accompanied by a vast numerical cost of assembling and solving numerous linear systems. In that context, this work proposes a model order reduction technique to mitigate the resulting computational cost of such analyses. First, a series expansion of the Green's function BEM kernel is leveraged to construct a series of frequency independent matrices. Next, in a model order reduction way, the arising matrices are projected on a reduced basis utilizing a Galerkin projection. By this off-line matrix projection, both the assembly and the solution of the BEM full-size linear systems degenerate into assembling and solving a reduced system for all frequencies. Significant speed-up factors can, thus, be achieved for both operations. The projection basis employed in this model reduction scheme is developed through an Arnoldi algorithm for the BEM systems on a grid of master frequencies. The method is based on Krylov subspaces recycling, as the subspaces produced at master frequencies are recycled to approximate the surface distribution of the acoustic variables on the whole frequency range of interest. Utilizing Krylov subspaces facilitates as well the definition of a robust error estimator that indicates the quality of the reduced system. The performance of the proposed method is assessed for both an exterior and an interior problem for a simple and more complicated geometry respectively.

Boundary Element Method frequency sweep analyses in acoustics are usually accompanied by a vast numerical cost of assembling and solving numerous linear systems. In that context, this work proposes a model order reduction technique to mitigate the resulting computational cost of such analyses. First, a series expansion of the Green's function BEM kernel is leveraged to construct a series of frequency independent matrices. Next, in a model order reduction way, the arising matrices are projected on a reduced basis utilizing a Galerkin projection. By this off-line matrix projection, both the assembly and the solution of the BEM full-size linear systems degenerate into assembling and solving a reduced system for all frequencies. Significant speed-up factors can, thus, be achieved for both operations. The projection basis employed in this model reduction scheme is developed through an Arnoldi algorithm for the BEM systems on a grid of master frequencies. The method is based on Krylov subspaces recycling, as the subspaces produced at master frequencies are recycled to approximate the surface distribution of the acoustic variables on the whole frequency range of interest. Utilizing Krylov subspaces facilitates as well the definition of a robust error estimator that indicates the quality of the reduced system. The performance of the proposed method is assessed for both an exterior and an interior problem for a simple and more complicated geometry respectively.

We present an automated performance evaluation framework that enables an automated workflow for testing and performance evaluation of software libraries. Integrating this component into an ecosystem enables sustainable software development, as a community effort, via a web application for interactively evaluating the performance of individual software components. The performance evaluation tool is based exclusively on web technologies, which removes the burden of downloading performance data or installing additional software. We employ this framework for the Ginkgo software ecosystem, but the framework can be used with essentially any software project, including the comparison between different software libraries. The Continuous Integration (CI) framework of Ginkgo is also extended to automatically run a benchmark suite on predetermined HPC systems, store the state of the machine and the environment along with the compiled binaries, and collect results in a publicly accessible performance data repository based on Git. The Ginkgo performance explorer (GPE) can be used to retrieve the performance data from the repository, and visualizes it in a web browser. GPE also implements an interface that allows users The authors would like to thank the BSSw and the xSDK community efforts, in particular Mike Heroux and Lois Curfman McInnes, for promoting guidelines for a healthy software life cycle.

We present an automated performance evaluation framework that enables an automated workflow for testing and performance evaluation of software libraries. Integrating this component into an ecosystem enables sustainable software development, as a community effort, via a web application for interactively evaluating the performance of individual software components. The performance evaluation tool is based exclusively on web technologies, which removes the burden of downloading performance data or installing additional software. We employ this framework for the Ginkgo software ecosystem, but the framework can be used with essentially any software project, including the comparison between different software libraries. The Continuous Integration (CI) framework of Ginkgo is also extended to automatically run a benchmark suite on predetermined HPC systems, store the state of the machine and the environment along with the compiled binaries, and collect results in a publicly accessible performance data repository based on Git. The Ginkgo performance explorer (GPE) can be used to retrieve the performance data from the repository, and visualizes it in a web browser. GPE also implements an interface that allows users The authors would like to thank the BSSw and the xSDK community efforts, in particular Mike Heroux and Lois Curfman McInnes, for promoting guidelines for a healthy software life cycle.

We develop a new energy-aware methodology to improve the energy consumption of a task-parallel preconditioned Conjugate Gradient iterative solver on a Haswell-EP Intel Xeon. <br> This technique leverages the power-saving modes of the processor and the frequency range of the {\em userspace} Linux governor, modifying the CPU frequency for some operations. <br> We demonstrate that its application during the main operations of the PCG solver can reduce its energy consumption.

We investigate the benefits that an energyaware implementation of the runtime in charge of the concurrent execution of ILUPACK-a sophisticated preconditioned iterative solver for sparse linear systems-produces on the time-power-energy balance of the application. Furthermore, to connect the experimental results with the theory, we propose several simple yet accurate power models that capture the variations of average power that result from the introduction of the energy-aware strategies as well as the impact of the P-states into ILUPACK's runtime, at high accuracy, on two distinct platforms based on multicore technology from AMD and Intel.

Special thanks to Dr. Raymond C. Rumpf for his support and guidance throughout this research endeavor. I would also like to express my great appreciation to Dr. Rodrigo Romero for his instrumental help in developing the software necessary to complete this task. Assistance provided by Alejandro Martinez in the area of parallel computing was greatly appreciated.

Most efficient linear solvers use composable algorithmic components, with the most common model being the combination of a Krylov accelerator and one or more preconditioners. A similar set of concepts may be used for nonlinear algebraic systems, where nonlinear composition of different nonlinear solvers may significantly improve the time to solution. We describe the basic concepts of nonlinear composition and preconditioning and present a number of solvers applicable to nonlinear partial differential equations. We have developed a software framework in order to easily explore the possible combinations of solvers. We show that the performance gains from using composed solvers can be substantial compared with gains from standard Newton-Krylov methods.

Matrix Chain Multiplication plays a key role in the training of deep learning models. They also appear in physics, computer graphics, image processing, etc. Matrix Multiplications often cause a bottleneck in terms of performance and energy because of the heavy costs in computations and memory operations. While the runtime performance has been studied for years, significantly less effort has been expended into optimizing energy efficiency. Most processors weren't built with energy in mind [3], only performance. For a supercomputer, the cost of power is the most significant expense, and even with the most affordable energy, the cost of a MW ends up at $1M. Moreover, the Exascale Computing Project [4] defines a 'capable' exascale system as one that operates in a power envelope of 20-30 MW. Thus, reducing the energy cost of these types of computations is a major challenge. Matrix Chain Multiplication Problem: Given a sequence of matrices {A 1, A 2, ... A n } with sizes {P 0, P 1, ... P n } , compute ς =1 • Multiplication order can significantly impact the performance of the algorithm • The Optimal Parenthesization (OP_Count) algorithm in "Cormen et al. Introduction to Algorithms." outputs an order of matrix multiplications that minimizes the total number of operations. Energy Efficient GPU implementation • Energy consumption of a GPU can be broken down into 2 major parts: • Energy of the operations executing on the GPU • Energy of the GPU itself • Memory operations dominate the executed operations Our goal : optimize the total energy consumption for Matrix Chain Multiplications

Krylov subspace iterative solvers are often the method of choice when solving large sparse linear systems. At the same time, hardware accelerators such as graphics processing units (GPUs) continue to offer significant floating point performance gains for matrix and vector computations through easy-to-use libraries of computational kernels. However, as these libraries are usually composed of a well optimized but limited set of linear algebra operations, applications that use them often fail to reduce certain data communications, and hence fail to leverage the full potential of the accelerator. In this paper we target the acceleration of Krylov subspace iterative methods for GPUs, and in particular the BiCGSTAB solver, showing that significant improvement can be achieved by reformulating the method to reduce data-communications through applicationspecific kernels instead of using the generic BLAS kernels, e.g., as provided by NVIDIA's cuBLAS library, and by designing a GPU-specific sparse matrix-vector product (SpMV) kernel that is able to more efficiently use the GPU's computing power. Furthermore, we derive a model estimating the performance improvement, and use experimental data to validate the expected runtime savings. Considering that the derived implementation achieves significantly higher performance, we assert that similar optimizations addressing algorithm structure, as well as SpMV, are crucial for the subsequent development of high-performance GPU-accelerated Krylov subspace iterative methods.

The finite cell method is a highly flexible discretization technique for numerical analysis on domains with complex geometries. By using a non-boundary conforming computational domain that can be easily meshed, automatized computations on a wide range of geometrical models can be performed. Application of the finite cell method, and other immersed methods, to large real-life and industrial problems is often limited due to the conditioning problems associated with these methods. These conditioning problems have caused researchers to resort to direct solution methods, which significantly limit the maximum size of solvable systems. Iterative solvers are better suited for large-scale computations than their direct counterparts due to their lower memory requirements and suitability for parallel computing. These benefits can, however, only be exploited when systems are properly conditioned. In this contribution we present an Additive-Schwarz type preconditioner that enables efficient and parallel scalable iterative solutions of large-scale multi-level hp-refined finite cell analyses.

Many machine learning methods involve iterative optimization and are amenable to a variety of alternate formulations. Many currently popular formulations for some machine learning methods based on core operations that essentially correspond to sparse matrix-vector products. A reformulation using sparse matrix-matrix products primitives can potentially enable significant performance enhancement. Sampled Dense-Dense Matrix Multiplication (SDDMM) is a primitive that has been shown to be usable as a core component in reformulations of many machine learning factor analysis algorithms such as Alternating Least Squares (ALS), Latent Dirichlet Allocation (LDA), Sparse Factor Analysis (SFA), and Gamma Poisson (GaP). It requires the computation of the product of two input dense matrices but only at locations of the result matrix corresponding to nonzero entries in a sparse third input matrix. In this paper, we address the development of cuSDDMM, a multi-node GPU-accelerated implementation for SDDMM. We analyze the data reuse characteristics of SDDMM and develop a model-driven strategy for choice of tiling permutation and tilesize choice. cuSDDMM improves significantly (upto 4.6x) over the best currently available GPU implementation of SDDMM (in the BIDMach Machine Learning library).

Abstract. For the simulation of industrial sheet forming processes, the time discretisation is one of the important factors that determine the accuracy and efficiency of the algorithm. For relatively small models, the implicit time integration method is preferred, because of its inherent equilibrium check. For large models, the computation time becomes prohibitively large and, in practice, often explicit methods are used. In this contribution a strategy is presented that enables the application of implicit finite element simulations for large scale sheet forming analysis. Iterative linear equation solvers are commonly considered unsuitable for shell element mod-els. The condition number of the stiffness matrix is usually very poor and the extreme reduction of CPU time that is obtained in 3D bulk simulations is not reached in sheet forming simulations. Adding mass in an implicit time integration method has a beneficial effect on the condition num-ber. If mass scaling is used—like in ...

Special thanks to Dr. Raymond C. Rumpf for his support and guidance throughout this research endeavor. I would also like to express my great appreciation to Dr. Rodrigo Romero for his instrumental help in developing the software necessary to complete this task. Assistance provided by Alejandro Martinez in the area of parallel computing was greatly appreciated.

Matrix Chain Multiplication plays a key role in the training of deep learning models. They also appear in physics, computer graphics, image processing, etc. Matrix Multiplications often cause a bottleneck in terms of performance and energy because of the heavy costs in computations and memory operations. While the runtime performance has been studied for years, significantly less effort has been expended into optimizing energy efficiency. Most processors weren't built with energy in mind [3], only performance. For a supercomputer, the cost of power is the most significant expense, and even with the most affordable energy, the cost of a MW ends up at $1M. Moreover, the Exascale Computing Project [4] defines a 'capable' exascale system as one that operates in a power envelope of 20-30 MW. Thus, reducing the energy cost of these types of computations is a major challenge. Matrix Chain Multiplication Problem: Given a sequence of matrices {A 1, A 2, ... A n } with sizes {P 0, P 1, ... P n } , compute ς =1 • Multiplication order can significantly impact the performance of the algorithm • The Optimal Parenthesization (OP_Count) algorithm in "Cormen et al. Introduction to Algorithms." outputs an order of matrix multiplications that minimizes the total number of operations. Energy Efficient GPU implementation • Energy consumption of a GPU can be broken down into 2 major parts: • Energy of the operations executing on the GPU • Energy of the GPU itself • Memory operations dominate the executed operations Our goal : optimize the total energy consumption for Matrix Chain Multiplications

In parallel linear iterative solvers, sparse matrix vector multiplication (SpMxV) incurs irregular point-to-point (P2P) communications, whereas inner product computations incur regular collective communications. These P2P communications cause an additional synchronization point with relatively high message latency costs due to small message sizes. In these solvers, each SpMxV is usually followed by an inner product computation that involves the output vector of SpMxV. Here, we exploit this property to propose a novel parallelization method that avoids the latency costs and synchronization overhead of P2P communications. Our method involves a computational and a communication rearrangement scheme. The computational rearrangement provides an alternative method for forming input vector of SpMxV and allows P2P and collective communications to be performed in a single phase. The communication rearrangement realizes this opportunity by embedding P2P communications into global collective communication operations. The proposed method grants a certain value on the maximum number of messages communicated regardless of the sparsity pattern of the matrix. The downside, however, is the increased message volume and the negligible redundant computation. We favor reducing the message latency costs at the expense of increasing message volume. Yet, we propose two iterative-improvementbased heuristics to alleviate the increase in the volume through one-to-one task-to-processor mapping. Our experiments on two supercomputers, Cray XE6 and IBM BlueGene/Q, up to 2,048 processors show that the proposed parallelization method exhibits superior scalable performance compared to the conventional parallelization method.

Understanding the impact of soft errors on applications can be expensive. Often, it requires an extensive error injection campaign involving numerous runs of the full application in the presence of errors. In this paper, we present a novel approach to arrive at the ground truth-the true impact of an error on the final output-for iterative methods by observing a small number of iterations to learn deviations between normal and error-impacted execution. We develop a machine learning based predictor for three iterative methods to generate groundtruth results without running them to completion for every error injected. We demonstrate that this approach achieves greater accuracy than alternative prediction strategies, including three existing soft error detection strategies. We demonstrate the effectiveness of the ground truth prediction model in evaluating vulnerability and the effectiveness of soft error detection strategies in the context of iterative methods.

Krylov subspace iterative solvers are often the method of choice when solving large sparse linear systems. At the same time, hardware accelerators such as graphics processing units continue to offer significant floating point performance gains for matrix and vector computations through easy-to-use libraries of computational kernels. However, as these libraries are usually composed of a well optimized but limited set of linear algebra operations, applications that use them often fail to reduce certain data communications, and hence fail to leverage the full potential of the accelerator. In this paper, we target the acceleration of Krylov subspace iterative methods for graphics processing units, and in particular the Biconjugate Gradient Stabilized solver that significant improvement can be achieved by reformulating the method to reduce data-communications through application-specific kernels instead of using the generic BLAS kernels, e.g. as provided by NVIDIA’s cuBLAS library, and b...

Understanding the impact of soft errors on applications can be expensive. Often, it requires an extensive error injection campaign involving numerous runs of the full application in the presence of errors. In this paper, we present a novel approach to arrive at the ground truth-the true impact of an error on the final output-for iterative methods by observing a small number of iterations to learn deviations between normal and error-impacted execution. We develop a machine learning based predictor for three iterative methods to generate groundtruth results without running them to completion for every error injected. We demonstrate that this approach achieves greater accuracy than alternative prediction strategies, including three existing soft error detection strategies. We demonstrate the effectiveness of the ground truth prediction model in evaluating vulnerability and the effectiveness of soft error detection strategies in the context of iterative methods.