Intel MKL Research Articles

Scale-Free Sparse Matrix-Vector Multiplication on Many-Core Architectures

Sparse matrix-vector multiplication (SpMV) is one of the most important kernels for many applications. In this paper, we study the implementation of SpMV for scale-free matrices on many-core architectures including graphic processing units and Xeon Phi coprocessors. We first propose a hardware oblivious implementation for heterogeneous many-core processors using OpenCL. Our OpenCL implementation uses a novel SpMV format called hybrid COO+CSR (HCC), which employs 2-D jagged partitioning to balance the workload among a large number of cores and improve the data locality. Moreover, the OpenCL implementation is designed to be parametric, which allows systematic performance tuning. We conduct experiments to evaluate the efficiency of our hardware oblivious implementation. Experiments show that it achieves comparable performance to the Intel MKL and state-of-the-art OpenCL-based ViennaCL library implementation. Although the OpenCL implementation provides functional portability for heterogeneous systems, it fails to take advantage of the low-level architectural features. To further improve the performance, we propose a hardware conscious implementation using the native parallel programming language. We use the Xeon Phi platform as a case study. In our hardware conscious implementation, we ensure that the HCC format efficiently utilizes the vector process units on Xeon Phi by employing low-level intrinsics, and improve the overall performance through locality-aware block mapping, and intrablock tiling. Experiments using a wide range of representative scale-free matrices demonstrate that compared with the OpenCL-based hardware oblivious implementation, the hardware conscious implementation achieves $2.2\boldsymbol {\times }$ speedup on average. Compared with MKL, the hardware conscious implementation achieves $3.1\boldsymbol {\times }$ speedup on Xeon Phi.

Asynchronous Task-Based Polar Decomposition on Single Node Manycore Architectures

This paper introduces the first asynchronous, task-based formulation of the polar decomposition and its corresponding implementation on manycore architectures. Based on a new formulation of the iterative QR dynamically-weighted Halley algorithm (QDWH) for the calculation of the polar decomposition, the proposed implementation replaces the original and hostile LU factorization for the condition number estimator by the more adequate QR factorization to enable software portability across various architectures. Relying on fine-grained computations, the novel task-based implementation is also capable of taking advantage of the identity structure of the matrix involved during the QDWH iterations, which decreases the overall algorithmic complexity. Furthermore, the artifactual synchronization points have been weakened compared to previous implementations, unveiling look-ahead opportunities for better hardware occupancy. The overall QDWH-based polar decomposition can then be represented as a directed acyclic graph (DAG), where nodes represent computational tasks and edges define the inter-task data dependencies. The StarPU dynamic runtime system is employed to traverse the DAG, to track the various data dependencies and to asynchronously schedule the computational tasks on the underlying hardware resources, resulting in an out-of-order task scheduling. Benchmarking experiments show significant improvements against existing state-of-the-art high performance implementations (i.e., Intel MKL and Elemental) for the polar decomposition on latest shared-memory vendors' systems (i.e., Intel Haswell/Broadwell/Knights Landing, NVIDIA K80/P100 GPUs and IBM Power8), while maintaining high numerical accuracy.

R PROFILING: STRATEGIES FORPERFORMANCE IMPROVEMENTS IN R APPLICATION

This paper is in continuation with my efforts on improving performance of R application. In my previous paper titled “A Study on CRAN R and MRAN R Interpreters”, I have demonstrated performance improvements by Microsoft R over core R out of the box. Microsoft R uses Intel MKL library to make R applications multithreaded. In this paper, I have highlighted some APIs using which performance will increase manifolds. This replacement requires little to no overhead in code change thus very effective to use by legacy apps.

GPU-Accelerated LOBPCG Method with Inexact Null-Space Filtering for Solving Generalized Eigenvalue Problems in Computational Electromagnetics Analysis with Higher-Order FEM

AbstractThis paper presents a GPU-accelerated implementation of the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) method with an inexact nullspace filtering approach to find eigenvalues in electromagnetics analysis with higher-order FEM. The performance of the proposed approach is verified using the Kepler (Tesla K40c) graphics accelerator, and is compared to the performance of the implementation based on functions from the Intel MKL on the Intel Xeon (E5-2680 v3, 12 threads) central processing unit (CPU) executed in parallel mode. Compared to the CPU reference implementation based on the Intel MKL functions, the proposed GPU-based LOBPCG method with inexact nullspace filtering allowed us to achieve up to 2.9-fold acceleration.

