Scalable parallel sparse LU factorization methods on shared memory multiprocessors

Parallel sparse LU factorization on different message passing platforms

Application demands and grand challenges in numerical simulation require for both highly capable computing platforms and efficient numerical solution schemes. Power constraints and further miniaturization of modern and future hardware give way for multi- and manycore processors with increasing fine-grained parallelism and deeply nested hierarchical memory systems -- as already exemplified by recent graphics processing units. Accordingly, numerical schemes need to be adapted and re-engineered in order to deliver scalable solutions across diverse processor configurations. Portability of parallel software solutions across emerging hardware platforms is another challenge. This work investigates multi-coloring and re-ordering schemes for block Gauß-Seidel methods and, in particular, for incomplete LU factorizations with and without fill-ins. We consider two matrix re-ordering schemes that deliver flexible and efficient parallel preconditioners. The general idea is to generate block decompositions of the system matrix such that the diagonal blocks are diagonal itself. In such a way, parallelism can be exploited on the block-level in a scalable manner. Our goal is to provide widely applicable, out-of-the-box preconditioners that can be used in the context of finite element solvers. We propose a new method for anticipating the fill-in pattern of ILU($p$) schemes which we call the power($q$)-pattern method. This method is based on an incomplete factorization of the system matrix $A$ subject to a predetermined pattern given by the matrix power $|A|^(p+1)$ and its associated multi-coloring permutation $. We prove that the obtained sparsity pattern is a superset of our modified ILU($p$) factorization applied to pi A pi^(-1). As a result, this modified ILU($p$) applied to multi-colored system matrix has no fill-ins in its diagonal blocks. This leads to an inherently parallel execution of triangular ILU($p$) sweeps. In addition, we describe the integration of the preconditioners into the HiFlow$^3$ open-source finite element package that provides a portable software solution across diverse hardware platforms. On this basis, we conduct performance analysis across a variety of test problems on multi-core CPUs and GPUs that proves efficiency, scalability and flexibility of our approach. Our preconditioners achieve a solver acceleration by a factor of up to 1.5, 8 and 85 for three different test problems. The GPU versions of the preconditioned solver are by a factor of up to 4 faster than an OpenMP parallel version on eight cores.

Supernode pivoting for unsymmetric matrices coupled with supernode partitioning and asynchronous computation can achieve high gigaflop rates for parallel sparse LU factorization on shared memory parallel computers. The progress in weighted graph matching algorithms helps to extend these concepts further and prepermutation of rows is used to place large matrix entries on the diagonal. Supernode pivoting allows dynamical interchanges of columns and rows during the factorization process. The BLAS-3 level efficiency is retained. An enhanced left-right looking scheduling scheme is uneffected and results in good speedup on SMP machines without increasing the operation count. These algorithms have been integrated into the recent unsymmetric version of the PARDISO solver. Experiments demonstrate that a wide set of unsymmetric linear systems can be solved and high performance is consistently achieved for large sparse unsymmetric matrices from real world applications.

Dense LU factorization has a high ratio of computation to communication and, as evidenced by the High Performance Linpack (HPL) benchmark, this property makes it scale well on most parallel machines. Nevertheless, the standard algorithm for this problem has non-trivial dependence patterns which limit parallelism, and local computations require large matrices in order to achieve good single processor performance. We present an alternative programming model for this type of problem, which combines UPC's global address space with lightweight multithreading. We introduce the concept of memory-constrained lookahead where the amount of concurrency managed by each processor is controlled by the amount of memory available. We implement novel techniques for steering the computation to optimize for high performance and demonstrate the scalability and portability of UPC with Teraflop level performance on some machines, comparing favourably to other state-of-the-art MPI codes.

A new method for sparse LU factorization is presented that combines a column pre-ordering strategy with a right-looking unsymmetric-pattern multifrontal numerical factorization. The column ordering is selected to give a good a priori upper bound on fill-in and then refined during numerical factorization (while preserving the bound). Pivot rows are selected to maintain numerical stability and to preserve sparsity. The method analyzes the matrix and automatically selects one of three pre-ordering and pivoting strategies. The number of nonzeros in the LU factors computed by the method is typically less than or equal to those found by a wide range of unsymmetric sparse LU factorization methods, including left-looking methods and prior multifrontal methods.

QR and LU decompositions are the most important matrix decomposition algorithms. Many studies work on accelerating these algorithms by FPGA or ASIC in a case by case style. In this paper, we propose a unified framework for the matrix decomposition algorithms, combining three QR decomposition algorithms and LU algorithm with pivoting into a unified linear array structure. The QR and LU decomposition algorithms exhibit the same two-level loop structure and the same data dependency. Utilizing the similarities in loop structure and data dependency of matrix decomposition, we unify a fine-grained algorithm for all four matrix decomposition algorithms. Furthermore, we present a unified coprocessor structure with a scalable linear array of processing elements (PEs), in which four types of PEs are same in the structure of memory channels and PE connections, but the only difference exists in the internal structure of data path. Our unified co-processor, which is IEEE 32-bit floating-point precision, is implemented and mapped onto a Xilinx Virtex5 FPGA chip. Experimental results show that our co-processors can achieve speedup of 2.3 to 14.9 factors compared to a Pentium Dual CPU with double SSE threads.

