Design and implementation of parallel approximate inverse classes using OpenMP

High performance finite element approximate inverse preconditioning

PurposeTo propose novel parallel/distributed normalized explicit finite element (FE) approximate inverse preconditioning for solving sparse FE linear systems.Design/methodology/approachThe design of suitable methods was the main objective for which several families of the normalized approximate inverse, based on sparse normalized approximate factorization, are produced. The main motive for the derivation of the new normalized approximate inverse FE matrix algorithmic techniques is that they can be efficiently used in conjunction with normalized explicit preconditioned conjugate gradient (NEPCG) – type schemes on parallel and distributed systems. Theoretical estimates on the rate of convergence and computational complexity of the NEPCG method are also derived.FindingsApplication of the proposed method on a three‐dimensional boundary value problem is discussed and numerical results for uniprocessor systems along with speed‐ups and efficiency for multicomputer systems are given. These results tend to become optimum, which are in qualitative agreement with the theoretical results presented for uniprocessor and distributed memory systems, using message passing interface (MPI) communication library.Research limitations/implicationsFurther parallel algorithmic techniques will be investigated in order to improve the speed‐ups and the computational complexity of the parallel normalized explicit approximate inverse preconditioning.Originality/valueThe proposed parallel/distributed normalized explicit approximate inverse preconditioning, using approximate factorization and approximate inverse algorithms, is an efficient computational method that is valuable for computer scientists and for scientists and engineers in engineering computations.

During the recent decades, multigrid methods have been extensively used in order to solve large scale linear systems derived from the discretization of partial differential equations using the finite difference method. The purpose of this work is the derivation of hybrid heterogeneous schemes based on time implicit approximating schemes, namely backward differences and Crank–Nicolson method, with the multigrid method in conjunction with explicit approximate inverses for the efficient iterative solution of unsymmetric linear systems. Thus, a new class of smoothers, based on approximate inverses, has been derived. Moreover, a dynamic relaxation scheme is proposed based on the dynamic over - under relaxation (DOUR) algorithm. In addition, a new class of Krylov subspace iterative methods in conjunction with multigrid preconditioning based on explicit approximate inverses is used for solving initial value problems. Finally, the applicability of the new proposed hybrid heterogeneous implicit time - explicit approximate inverses combined with multigrid methods for parabolic PDEs is discussed by solving characteristic time-dependent problems in two dimensions and the numerical results are given.

A new class of Java multithreading based parallel approximate inverse preconditioning is introduced, for efficiently solving sparse arrow-type linear systems. The Parallel Explicit Preconditioned Biconjugate Conjugate Gradient - STAB method for shared memory systems is presented in order to examine the parallel behavior of this scheme using explicit approximate inverses as the suitable preconditioner. Design and implementation issues of Java's multithreading techniques are also discussed. The performance in terms of speedups and parallel efficiencies of the method, using Java multithreading, is illustrated by solving sparse arrow-type linear systems. Static and dynamic workload scheduling systems implemented in Java and the results of their use are presented and discussed.

During the last decades, research efforts have been focused on the derivation of effective preconditioned iterative methods. The preconditioned iterative methods are mainly categorized into implicit preconditioned methods and explicit preconditioned methods. In this manuscript we review implicit preconditioned methods, based on incomplete and approximate factorization, and explicit preconditioned methods, based on sparse approximate inverses and explicit approximate inverses. Modified Moore-Penrose conditions are presented and theoretical estimates for the sensitivity of the explicit approximate inverse matrix of the explicit preconditioned method are derived. Finally, the performance of the preconditioned iterative methods is illustrated by solving characteristic 2D elliptic problem and numerical results are given indicating a qualitative agreement with the theoretical estimates.

During the last decades, explicit approximate inverse preconditioning methods have been used for efficiently solving sparse linear systems on multiprocessor systems. The effectiveness of explicit approximate inverse preconditioning schemes relies on the use of efficient preconditioners that are close approximants to the coefficient matrix and are fast to compute in parallel. A new parallel computational technique is proposed for the parallelization of the explicit preconditioned conjugate gradient type method on a Graphics Processing Unit (GPU). The proposed parallel methods have been implemented using Compute Unified Device Architecture (CUDA) developed by NVIDIA. The inherently parallel operations between vectors and matrices involved in the explicit preconditioned biconjugate conjugate gradient type schemes exhibit significant amounts of loop-level parallelism because of the matrix–vector and the vector–vector products that can lead to high performance gain on the GPU systems, specifically designed for such computations. Finally, numerical results for the performance of the explicit preconditioned biconjugate conjugate gradient type method for solving characteristic two dimensional boundary value problems, using the finite difference method, on a massive multiprocessor interface on a GPU are presented. The CUDA implementation issues of the proposed method are also discussed.