Effective Implementation of Matrix–Vector Multiplication on Intel's AVX multicore Processor

Matrix–vector multiplication kernel is one of the most important and common computational operations which form the core of varied important application areas such as scientific and engineering applications. Therefore, it is substantial to optimize and accelerate its implementation. This paper proposes an optimized algorithm for single-precision matrix–vector multiplication (SGEMV) on the Intel core i7 processor. An overview of the Intel's advanced vector extension instructions in implementing dense matrix–vector multiplication kernels in parallel has been comprehensively addressed. Also, a variety of performance optimization techniques using Intel's advanced vector extension (AVX) instruction sets, memory access optimization, and OpenMP parallelization has been designed. Additionally, the performance of the proposed algorithms is evaluated in compared to the latest version of Intel Math Kernel Library SGEMV 2017 subroutines because Intel Math Kernel Library subroutines also consider the same optimization methods that are used in this paper. In this paper, we have introduced an overview of the optimization techniques, have explained the specific details of handling them in the proposed algorithm, and also have showed the advantages and the challenges of combining them together in contrast to the previous works which usually have concentrated on a single technique and the performance achieved by it. The guidelines of parallel implementation of the proposed algorithm and the characteristics of the target architecture that should be considered when implementing this algorithm have been investigated. An overview of the Intel's advanced vector extension instructions in implementing dense matrix–vector multiplication kernels in parallel has been comprehensively addressed. A comparative study of the two most popularly used C++ compilers: Intel C++ compiler 17.0 in Intel Parallel Studio XE 2017 against Microsoft Visual Studio C++ compiler 2015 has been investigated. Finally, the comparison between two primary ways of utilizing AVX instructions: inline assembly and intrinsic functions, and the comparison between single-core and multi-core platforms have introduced. The results are evaluated in Intel Core i7-5600U processor of 2.6 GHz with 128 KB L1 cache, 512 KB L2 cache, and 4MB L3 cache running on windows 10 operating system and on a Broadwell system. The obtained results of the proposed optimized algorithm are implemented on square matrices of different large sizes range from 1024 to 19456. The results indicate a performance improvement of 18.2% and 14.1% for (y = A. x) and (y = AT. x) respectively in compared with the results which are obtained using the latest version of Intel Math Kernel Library 2017(SGEMV) subroutines on multi-core platform.

Parallel Solution of Hierarchical Symmetric Positive Definite Linear Systems

Abstract We present a prototype task-parallel algorithm for the solution of hierarchical symmetric positive definite linear systems via the ℋ-Cholesky factorization that builds upon the parallel programming standards and associated runtimes for OpenMP and OmpSs. In contrast with previous efforts, our proposal decouples the numerical aspects of the linear algebra operation from the complexities associated with high performance computing. Our experiments make an exhaustive analysis of the efficiency attained by different parallelization approaches that exploit either task-parallelism or loop-parallelism via a runtime. Alternatively, we also evaluate a solution that leverages multi-threaded parallelism via the parallel implementation of the Basic Linear Algebra Subroutines (BLAS) in Intel MKL.

Cholesky factorization on SIMD multi-core architectures

Abstract Many linear algebra libraries, such as the Intel MKL, Magma or Eigen, provide fast Cholesky factorization. These libraries are suited for big matrices but perform slowly on small ones. Even though State-of-the-Art studies begin to take an interest in small matrices, they usually feature a few hundreds rows. Fields like Computer Vision or High Energy Physics use tiny matrices. In this paper we show that it is possible to speed up the Cholesky factorization for tiny matrices by grouping them in batches and using highly specialized code. We provide High Level Transformations that accelerate the factorization for current multi-core and many-core SIMD architectures ( SSE, AVX 2, KNC, AVX512 , Neon, Altivec). We focus on the fact that, on some architectures, compilers are unable to vectorize and on other architectures, vectorizing compilers are not efficient. Thus hand-made SIMD ization is mandatory. We achieve with these transformations combined with SIMD a speedup from × 14 to × 28 for the whole resolution in single precision compared to the naive code on a AVX 2 machine and a speedup from × 6 to × 14 on double precision, both with a strong scalability.