LU factorization for sparse matrices is an important computing step for many engineering and scientific problems such as circuit simulation. There have been many efforts toward parallelizing and scaling this algorithm, which include the recent efforts targeting the GPUs. However, it is still challenging to deploy a complete sparse LU factorization workflow on a GPU due to high memory requirements and data dependencies. In this paper, we propose the first complete GPU solution for sparse LU factorization. To achieve this goal, we propose an out-of-core implementation of the symbolic execution phase, thus removing the bottleneck due to large intermediate data structures. Next, we propose a dynamic parallelism implementation of Kahn's algorithm for topological sort on the GPUs. Finally, for the numeric factorization phase, we increase the parallelism degree by removing the memory limits for large matrices as compared to the existing implementation approaches. Experimental results show that compared with an implementation modified from GLU 3.0, our out-of-core version achieves speedups of 1.13--32.65X. Further, our out-of-core implementation achieves a speedup of 1.2--2.2 over an optimized unified memory implementation on the GPU. Finally, we show that the optimizations we introduce for numeric factorization turn out to be effective.

Virtual earthquake engineering laboratory with physics-based degrading materials on parallel computers

In this paper we address the problem of matrix factorization on compressively-sampled measurements which are obtained by random projections. While this approach improves the scalability of matrix factorization, its performance is not satisfactory. We present a matrix co-factorization method where compressed measurements and a small number of uncompressed measurements are jointly decomposed, sharing a factor matrix. We evaluate the performance of three matrix factorization methods in terms of Cramer-Rao bounds, including: (1) matrix factorization on uncompressed data (MF); (2) matrix factorization on compressed data (CS-MF); (3) matrix co-factorization on compressed and uncompressed data (CS-MCF). Numerical experiments demonstrate that CS-MCF improves the performance of CS-MF, emphasizing the useful behavior of exploiting side information (a small number of uncompressed measurements).

The 1990s are seeing the explosive growth of the Internet and Web-based information sharing and dissemination systems. The Internet is also showing a potential of forming of a supercomputing resource out of networked computers. Parallel computing on the Internet often works in a machine-centric “pull” execution model. That is, a coordinator machine maintains a pool of tasks and distributes the tasks to other participants on demands. This paper proposes a novel mobile agent based “push” methodology from the perspective of applications. In the method, users declare their computation-bound jobs as autonomous agents. The computational agents will roam on the Internet to find servers to run. Since the agents can be programmed to satisfy their goals, even if they move and lose contact with their creators, they can survive intermittent or unreliable network connection. During their lifetime, the agents can also move themselves autonomously from one machine to another for load balancing, enhancing data locality, and tolerating faults. We present an agent-oriented programming and resource brokerage infrastructure, TRAVELER, in support of global parallel computing. The TRAVELER provides a mechanism for clients to wrap their parallel applications as mobile agents. The agents are dispatched to a resource broker. The broker forms a parallel virtual machine atop available servers to execute the agents. TRAVELER relies on an integrated distributed shared array (DSA) runtime system to support inter-agent communication and synchronization on clusters of servers. We demonstrate the feasibility of the TRAVELER in parallel sorting and LU factorization problems.

We present algorithms to determine the number of nonzeros in each row and column of the factors of a sparse matrix, for both the QR factorization and the LU factorization with partial pivoting. The algorithms use only the nonzero structure of the input matrix, and run in time nearly linear in the number of nonzeros in that matrix. They may be used to set up data structures or schedule parallel operations in advance of the numerical factorization. The row and column counts we compute are upper bounds on the actual counts. If the input matrix is strong Hall and there is no coincidental numerical cancellation, the counts are exact for QR factorization and are the tightest bounds possible for LU factorization. These algorithms are based on our earlier work on computing row and column counts for sparse Cholesky factorization, plus an efficient method to compute the column elimination tree of a sparse matrix without explicitly forming the product of the matrix and its transpose.

The sparse matrix solver has become the bottleneck in a Simulation Program with Integrated Circuit Emphasis circuit simulator. It is difficult to parallelize the sparse matrix solver because of the high data dependence during the numerical LU factorization. In this brief, a parallel LU factorization algorithm is developed on shared-memory computers with multicore central processing units, based on KLU algorithms. An Elimination Scheduler (EScheduler) is proposed to represent the data dependence during the LU factorization. Based on the EScheduler, the parallel tasks are scheduled in two modes to achieve a high level of concurrence, i.e., cluster mode and pipeline mode . The experimental results on 26 circuit matrices reveal that the developed algorithm can achieve speedup of 1.18-4.55t (on geometric average), as compared with KLU, with 1-8 threads. The result analysis indicates that for different data dependence, different parallel strategies should be dynamically selected to obtain optimal performance.