Modeling of the heart ventricles is one of the most challenging tasks in soft tissue mechanics because cardiac tissue is a strongly anisotropic incompressible material with an active component of stress. In most current approaches with active force, the number of degrees of freedom (DOF) is limited by the direct method of solution of linear systems of equations. We develop a new approach for high-resolution heart models with large numbers of DOF by: (1) developing a hex-dominant finite element mixed formulation and (2) developing a Krylov subspace iterative method that is able to solve the system of linearized equations for saddle-point problems with active stress. In our approach, passive cardiac tissue is modeled as a hyperelastic, incompressible material with orthotropic properties, and mixed pressure-displacement finite elements are used to enforce incompressibility. Active stress is generated by a model with force dependence on length and velocity of muscle shortening. The ventricles are coupled to a lumped circulatory model. For efficient solution of linear systems, we use Flexible GMRES with a nonlinear preconditioner based on block matrix decomposition involving the Schur complement. Three methods for approximating the inverse of the Schur complement are evaluated: inverse of the pressure mass matrix; least squares commutators; and sparse approximate inverse. The sub-matrix corresponding to the displacement variables is preconditioned by a V-cycle of hybrid geometric-algebraic multigrid followed by correction with several iterations of GMRES preconditioned by sparse approximate inverse. The overall solver is demonstrated on a high-resolution two ventricle mesh based on a human anatomy with roughly 130K elements and 1.7M displacement DOF. Effectiveness of the numerical method for active contraction is shown. To the best of our knowledge, this solver is the first to efficiently model ventricular contraction using an iterative linear solver for the mesh size demonstrated and therefore opens the possibility for future very high-resolution models. In addition, several relatively simple benchmark problems are designed for a verification exercise to show that the solver is functioning properly and correctly solves the underlying mathematical model. Here, the output of the newly designed solver is compared to that of the mechanics component of Chaste ('Cancer, Heart and Soft Tissue Environment'). These benchmark tests may be used by other researchers to verify their newly developed methods and codes.

In this dissertation we develop numerical algorithms and software tools to facilitate parallel solutions of nonlinear programming (NLP) problems. In particular, we address large-scale, block-structured problems with an intrinsic decomposable configuration. These problems arise in a great number of engineering applications, including parameter estimation, optimal control, network optimization, and stochastic programming. The structure of these problems can be leveraged by optimization solvers to accelerate solutions and overcome memory limitations, and we propose variants to two classes of optimization algorithms: augmented Lagrangian (AL) schemes and Schur-complement interior-point methods. The convergence properties of augmented Lagrangian decomposition schemes like the alternating direction method of multipliers (ADMM) and progressive hedging (PH) are well established for convex optimization but convergence guarantees in non-convex settings are still poorly understood. In practice, however, ADMM and PH often perform satisfactorily in complex non-convex NLPs. In this work, we study connections between the method of multipliers (MM), ADMM, and PH to derive benchmarking metrics that explain why PH and ADMM work in practice. We illustrate the concepts using challenging dynamic optimization problems. Our exposition seeks to establish more formalism in benchmarking ADMM, PH, and AL schemes and to motivate algorithmic improvements. The effectiveness of nonlinear interior-point solvers for solving large-scale problems relies quite heavily on the solution of the underlying linear algebra systems. The schur-complement decomposition is very effective for parallelizing the solution of linear systems with modest coupling. However, for systems with large number of coupling variables the schur-complement method does not scale favorably. We implement an approach that uses a Krylov solver (GMRES) preconditioned with ADMM to solve block-structured linear systems that arise in the interior-point method. We show that this ADMM-GMRES approach overcomes the well-known scalability issues of Schur decomposition. One important drawback of using decomposition approaches like ADMM and PH is their convergence rate. Unlike Schur-complement interior-point algorithms that have super-linear convergence, augmented Lagrangian approaches typically exhibit linear and sublinear rates. We exploit connections between ADMM and the Schur-complement decomposition to derive an accelerated version of ADMM. Specifically, we study the effectiveness of performing a Newton-Raphson algorithm to compute multiplier estimates for augmented Lagrangian methods. We demonstrate using two-stage stochastic programming problems that our multiplier update achieves convergence in fewer iterations for MM on general nonlinear problems. In the case of ADMM, the newton update significantly reduces the number of subproblem solves for convex quadratic programs (QPs). Moreover, we show that using newton multiplier updates makes the method robust to the selection of the penalty parameter. Traditionally, state-of-the-art optimization solvers are implemented in low-level programming languages. In our experience, the development of decomposition algorithms in these frameworks is challenging. They present a steep learning curve and can slow the development and testing of new numerical algorithms. To mitigate these challenges, we developed PyNumero, a new open source framework implemented in Python and C++. The package seeks to facilitate development of optimization algorithms for large-scale optimization within a high-level programming environment while at the same time minimizing the computational burden of using Python. The efficiency of PyNumero is illustrated by implementing algorithms for problems arising in stochastic programming and optimal control. Timing results are presented for both serial and parallel implementations. Our computational studies demonstrate that with the appropriate balance between compiled code and Python, efficient implementations of optimization algorithms are achievable in these high-level languages.