Evaluation of the Intel Xeon Phi 7120 and NVIDIA K80 as accelerators for two-dimensional panel codes.

To optimize the geometry of airfoils for a specific application is an important engineering problem. In this context genetic algorithms have enjoyed some success as they are able to explore the search space without getting stuck in local optima. However, these algorithms require the computation of aerodynamic properties for a significant number of airfoil geometries. Consequently, for low-speed aerodynamics, panel methods are most often used as the inner solver. In this paper we evaluate the performance of such an optimization algorithm on modern accelerators (more specifically, the Intel Xeon Phi 7120 and the NVIDIA K80). For that purpose, we have implemented an optimized version of the algorithm on the CPU and Xeon Phi (based on OpenMP, vectorization, and the Intel MKL library) and on the GPU (based on CUDA and the MAGMA library). We present timing results for all codes and discuss the similarities and differences between the three implementations. Overall, we observe a speedup of approximately 2.5 for adding an Intel Xeon Phi 7120 to a dual socket workstation and a speedup between 3.4 and 3.8 for adding a NVIDIA K80 to a dual socket workstation.

Basker: Parallel sparse LU factorization utilizing hierarchical parallelism and data layouts

Abstract Transient simulation in circuit simulation tools, such as SPICE and Xyce, depend on scalable and robust sparse LU factorizations for efficient numerical simulation of circuits and power grids. As the need for simulations of very large circuits grow, the prevalence of multicore architectures enable us to use shared memory parallel algorithms for such simulations. A parallel factorization is a critical component of such shared memory parallel simulations. We develop a parallel sparse factorization algorithm that can solve problems from circuit simulations efficiently, and map well to architectural features. This new factorization algorithm exposes hierarchical parallelism to accommodate irregular structure that arise in our target problems. It also uses a hierarchical two-dimensional data layout which reduces synchronization costs and maps to memory hierarchy found in multicore processors. We present an OpenMP based implementation of the parallel algorithm in a new multithreaded solver called Basker in the Trilinos framework. We present performance evaluations of Basker on the Intel SandyBridge and Xeon Phi platforms using circuit and power grid matrices taken from the University of Florida sparse matrix collection and from Xyce circuit simulation. Basker achieves a geometric mean speedup of 5.91× on CPU (16 cores) and 7.4× on Xeon Phi (32 cores) relative to state-of-the-art solver KLU. Basker outperforms Intel MKL Pardiso solver (PMKL) by as much as 30× on CPU (16 cores) and 7.5× on Xeon Phi (32 cores) for low fill-in circuit matrices. Furthermore, Basker provides 5.4× speedup on a challenging matrix sequence taken from an actual Xyce simulation.

Evaluation of MATPOWER and OpenDSS load flow calculations in power systems using parallel computing

This study presents the work carried out by the authors to apply Intel MKL PARDISO (a parallel sparse matrix solver) to the load flow solution algorithms of MATPOWER and OpenDSS. The goal is to explore the potential execution time reduction obtained when working with large power systems and multi-core installations. Test systems of different sizes were solved in order to observe the time reduction as function of the system size and the number of cores used in the parallel execution. Results show that except for the full Newton–Raphson algorithm, one should not expect a reduction in execution times when using a parallel routine. The use of a parallel sparse matrix solver is only justified when the sparse matrix must be recalculated at every step of the solution process or the systems under study are much larger than those analysed here. This study also presents how a parallel computing solution can be implemented in different applications by using available high-performance parallel libraries.

Multi-dimensional Homomorphisms and Their Implementation in OpenCL