Recent trends have led hardware manufacturers to place multiple processing cores on a single chip, making parallel programming the intended way of taking advantage of the increased processing power. However, bringing concurrency to average programmers is considered to be one of the major challenges in computer science today. The difficulty lies not only in writing correct parallel programs, but also in achieving the required efficiency and performance. For parallel programs, performance is not only about obtaining low execution times on a fixed number of cores, but also about maintaining efficiency as the number of available cores is increased. Ideally, programmers should have in their toolkit techniques that can be used when designing parallel programs, before any code is available for testing on production hardware, and which are then able to predict scalability once the program is implemented. Existing methods are either unreliable at predicting scalability, such as the case of disjoint-access parallelism, or do not apply to lock-based programs, which currently make up a large part of existing concurrent programs. Furthermore, using some of these techniques is so complicated that it outweighs the time required to implement, debug and test the program on real hardware. In this thesis we study the problem of predicting the scalability of concurrent programs without implementing them. This allows programmers in the design phase of a concurrent algorithm to choose only one or a few promising solutions that will be implemented, debugged and tested on production hardware. We first consider disjoint-access parallelism, an existing property that applies only to a very restricted class of programs. After an extensive practical evaluation spanning across a variety of scenarios, we find it to be ineffective at predicting scalability. For predicting the scalability of more general concurrent algorithms, we propose the obstruction degree, a new scalability metric based on the consistency requirements of algorithms. It applies to programs using locks, invalidation primitives and transactional memory. Our metric allows programmers to compare two given algorithms as well as predict their scalability limit, the maximum number of processors to which they can scale, thus allowing programmers to choose the appropriate size hardware for running their programs. We also examine the composition of relaxed memory transactions in order to combine the ease of programming offered by transactional memory with the increased scalability of transactions that circumvent the traditional transactional model. We present outheritance, a property we show to be both necessary and sufficient for ensuring the correct composition of relaxed transactions, and we show how to calculate the obstruction degree of compositions that use this new property. We use outheritance to build OE-STM, a new software transactional memory algorithm having elastic transactions that correctly compose.

In a semiconductor market dominated by portable consumer applications, embedded flash memory technology has experienced a rapid diffusion. It is now considered the preferred solid-state memory solution for its non-volatile characteristics, high read and write speeds and scalability properties. As technology scales down in the nanometer range, new accurate physical tools should be made available to circuit designers, to support the development and optimization of high-voltage circuit blocks. A surface potential-based model for the flash memory cell has been developed with the purpose of providing a comprehensive physical understanding of the device operation, suitable for performance optimization in memory circuit design. An accurate validation methodology also takes into account charge balance effects on the isolated floating gate node and parasitic couplings inside and between the memory cells. The compact model supports DC, AC and transient analyses, including program/erase bias scalability, temperature effects, process corners and statistical variations. The results have been compared to Technology Computer-Aided Design (TCAD) simulations demonstrating that short channel effects, overlap capacitances and velocity saturation dominate over the intrinsic behaviour of the cell in ultrascaled devices. The approach includes drain disturb and memory endurance degradation models due to oxide aging. These effects are becoming dominant in ultrascaled devices. The model has been implemented using the Verilog-A language for portability into common circuit simulators. Validation has been performed on measurement results of test structures integrated in a 65nm derivative NOR CMOS technology. The compact model development has been based on a rigorous modeling approach combining conventional TCAD simulation tools with physically-based analyses. A new TCAD tool has been proposed for the investigation of advanced quantum effects, band structure models, quantum tunneling and multiphonon-assisted charge trapping effects in dielectrics. The effects of charge trapping in oxide layers and Si/SiO2 interfaces have been studied, specifically focusing on flash technology, where high voltage biases represent a major issue for dielectric degradation. A multiphonon-assisted model has been coupled with a Poisson-Schrödinger quantum solver. A novel impedance calculation method has been applied to the analysis of DC and AC MOS characteristics. This approach permits the physical modeling of trap filling, frequency response and device electrostatics. Transient effects of trap filling and trap-assisted tunneling through the gate have also been investigated. The adoption of such a multilevel approach permits to apply the methodology to flash memory cells. This enabled the investigation of the role of defects on electrostatics and program/erase efficiency reduction. The flash compact model has been applied in a process development and IC memory design perspective. Technology development requires a profound understanding of trade-offs in flash devices, which affect DC, transient and long-term performances. The design, integration and characterization of a 40 KB memory sector for smart card applications has been performed to demonstrate the capabilities of the compact model.

We present parallel and sequential dense QR factorization algorithms that are both optimal (up to polylogarithmic factors) in the amount of communication they perform and just as stable as Householder QR. We prove optimality by deriving new lower bounds for the number of multiplications done by “non-Strassen-like” QR, and using these in known communication lower bounds that are proportional to the number of multiplications. We not only show that our QR algorithms attain these lower bounds (up to polylogarithmic factors), but that existing LAPACK and ScaLAPACK algorithms perform asymptotically more communication. We derive analogous communication lower bounds for LU factorization and point out recent LU algorithms in the literature that attain at least some of these lower bounds. The sequential and parallel QR algorithms for tall and skinny matrices lead to significant speedups in practice over some of the existing algorithms, including LAPACK and ScaLAPACK, for example, up to 6.7 times over ScaLAPACK. A performance model for the parallel algorithm for general rectangular matrices predicts significant speedups over ScaLAPACK.