Intelligent optimization algorithm (IOA) has been widely studied and applied to solve various optimization problems. When scholars improve IOA with mathematical methods, they also want to seek an effective method to implement algorithms with higher real time, especially for a complex problem. Parallel design is an effective method to improve the real time of IOA. Currently, the parallel programming based on open multi-processing (OpenMP) and compute unified device architecture (CUDA) are two popular methods. To find and develop a new IOA parallel method, in this paper, a parallel design and implementation method based on field programmable gate array (FPGA) is explored. In order to validate the proposed method, parallel genetic algorithm (GA) and parallel particle swarm optimization (PSO) algorithm are realized by the proposed method. Furthermore, the performance and advantage of the proposed FPGA-based parallel IOA method are tested by comparing with OpenMP-based parallel programming and CUDA-based parallel programming, the final results show that the proposed method with highest real-time performance in IOA parallel implementation. A case study by using FPGA-based parallel simulate annealing (SA) to address job shop scheduling problem (JSSP) to illustrate the proposed method has high potential in industrial applications.

The efficient numerical solution of the large linear systems of fractional differential equations is considered here. The key tool used is the short---memory principle. The latter ensures the decay of the entries of the inverse of the discretized operator, whose inverses are approximated here by a sequence of sparse matrices. On this ground, we propose to solve the underlying linear systems by these approximations or by iterative solvers using sequence of preconditioners based on the above mentioned inverses.

A hybrid parallel scheme for the ensemble square root filter (EnSRF) suitable for parallel assimilation of multiscale observations, including those from dense observational networks such as those of radar, is developed based on the domain decomposition strategy. The scheme handles internode communication through a message passing interface (MPI) and the communication within shared-memory nodes via Open Multiprocessing (OpenMP) threads. It also supports pure MPI and pure OpenMP modes. The parallel framework can accommodate high-volume remote-sensed radar (or satellite) observations as well as conventional observations that usually have larger covariance localization radii. The performance of the parallel algorithm has been tested with simulated and real radar data. The parallel program shows good scalability in pure MPI and hybrid MPI–OpenMP modes, while pure OpenMP runs exhibit limited scalability on a symmetric shared-memory system. It is found that in MPI mode, better parallel performance is achieved with domain decomposition configurations in which the leading dimension of the state variable arrays is larger, because this configuration allows for more efficient memory access. Given a fixed amount of computing resources, the hybrid parallel mode is preferred to pure MPI mode on supercomputers with nodes containing shared-memory cores. The overall performance is also affected by factors such as the cache size, memory bandwidth, and the networking topology. Tests with a real data case with a large number of radars confirm that the parallel data assimilation can be done on a multicore supercomputer with a significant speedup compared to the serial data assimilation algorithm.

New parallel computational techniques are introduced for the parallelization of explicit finite difference approximate inverse methods, using higher order approximation schemes, based on Portable Operating System Interface for UniX (POSIX) Threads, for multicore systems. Parallelization of the Optimized Banded Generalized Approximate Inverse Finite Element Matrix algorithm is achieved based on the concept of the fish bone approach with the use of a thread pool pattern. Theoretical estimates on speedups and efficiency are also presented. Finally, numerical results for the performance of the POBGAIFEM algorithm as well as the loop level and functional level parallelized Explicit Preconditioned Bi-conjugate Gradient Gradient-STAB method for solving two dimensional boundary value problems on multicore computer systems are presented.

In this paper we present new parallel inverse arrow-type matrix algorithms based on the concept of sparse factorization procedures, for computing explicitly exact and approximate inverses, on symmetric multiprocessor systems. The parallel implementation of the new inversion algorithms is discussed and numerical results are presented, using the simulation tool of Multi-Pascal.

The critical path of the hardware for the global elimination algorithm (GEA) is too long to meet the real-time constraints for high-end applications. In this paper, we propose a new parallel GEA and its corresponding architecture. By dividing candidate blocks into independent groups and finding the most probable candidates of each group in parallel, instead of sequentially searching within the whole search range, parallel design can be developed as an array of GEA processing elements with much shorter critical path. Besides, the GEA processing element is optimized to reduce 30% of the gates, and the 2D data reuse is organized to save 80% of the SRAM bandwidth, which also reduces a lot of power. Simulation results show that our implementation can achieve real time processing of D1 30 Hz video with search range as H[-64, +63.5] V[-32, +31.5] while the operating frequency is 70 MHz, and the gate count is 113 K. Compared with full search, our gate count is six times smaller under the same frequency, and the PSNR loss is at most 0.1-0.2 dB.

Computer Solution of Large Linear Systems