Homomorphisms (traditionally defined on lists) are functions that can be parallelized by the divide-and-conquer paradigm. In this paper, we introduce an extension of the traditional homomorphism concept—multi-dimensional homomorphisms (MDHs)—which capture parallelism on multi-dimensional arrays. We propose md_hom—a new parallel pattern (a.k.a. algorithmic skeleton), based on the MDH concept, to simplify parallel programming for a broad class of applications. The md_hom pattern is general enough to subsume common parallel patterns such as map and reduce, and also more complex functions built by composing and nesting several patterns. We present a generic implementation schema for md_hom in form of an efficient, correct-by-construction OpenCL pseudocode that targets various parallel architectures such as multi-core CPU and graphics processing unit (GPU). We develop our pseudocode schema as parametrized in tuning parameters: these allow to optimize the code for different devices and input sizes by performing an automated search on the parameter space. We evaluate the schematically generated, executable OpenCL code using the example of general matrix–vector multiplication (GEMV)—an important linear algebra routine which has gained more attention recently due to its use in the application area of deep learning—on two parallel architectures—Intel CPU and NVIDIA GPU. Our performance results are competitive and in some cases even better than the hand-tuned GEMV implementations provided by the state-of-the-art libraries Intel MKL and NVIDIA cuBLAS, as well as the auto-tunable OpenCL BLAS library CLBlast.

LU factorization on heterogeneous systems: an energy-efficient approach towards high performance

Dense lower–upper (LU) factorization (hereafter referred to as LU) is a critical kernel that is widely used to solve dense linear algebra problems. Hybrid LU algorithms have been well designed to exploit the full capacity of heterogeneous systems. However, existing heterogeneous implementations are typically CPU-centric, which rely highly on CPU cores and suffer from a large amount of data transfers via the PCIe bus, and thus reduce the overall energy efficiency of the entire computer system. In this paper, we provide a coprocessor-resident implementation of LU for a heterogeneous platform to improve energy efficiency by relieving the CPUs from performing heavy load computations and avoiding excessive data transfers via PCIe. To maintain the performance, we conduct optimizations to pipeline the CPU computation, coprocessor computation, MPI communication, and PCIe transfer between the CPUs and coprocessors. The experiments on the Tianhe-2 supercomputer show that our LU implementation can compete with the highly optimized Intel MKL implementation in performance and overcome the limitations of energy efficiency.

A performance analysis of computing the LU and the QR matrix decompositions on the CPU and the GPU

We present an analysis of time efficiency of five different implementations of the LU and the QR decomposition of matrices performed on central processing unit (CPUs) and graphics processing units (GPUs). Three of the considered implementations, developed using the Eigen C++ library, Intel MKL, and MATLAB are executed on a multi-core CPU. The remaining two implementations are processed on a GPU and employ MATLAB's Parallel Computing Toolbox and Nvidia CUDA augmented with the cuSolver library. Computation times are compared using randomly generated single- and double-precision floating-point matrices. The experiments for the LU decomposition show that the two GPU implementations offer best performance for matrices that can fit into the GPU global memory. For larger LU decomposition problem instances, Intel MKL on the CPU is found to be the fastest approach. Furthermore, Intel MKL also proves to be the fastest method for computing QR decomposition for all considered sizes of matrices.

Analysis of a Splitting Approach for the Parallel Solution of Linear Systems on GPU Cards

We discuss an approach for solving sparse or dense banded linear systems $Ax = b$ on a graphics processing unit (GPU) card. The matrix ${A} \in {\mathbb{R}}^{N \times N}$ is possibly nonsymmetric and moderately large, i.e., ${10,000}{} \leq N \leq {500,000}{}$. The split and parallelize (SaP) approach seeks to partition the matrix $A$ into diagonal subblocks $A_i$, $i=1,\ldots,P$, which are independently factored in parallel. The solution may choose to consider or to ignore the matrices that couple the diagonal subblocks $A_i$. This approach, along with the Krylov-subspace-based iterative method that it preconditions, are implemented in a solver called SaP::GPU, which is compared in terms of efficiency with three commonly used sparse direct solvers: PARDISO, SuperLU, and MUMPS. SaP::GPU, which runs entirely on the GPU except for several stages involved in preliminary row and column permutations, is robust and compares well in terms of efficiency with the aforementioned direct solvers. In a comparison against Intel's MKL, SaP::GPU also fared well when used to solve dense banded systems that are close to being diagonally dominant. SaP::GPU is publicly available and distributed as open source under a permissive BSD-3 license.

Task-based Cholesky decomposition on Xeon Phi architectures using OpenMP

The increasing number of computational cores in modern many-core processors, as represented by the Intel Xeon Phi architectures, has created the need for an open-source, high performance and scalable task-based dense linear algebra package that can efficiently use this type of many-core hardware. In this paper, we examined the design modifications necessary when porting PLASMA, a task-based dense linear algebra library, run effectively on two generations of Intel's Xeon Phi architecture, known as knights corner (KNC) and knights landing (KNL). First, we modified PLASMA's tiled Cholesky decomposition to use OpenMP tasks for its scheduling mechanism to enable Xeon Phi compatibility. We then compared the performance of our modified code to that of the original dynamic scheduler running on an Intel Xeon Sandy Bridge CPU. Finally, we looked at the performance of the OpenMP tiled Cholesky decomposition on knights corner and knights landing processors. We detail the optimisations required to obtain performance on these platforms and compare with the highly tuned Intel MKL math library.

Semi-External Memory Sparse Matrix Multiplication for Billion-Node Graphs

Sparse matrix multiplication is traditionally performed in memory and scales to large matrices using the distributed memory of multiple nodes. In contrast, we scale sparse matrix multiplication beyond memory capacity by implementing sparse matrix dense matrix multiplication (SpMM) in a semi-external memory (SEM) fashion; i.e., we keep the sparse matrix on commodity SSDs and dense matrices in memory. Our SEM-SpMM incorporates many in-memory optimizations for large power-law graphs. It outperforms the in-memory implementations of Trilinos and Intel MKL and scales to billion-node graphs, far beyond the limitations of memory. Furthermore, on a single large parallel machine, our SEM-SpMM operates as fast as the distributed implementations of Trilinos using five times as much processing power. We also run our implementation in memory (IM-SpMM) to quantify the overhead of keeping data on SSDs. SEM-SpMM achieves almost 100% performance of IM-SpMM on graphs when the dense matrix has more than four columns; it achieves at least 65% performance of IM-SpMM on all inputs. We apply our SpMM to three important data analysis tasks--PageRank, eigensolving, and non-negative matrix factorization--and show that our SEM implementations significantly advance the state of the art.

Explicit Fourth-Order Runge\u2013Kutta Method on Intel Xeon Phi Coprocessor

This paper concerns an Intel Xeon Phi implementation of the explicit fourth-order Runge–Kutta method (RK4) for very sparse matrices with very short rows. Such matrices arise during Markovian modeling of computer and telecommunication networks. In this work an implementation based on Intel Math Kernel Library (Intel MKL) routines and the authors’ own implementation, both using the CSR storage scheme and working on Intel Xeon Phi, were investigated. The implementation based on the Intel MKL library uses the high-performance BLAS and Sparse BLAS routines. In our application we focus on OpenMP style programming. We implement SpMV operation and vector addition using the basic optimizing techniques and the vectorization. We evaluate our approach in native and offload modes for various number of cores and thread allocation affinities. Both implementations (based on Intel MKL and made by the authors) were compared in respect of the time, the speedup and the performance. The numerical experiments on Intel Xeon Phi show that the performance of authors’ implementation is very promising and gives a gain of up to two times compared to the multithreaded implementation (based on Intel MKL) running on CPU (Intel Xeon processor) and even three times in comparison with the application which uses Intel MKL on Intel Xeon Phi.

OptQC v1.3: An (updated) optimized parallel quantum compiler

We present a revised version of the OptQC program of Loke et al. (2014) [1]. We have removed the simulated annealing process in favour of a descending random walk. We have also introduced a new method for iteratively generating permutation matrices during the random walk process, providing a reduced total cost for implementing the quantum circuit. Lastly, we have also added a synchronization mechanism between threads, giving quicker convergence to more optimal solutions. New version program summaryProgram title: OptQC v1.3Catalogue identifier: AEUA_v1_3Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEUA_v1_3.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 240903No. of bytes in distributed program, including test data, etc.: 632395Distribution format: tar.gzProgramming language: Fortran, MPI.Computer: Any computer with Fortran compiler (not gfortran4.9 or earlier) and MPI library.Operating system: Linux.Classification: 4.15.Catalogue identifier of previous version: AEUA_v1_3Journal reference of previous version: Comput. Phys. Comm. 185(2014)3307External routines: Intel MKL LAPACK routines and MPI routines.Does the new version supersede the previous version?: YesNature of problem:It aims to minimize the number of quantum gates required to implement a given unitary operation.Solution method:It utilizes a descending random walk to select permutation matrices P and Q for a given unitary matrix U such that the number of gates in the quantum circuit of U=QTPTU′PQ is minimized, where U′ is equivalent to U up to a permutation. The decomposition of a unitary operator is performed by recursively applying the cosine–sine decomposition.Reasons for new version:Simulated annealing process was found to give suboptimal results compared to a normal descending random walk. Computation time was also bloated by the necessity of running the CS decomposition thrice (for U′,P and PT) for each iteration of the optimization process.Summary of revisions:•Simulated annealing process was replaced by a descending random walk (equivalent to setting the threshold value β to 0), due to poor convergence of the simulated annealing process, and due to the lack of pathological instances of local minima in this particular search space.•Introduced an iterative method for generating the general permutation matrix P, in which we start with P=I and build up the quantum circuit (and modify P correspondingly) by adding gates onto the existing circuit. This removes the requirement of running the CS decomposition on P and PT to find the quantum circuit implementation, since it is already known by construction.•Added a synchronization mechanism between threads (after some prescribed number of iterations) in which the current state of the top 10% processes with the fittest solutions is copied over to the remaining 90%. This works so as to discard the less fit solutions and focuses the searching algorithm in the state space with the fittest solutions.Additional comments:The program contains some Fortran2003 features and will not compile with gcc4.9 or earlier.Running time:As before, running time increases with the size of the unitary matrix, as well as the prescribed maximum number of iterations for qubit permutation selection and the descending random walk. All simulation results presented in this paper are obtained from running the program on the Fornax supercomputer managed by iVEC@UWA with Intel Xeon X5650 CPUs. A comparison of running times is also given in Table 1.

An improved divide-and-conquer algorithm for the banded matrices with narrow bandwidths

In this paper we propose a novel divide-and-conquer (DC) algorithm to compute the SVD of banded matrices, and further accelerate it by using rank-structured matrix techniques, especially the hierarchically semiseparable (HSS) matrix. The DC algorithm for the symmetric banded eigenvalue problem can also be accelerated similarly. For matrices with few deflations, the banded DC algorithms require more flops than the classical DC algorithm, and thus they are suitable for narrowly banded matrices. While, if there exist many deflations, the banded DC algorithms can be faster than the classical ones for matrices with relatively large bandwidths. Numerous experiments have been done to test the proposed algorithms. Some of the tested matrices are from construction and some are from real applications. Comparing with the DC algorithm in Intel MKL, our proposed algorithms can be hundreds times faster for matrices with narrow bandwidths or many deflations.

New fast divide‐and‐conquer algorithms for the symmetric tridiagonal eigenvalue problem

SummaryIn this paper, two accelerated divide‐and‐conquer (ADC) algorithms are proposed for the symmetric tridiagonal eigenvalue problem, which cost O(N2r) flops in the worst case, where N is the dimension of the matrix and r is a modest number depending on the distribution of eigenvalues. Both of these algorithms use hierarchically semiseparable (HSS) matrices to approximate some intermediate eigenvector matrices, which are Cauchy‐like matrices and are off‐diagonally low‐rank. The difference of these two versions lies in using different HSS construction algorithms, one (denoted by ADC1) uses a structured low‐rank approximation method and the other (ADC2) uses a randomized HSS construction algorithm. For the ADC2 algorithm, a method is proposed to estimate the off‐diagonal rank. Numerous experiments have been carried out to show their stability and efficiency. These algorithms are implemented in parallel in a shared memory environment, and some parallel implementation details are included. Comparing the ADCs with highly optimized multithreaded libraries such as Intel MKL, we find that ADCs could be more than six times faster for some large matrices with few deflations. Copyright © 2016 John Wiley & Sons, Ltd.