Optimizing Communication by Compression for Multi-GPU Scalable Breadth-First Searches

TL;DR

This paper introduces a compression scheme to optimize communication in distributed BFS algorithms on multi-GPU systems, aiming to improve performance by reducing inter-processor communication overhead.

Contribution

It proposes a novel compression method tailored for GPU-based distributed BFS, addressing communication bottlenecks in parallel graph processing.

Findings

Reduced communication volume in multi-GPU BFS

Improved scalability of BFS algorithms on distributed systems

Potential for faster graph analytics in GPU clusters

Abstract

The Breadth First Search (BFS) algorithm is the foundation and building block of many higher graph-based operations such as spanning trees, shortest paths and betweenness centrality. The importance of this algorithm increases each day due to it is a key requirement for many data structures which are becoming popular nowadays. When the BFS algorithm is parallelized by distributing the graph between several processors the interconnection network limits the performance. Hence, improvements on this area may benefit the overall performance of the algorithm. This work presents an alternative compression scheme for communications in distributed BFS processing. It focuses on BFS processors using General-Purpose Graphics Processing Units.

Tables

Table 1. Table 2: Problem sizes for a Graph 500 BFS Benchmark

Problem class	Scale	Edge factor	Approx. storage size in TB
Toy (level 10)	26	16	0.0172
Mini (level 11)	29	16	0.1374
Small (level 12)	32	16	1.0995
Medium (level 13)	36	16	17.5922
Large (level 14)	39	16	140.7375
Huge (level 15)	42	16	1125.8999

Table 2. Table 1: Other implementations

Optimizations	SC11	ISC12	SC12	ISC14
2D decomposition	✓	✓	✓	✓
Vertex sorting	✓
Direction optimization				✓
Data compression	✓	✓	✓
Sparse vector with pop counting				✓
Adaptive data representation				✓
Overlapped communication	✓	✓	✓	✓
Shared memory				✓
GPGPU		✓	✓
NUMA-aware				✓

Table 3. Table 3: Configuration used in both the real sample of a Frontier Queue transfer and the TREC GOV2 Corpus. Both tests have been performed in the “Creek” experiment platform (Table 1 ).

	Extracted Frontier Queue from our BFS
Experiment platform	“Creek”
Prob. Distribution	Uniform (slightly skewed)
Empirical Entropy	14.58845 bit
Skewness	-0.003576766
Integer Range	[1-65532]
Total integers in sample	29899
	TREC GOV2 Corpus
Experiment platform	“Creek”
Prob. Distribution	Power Law
Total integers in sample	23918M

Equations

$$P=C\ast V^2\ast F$$

$$A=\left[\begin{array}[]{cc}A_{1}\\ \hline\cr A_{2}\\ \hline\cr\vdots\\ \hline\cr A_{k}\end{array}\right]$$

$$A=\left[\begin{array}[]{ccc}A_{1,1}^{(1)}&\cdots&A_{1,C}^{(1)}\\ \hline\cr\vdots&\ddots&\vdots\\ \hline\cr A_{R,1}^{(1)}&\cdots&A_{R,C}^{(1)}\\ \hline\cr\vdots&\ddots&\vdots\\ \vdots&\ddots&\vdots\\ A_{1,1}^{(C)}&\cdots&A_{1,C}^{(C)}\\ \hline\cr\vdots&\ddots&\vdots\\ \hline\cr A_{R,1}^{(C)}&\cdots&A_{R,C}^{(C)}\end{array}\right]$$

$$vertices=2^{scale}\wedge edges=2^{scale}\times edgefactor$$

$$H(X)=-\sum _{i=0}^{N-1}{p}_i{\log }_2{p}_i$$

$$bits=\big{[}H(X)\big{]}$$

$$ratio$$

$$\exists min\in B\cdot \forall i\in B\cdot (i^{‘}=i-min)$$

$$size=\big{[}\log_{2}{(max+1-min)}\big{]}$$

$$\big{[}\log_{2}(2127+1-2000)\big{]}=7$$

$$5,120,300,530,990,25,980,799273656454$$

$$\big{[}\log_{2}(990+1-5)\big{]}=9,9$$

$$ratio\%=\frac{1}{compression\_ratio}\ast 100$$

Full text

1. 0 Introduction

2. 1 Introduction

3. 2 Motivation

4. 3 Contributions

5. 4 Structure of the thesis

6. 1 Background

7. 1 Milestones in supercomputing

8. 2 Global architectural concepts

9. 3 Clusters and High Performance Computing

10. 1 Hardware architectures

11. 2 Interconnection networks in HPC

12. 3 The classification of supercomputers used in HPC

13. 4 The popularization of the multicore architectures

14. 5 Energy efficiency

15. 4 General Purpose Graphic Proccessor Units

16. 1 GPGPU architectures

17. 2 Regular and non-regular problems

18. 5 Message Passing Interface

19. 6 Graph computations

20. 1 Graphs

21. 2 Graphs partitioning

22. 3 Graph 500 challenge

23. 7 Input data

24. 1 Synthetic data

25. 2 Real world graphs (datasets)

26. 2 Related work

27. 1 Optimizations

28. 2 Other implementations

29. 3 Compression

30. 3 Problem Analysis

31. 1 Our initial implementation

32. 1 Optimizations table in “Baseline” implementation

33. 2 Our “Baseline” algorithm

34. 3 General communication algorithm

35. 4 Data and communications

36. 2 Instrumentation

37. 1 Instrumented zones

38. 3 Analysis

39. 1 Communication overhead

40. 2 Instruction overhead

41. 4 Optimizing data movement by compression

42. 1 Concepts

43. 2 Compression algorithms (codecs or schemes)

44. 3 Compression libraries

45. 1 Lemire et al.

46. 2 Turbo-PFOR

47. 3 Alenka CUDA

48. 4 Integration of the compression

49. 1 Experimental performance comparison

50. 2 Observations about experimental performance

51. 3 Usage of compression thresholds

52. 4 Communication before and after compression

53. 5 Optimizing instruction overhead

54. 1 Scalar optimizations

55. 1 Strength reduction

56. 2 Vectorization

57. 1 Compiler optimizations

58. 3 Thread parallelism

59. 1 Usage of fork-join model when suitable

60. 2 Optimal thread scheduling

61. 3 Thread contention prevention

62. 4 NUMA control / NUMA aware

63. 4 Memory access

64. 1 Locality of data access

65. 2 Merge subsequent loops

66. 3 Variable caching in the outer loop

67. 5 Communication patterns

68. 1 Data transfer grouping

69. 2 Reduce communication overhead by overlapping

70. 6 Summary of implemented optimizations

71. 6 Final Results

72. 1 Experiment platforms

73. 2 Results

74. 1 Scalability analysis

75. 2 Overhead of the added compression

76. 3 Instruction overhead analysis

77. 4 Compression analysis

78. 7 Conclusions

79. 8 Future work

Chapter 0 Introduction

1 Introduction

The study of large Graphs datasets has been a very popular subject for some years. Fields such as Electronic engineering, Bioscience, the World Wide Web, Social networks, Data mining, and the recent expansion of “Big Data” result in larger amounts of data each time, which need to be explored by the current and state-of-the-art graph algorithms. One result of this recent expansion of large-sized graph is the establishment of benchmarks that engage organisations in the development of new techniques and algorithms that solve this new challenge.

Graph 500111http://www.graph500.org is a new benchmark created in 2010 which serves this purpose. It uses the Breadth-first Search (BFS) algorithm as a foundation. This benchmark contrasts with other well-known and more mature tests like LINPACK (used in the Top 500222http://www.top500.org/ Challenge), in that the latter executes linear equations in a computation-intensive fashion. In contrast, the former makes use of graph based data-intensive operations to traverse the graph.

In this work we focus on possible optimizations of the Breadth-first Search (BFS) algorithm.

2 Motivation

As stated before, as the growth of data volumes increases at a very high rate and the structures containing this data use an internal graph representation, the research in graph algorithms also increases.

Some of the current research done on the area of the BFS algorithms is framed under the Graph 500 challenge.

The work done on improving our previous implementation of the garph500 benchmarch (‘‘Baseline’’), as well as the research of other authors in this area, shows that there is a severe impact in the overall performance of the distributed version of the algorithm due to the latency introduced by the the data movement between processors.

3 Contributions

We believe that the conclusions presented in this work may benefit a broad variety of Breadth-first Search implementations. Compressing the data movement may alleviate the main bottleneck of the distributed version of the algorithm: the communication. In this work we make use of a Compressed Sparse Row (CSR) matrix, 2D data partitioning and Sparse vector multiplications (SpMVM). By the analysis of the transmitted data (a SpMV vector represented as an integer sequence) we have been able to successfully achieve over 90% in terms of data transfer reduction and over an 80% of communication time reduction. All this is achieved thanks to specific numerical proprieties the the data of the analyzed graph (Section 6).

Furthermore, we will also show how different techniques to reduce the instruction overhead may improve the performance, in terms of Traversed Edges per Second TEPS (Section 3) and time, of a graph500 implementation. Finally, we will discuss further optimizations. Here follows the main contributions of this work:

As mayor goal, we have used a Bit Packing compression scheme (Section 2 with delta compression on top, optimized for the SIMD instruction set of Intel™ x86 architecture. The reasoning and criteria of the used algorithm are discussed in further detail on Section 4. Also, a scalability study before and after integrating the compression has been performed in sections 4 and 6. 2. 2.

In addition to the added compression library (which contains several compression codecs), two more libraries have been integrated (partially). With this we enable our graph500 application to test more codecs for other different graphs and problem sizes. This is described in further detail in section 3. The newly added implementations are not fully operative yet. These new packages offer benefits such as new compression codecs, different input vector sizes (which would allow bigger problem sizes and the avoidance of the performance penalty imposed by the conversion of the full vector, required by the current library), and a GPU-based implementation. This converts our graph500 implementation on a modular compression test-bench. 3. 3.

As part of this work, the instruction overhead has been reduced. With this, we have tested the effects of some common instruction transformation techniques at different scales in our BFS implementation. The used techniques have been described more in detail on Section 5 (Optimizing instruction overhead) 4. 4.

An external instrumentation (Score-P) has been added to help with the detection of the bottlenecks, optimization of our code and validation of the final results. The reasoning for the profiler selection criteria is discussed in Section 2. As a summary, the selected profiler provides a fine-grain way to gather a high level of details with low overhead 333http://www.vi-hps.org/projects/score-p/. 5. 5.

Some other tasks have been performed to our Graph 500 implementation. Briefly, these have been:

• –

The application’s built process has been improved to reduce the building time. Also, a pre-processing of the compiled code allows a more fine-grain tuning. For this purpose we have chosen Maketools 444https://www.gnu.org/software/make/. As result, we detect the maximum capabilities of the target system at pre-compile time.

• –

We have created scripts to automate the visualization and interpretation of the results. For this we combine c/c++ with R.

4 Structure of the thesis

The thesis is formed by 4 main parts (i) a background section about the involved elements of the thesis. (ii) An analysis of our problem with two possible solutions: compression and instruction overhead reduction. (iii) Finally in the last part, we discuss the result and see how much do they mach with the purposed solutions for our problem.

Chapter 1 Background

As the background of optimizations made to a Graph 500 application, a short timeline of the history of Supercomputing will be listed to place the cited concepts on the time. Some of the core optimizations technologies, used in both the compression and the main Graph 500 implementation will be reviewed.

Regarding Supercomputing, it will be provided a background on some of its core concepts: High Performance Computing (HPC), General Processing Graphic Processing Units (GPGPU) and Message Passing Interface (MPI).

*Last, concepts about graphs (in a Linear Algebra context), and the datasets used often on graphs literature will also be described. *

1 Milestones in supercomputing

As an introduction to Supercomputing, some of the personal names and milestones of its history111http://www.computerhistory.org are listed below. Many of these will be referenced along in this work.

• •

1930’s - Fist theoretical basis of computing.

• –

1936: Alan Turing develops the notion of a “Universal machine” through his paper “On Computable Numbers” [50]. It is capable of computing anything that is computable.

• •

1940’s - The foundation of computing.

• –

1945: John von Neumann introduces the concept of a stored program in the draft report on the EDVAC design. This way he creates the design of a sequential computer [53] (known as the von Neumann architecture).

• –

1947: Bell Laboratories invents the Transistor.

• –

1948: Claude Shannon writes “A Mathematical Theory of Communication” (Sec. 1 ) [45]

• •

1950’s - The first commercial general-purpose computers are born.

• –

1951: UNIVAC-1 is launched. This is the very first commercial general-purpose computer. Even though it is general-purpose, it is focus on home use.

• –

1953: IBM announces Model 650. Stores the data on rotating tape. This is the first mass-produced computer.

• –

1953: IBM announces Model 701. The first Mainframe222These computers, also referred as “Big Iron”, receive their name for the shape of their cabinet structure. They are often used in statistical computation, banking transactions, etc.

• –

1956: UNIVAC announces an UNIVAC built with transistors.

• •

1960’s - During this period the IBM Corporation dominated the general purpose early computing industry. This was the era of the Mainframes. There are new additions to the foundation of computing.

• –

1965: Gordon Moore coins the Moore’s Law. This predicts the number of in-chip transistors throughout time.

• –

1966: Michael Flynn proposes a classification of computer architectures. This is known as Flynn’s Taxonomy. It is still used nowadays.

• –

1968: CDC 7600 is launched by Seymour Cray. It is considered by many the first true Supercomputer.

• •

1970’s - First Supercomputers enter in scene. They increase their performance through Vector Processors Units among other methods.

• –

1971: Intel Corporation launches its first chip: Intel 4004

• –

1972: Cray Research Inc. is created by Seymour Cray. He is consider as the father of Supercomputing.

• –

1975: Cray Research Inc completes the development of Cray 1 - The first Vector-processor Supercomputer.

• –

1978: Intel Corporation introduces the first 16-bit processor, the 8086.

• •

1980’s - Many technological improvements. First IBM/Intel personal computers (PC). First Connection Machine Supercomputers.

• –

1982: Cray Research Inc introduces the Cray X-MP Supercomputer. This version uses Shared Memory and a vector-processor. It is a ’cleaned-up’ version of the CRAY-1. Successive versions of this Supercomputer increase the number of CPUs, raise the clock frequency and expand the instruction size from 24 bits to 32.

• –

1983: Thinking Machines introduces CM-1, the first Connection Machine (CM). It is based upon the SIMD classification.

• –

1984: IBM introduces the IBM PC/AT based on the chip Intel 80286. The chip works at 16 bits.

• –

1985: Intel Corporation introduces the 80386 chip with 32-bit processing and on-chip memory management.

• •

1990’s - Next step in the Supercomputer era. First Massively Parallel architectures. The home computer industry is based on ’clones’ of the original IBM PC.

• –

1991: Thinking Machines introduces CM-5, a new version of their Connection Machine running on an RISC SPARC, and replacing the connection network of the previous CM-2. The CM-5 changes, this way, to a MIMD design. This is the first NUMA architecture.

• –

1992: Touchstone Delta Is an experiment carried out in Caltech with 64x Intel 8086 microprocessors. It opens a door to a new era of parallelism. This is later referred as “The attack of the killer micros”. Here, a large number of Commercial off-the-shelf (COTS) microprocessors, invaded a world dominated by “Big Iron” Vector computers.

• –

1993: Intel Corporation releases “Pentium” chip. Personal computers continue to grow.

• –

1993: top500.org ranking is created. It uses a linear algebra LINPACK benchmark.

• –

1997: The Intel ASCI Red Supercomputer was developed based on the Touchstone Delta experiment. This Massively parallel Supercomputer was the fastest in the world until the early 2010’s.

• –

1999: IBM PowerPC 440 microprocessor is launched. This 32-bit RISC high performance core will be the main processor of many future HPP architectures like the Blue Gene/L or Cray XT3

• •

2000’s - The fastest Supercomputers are Massively Parallel architectures.

• –

2001: General Purpose GPU processing (GPGPU) begins its development by the advent of the programmable shaders and floating point units on graphics processors.

• –

2004: Intel Corporation announces the cancelation of two of their processors [36]. This is known as the End of Frequency Scaling. The new speedups are based upon parallel techniques developed in the previous Supercomputing era. With the Multicore the personal computer industry also go parallel.

• •

2010’s - The only Supercomputer architectures in Top500333http://www.top500.org/statistics/overtime/ are Massively Parallel: MPPs and Clusters.

• –

2010: graph500.org challenge is created. It uses graph-based data insensitive computation over a Breadth-first Search (BFS) algorithm.

• –

2012: green.graph500.org challenge is created. The current energy awareness leads to this energy focused variation of the Graph 500 Challenge. This new benchmark is also based upon the BFS algorithm but reflexes the power efficiency in the results.

2 Global architectural concepts

Continuing with the background about High Performance Computing (HPC) and related technologies like GPGPUs, we will introduce in this section some concepts which will be referenced in the document. This concepts encompass the internal methods used in the compression algorithm (SIMD / Vectorization [28]), other optimizations made to the selected compression algorithm (Super-Scalar optimizations, improvements on the Pipelinig [68]). These concept will also be referenced in this Background section for a better understanding of HPC, GPGPUs and Supercomputing in general.

The Von Neumann design

It was in 1945 when the physicist and mathematician John von Neumann and others, designed the first programmable architecture. It was called Von Neumann architecture and had an intrinsically serial design. In this, there is only one processor executing a series of instructions and the execution flow occur in the same order as it appears in the original program [53].

Moore’s forecast

Gordon Moore, in 1965, created a forecast for the number of in-chip transistors on a silicon microprocessor [35]. It was able to predicted the number of transistors on a chip troughtout the time: “In-chip transistror number would double every 2 years” (Figure 1).

This prediction was true for a long time and the clock frequency of microprocessors was increasing due to deeper Pipelines and more transistors. In the early 2000’s the micro processors reached a temperature limit [35]. Because of this, clock frequency has not been since 2005 (Sec. 5).

This, has ultimately led to Multi-core technology (a chip with multiple cores) allowing a similar speed. As a result, nowadays software applications are required to have parallelism to benefit from this [34].

Pipeline instruction parallelism

The Instruction Pipelining technique uses instruction-level parallelism inside a single processor. The basic instruction cycle is chopped into a series which is called a Pipeline. This allows a quicker execution throughput. Put simply, instead of processing each instruction sequentially (completing the instruction before the next one), each instruction is split into sub steps. These sub steps can be executed in parallel.

Superscalar architectures

A Superscalar processor is an architecture in which the processor is able to issue multiple instructions in a single clock. This is achieved using redundant facilities to execute an instruction.

Each of the superscalar replicated execution units is a resource within a single CPU, such units may be an Arithmetic Logic Unit, a Vector Processing Unit, bitwise operators, or sets of CPU registers (Intel™ Multithreading). This excludes the case of the replicated unit being a separate processing unit (for example other core in a Multi-core processor).

Note that, while a superscalar design is typically also pipelined, Pipelining and Superscalar architectures are different concepts.

Flynn’s taxonomy

In 1996, Michael Flynn creates the earliest classification systems for parallel and sequential computers [20]. The categories in the taxonomy may be seen in Table 1

He classified systems no matter if they were functioning using a single set or multiple sets of instructions. Also to whether (or not) the instructions were using a single set or multiple sets of data. Flynn’s Taxonomy is still in use today.

According to the authors Patterson and Hennessy, “Some machines are hybrids of these categories, but this classic model has survived because it is simple, easy to understand, and gives a good approximation” [26].

Two terms in Flynn’s Taxonomy relevant to this work: SIMD and MIMD are described below.

Simple Instruction Multiple Data (SIMD)

In a SIMD architecture, the parallelism is in the data. There is only one program counter, and it moves through a set of instructions. One of these instructions may operate on multiple data elements in parallel, at a time.

Examples in this category are a modern GPGPU, or Intel SSE and AltiVec Instruction-sets. 2. 2.

Multiple Instruction Multiple Data (MIMD)

MIMD are autonomous processors simultaneously executing different instructions on different data.

Distributed systems are generally categorized as MIMD architectures. In these systems each processor has an individual memory. Also, each processor has no knowledge about the memory in other processors. In order to share data, this must be passed as a message from one processor to another.

Examples of this category are a Multicore-superscalar processor, a modern MPP, or a Cluster.

As of 2015, all the Top500 Supercomputers are within this category.

Types of parallelism according to their level of synchronization

Other way to classify the parallelism is by the concept of Granularity. In this context we relate the amount of computation with the amount of communication due to synchronization. This classification would be as follows:

• •

fine-grained parallelism In this class of granularity the tasks have small computation and high amount of synchronization. As a downside of this level of synchronization the overhead due to communication is bigger.

• •

coarse-grained parallelism In this type of granularity, the tasks are bigger, with less communication between them. As a result the overhead due to communication is smaller.

In order to attain a good parallel performance, the best balance between load and communication overhead needs to be found. If the granularity is too fine, the application will suffer from the overhead due to communication. On the other side, if the granularity is too coarse, the performance can suffer from load imbalance.

3 Clusters and High Performance Computing

High Performance Computing (HPC), or supercomputing, is a kind of computation focused in high consuming computational operations where specific parallel paradigms and techniques are required. Some examples of these are predictive models, scientific computation, or simulations. The types of Supercomputers (and their features), and how they evolved in history until today, are defined in this section.

1 Hardware architectures

Based upon the taxonomy purposed by Flynn we show a general classification of the HPC architectures. This can be see in Figure 2.

2 Interconnection networks in HPC

As this topic is very broad we will only review a concept related to networking, GPGPU, hardware architectures and in an extend, to this work.

The Remote Direct Memory Access (RDMA) is a technique created to solve a problem in networking applications. First, the problem is the performance penalty imposed when a system (sometimes, a processor unit within the main system) accesses the memory of another system (sometimes the original one) involving the Operating system (e.g. CPU, caches, memory, context switches) of the first one. The solution for this problem consists of the access via the network device to the peer system.

Some implementations of this technique use hardware specific vendors as Infiniband™ , and others give support to already existing technologies as Ethernet, such as RoCE and iWarp.

The specific BFS implementation based on GPGPU, in which multiple GPUS might be allocated within the same system, would suffer a high penalty without this technique. In our tests, the RDMA support is enabled and managed by the Message Passing Interface MPI implementation (Section 5)

3 The classification of supercomputers used in HPC

The first computers in the market were mainly built for governmental, scientific or military purposes. It could be stated that computing started with Supercomputing. Here we list a classification of supercomputers according to similar design characteristics and sorted chronologically.

Vector supercomputers

These were the first to appear [62]. In these architectures, the parallelism was achieved with the help of Vectorizing Compilers. The applications needed to be adapted for their use with these compilers.

In this scenario some loops could sometimes be vectorized through annotations in the language. Other loops needed to be arranged in order to remove the dependencies.

The main benefits of these Supercomputers were for programs that would involve a lot of array processing.

They enter in the Flynn’s Taxonomy (Table 1 page 1) under the SIMDs category and unlike other SIMD architectures like the (described in the next section) Connection Machines, the modern GPGPUs or the Intel SSE / AltiVec instruction sets these used larger and variable vector sizes.

Data-parallel architectures

This was a brief moment in time. In this category enter the initial Connection Machines. These architectures also enter in the Flynn’s SIMD category.

Here, the data-parallel language will distribute the work. The algorithms were very clever and involved data-parallel operations. Examples of those operations were segmented scans, Sparse Matrix operations or reductions.

Some parallel languages were C*, CMF (Connection Machine Fortran) or *Lisp [62].

Shared Memory

In these systems, several processors share a common memory space. The programming model in these systems is based on Threads and Locks.

The algorithms were some variation of Parallel Random Access Machine (PRAM) algorithm, or hardware-based algorithms like cache coherence implementations. This latter were supported by hundreds of technical papers. In this classification enter the sub-taxonomy listed in Figure 2. Two of this types will be described below.

in a UMA model all the processors share the physical memory uniformly and access time to a memory location is independent of which processor makes the request. 2. 2.

in a NUMA model the memory access time depends on the memory location relative to the processor: the local memory of a processor will be accessed faster than its remote memory (memory local to another processor or memory shared between other processors). This type of architectures are common in cluster programming and have relevance in this work due to the new advances made over the distributed BFS algorithm for big cluster supercomputers. A diagram of the local and remote banks of memory may be seen in figure 3. This architecture will be referenced in section 2.

Multicomputers

The Multicomputers replaced the Vector Supercomputers. The Shared memory machines were still popular but in contrast with them, these new architectures provided thousands of processes rather than just tens. These two models are not incompatible and are often combined. The programming model was made using Message Passing. In its early days each organization used its own message-passing implementation. Later, the scientific community created an standard, Message Passing Interface or simply MPI [62]. From an algorithm point of view, there had to be made complete changes to the code. This was called domain decomposition and required a complete re-structuring of the program, which needed to be split.

These architecture designs enter in Flynn’s MIMD category. In this group there are two main classes: Massively Parallel Processor (MPP) and Commodity Clusters (or simple Clusters). Both are described below.

MPP Massively Parallel Processor (MPP) are complex and expensive systems. An MPP is a single computer system which sometimes is formed by many nodes. The nodes contain many connected processors which are tightly coupled by specialized interconnect networks (e.g: Hypercube, Mesh or Torus [1]) [26, 17]. Examples of MPPs are the IBM Blue Gene series or the Cray Supercomputers XT30 and XC40 among others. As it will be difined in the section Graph computations (Section 6), these systems are referred as Lightweight by the Graph 500 BoF444http://www.graph500.org/bof. 2. 2.

Clusters systems are usually cited as: “a home-made version of an MPP for a fraction of the price” [47]. However this is not always true, as it can be the example of the Tianhe-2, in China. These, are based upon independent systems (Usually Shared Memory Systems) connected through high speed networks (e.g. Infiniband, Myrinet, Gigabit Ethernet, etc).

As of November 2015, the Top10 architectures in the Top500 555http://top500.org/lists/2015/11/ list are 2 Clusters and 8 MPPs (being the fastest, the Tianhe-2 Cluster in China)

In the Top10 of the Graph500 666http://www.graph500.org/results_nov_2015 list, as of November 2015, there is 1 Cluster and 9 MPPs (being the fastest, the K-computer MPP in Japan)

4 The popularization of the multicore architectures

Due to the change in the manufacturing process of the silicon-chips, the rest of the world followed suit.

The required power to enhance the performance (by increasing the area of the circuit block) varies at a higher rate than the performance does. (Figure 4) [62].

One direct effect of the latter, is that the generated heat increases (proportional to the consumed power). As a result, alleviating the effects of this heat, and keeping the silicon chips within an operative range is difficult. The physical explanation of this is:

The generated heat in a chip is proportional to its Power consumption $P$ [41], which is given by (4). On this, when the frequency increases also the power consumption and the heat rise.

$\textbf{C}\longrightarrow$ The Capacitance being switched per clock cycle (proportional to the number of transistors)

$\textbf{V}\longrightarrow$ Voltage

$\textbf{F}\longrightarrow$ The processor frequency

[TABLE]

This is what ultimately led to Intel’s cancellation of several microprocessors in 2004 [36]. Also, this moment in time is generally cited as the “End of Frequency Scaling”. As this stopped the increasement of the frequency, many cores started to be added into one single chip. There are also more reasons (economical) which affected the evolutions in this chip technology. For example, the relationship between the profits and the technology costs (smaller technologies are more expensive). This ultimately led to an slow down to new smaller technologies [62] (Fig. 5).

From the point of view of the software, this evolution made threads a popular programming model in the last years (as it had been learned a lot about of them with Shared Memory architectures, used in Supercomputing). Also, other programming models like the less restricted “fork-join” threads used in CILK, started to be more popular and practical.

5 Energy efficiency

After going through the limitations that heat imposes to chips, and how the energy consumption rises with the in-chip number of transistors, it can be better understood why energy efficiency is becoming today one of the main design constrains in High Performance Computing.

To illustrate the importance of this, a real example (using modern Supercomputers) will be is used.

As of 2015, in the United States of America, the cost per megawatt is around 1 million dollars [62] and the consumed energy usually follows a chip scaling. With this in mind, in the scenario that we are going to describe, in 2008 a 1 petaflop system used 3 megawatts to run. Following this scaling, in 2018 a 1 exaflop Supercomputer will require 20 megawatts.

To contrast the importance of the energy saving plans, a 20 megawatts supercomputer would require an approximate budget of 20 Million dollars which was a quarter of the yearly budget destined to nursing in that country in 2015. [62]).

The current energy cost in super computers is serious and not sustainable. In the last years researching on energy efficiency has moved also to algorithms and software design. As it will be illustrated in Section 6, new benchmarks (like the green.graph500777http://green.graph500.org/) are being created to measure this Energy impact of computation.

4 General Purpose Graphic Proccessor Units

*Overview *

1 GPGPU architectures

GPUS are widely used as commodity components in modern-day machines. In the case of NVIDIA™ architectures, a GPU consists of many individual execution units (SMs/ SMXs), each of which executes in parallel with the others (Figure 7). During runtime, threads on each execution unit are organized into thread blocks, and each block consists of multiple 32-thread groups, called Warp. If threads within a Warp are set to execute different instructions, they are called “diverged”. Computations in diverged threads are only partially parallel, thus reducing the overall performance significantly.

The GPU includes a large amount of device memory with high bandwidth and high access latency, called global memory. In addition, there is a small amount of shared memory on each execution unit, which is essentially a low latency, high bandwidth memory running at register speeds. Due to such massive amounts of parallelism, GPUs have been adopted to accelerate data and graph processing. However, this maximum potential parallelism is many times difficult to achieve.

2 Regular and non-regular problems

• •

In regular code, control flow and memory references are not data dependent, for example SpMVM is a good example. Knowing only the source code, the input size, and the starting addresses of the matrix and vectors we can predict the program behavior on an in-order processor.

• •

In irregular code, both control flow and memory addresses may be data dependent. The input values determine the program’s runtime behavior, which therefore cannot be statically predicted. For example, in a binary-search-tree implementation, the values and the order in which they are processed affect the control flow and memory references. Processing the values in sorted order will generate a tree with only right children whereas the reverse order will generate a tree with only left children.

Graph-based applications in particular tend to be irregular. Their memory-access patterns are generally data dependent because the connectivity of the graph and the values on nodes and edges determine which graph elements are computed next accessed, but values are unknown before the input graph is available and may change dynamically.

5 Message Passing Interface

The Message Parsing Interface (MPI) resulted as a standardisation of the multiple implementations of communication protocols in early supercomputing.

MPI only defines an interface of calls, datatypes or restriction that are separately implemented by other organizations. Two popular implementations of MPI are OpenMPI and MPICH.

On MPI exist two types of communication calls:

• •

Point-to-Point calls. Allow the direct data send and receive from / to other node. Some examples are MPI_send, MPI_isend (non-blocking) or MPI_Recv.

• •

Collective communications allow us to send data from / to multiple nodes at a time. Some examples are MPI_Reduce, MPI_Allgatherv, MPI_scatter.

6 Graph computations

As examined in the Introduction section (Section 1) Graphs and Graphs computations are becoming more important each day.

In a first part of this section it will be reviewed the Graphs structures and their characteristics.

On a second part of the section it will be discussed the Graph 500 challenge. About the latter it will be described what is it main motivation, what parts conform a Graph 500 application, what is the basic algorithm in a non parallelized application. Lastly it will be discussed why the parallelization of this algorithm (BFS) is theoretically and computationally difficult.

Also other recent benchmark (Green Graph 500), which is very related with this one, will be briefly described.

1 Graphs

The graph algorithms encompass two classes of algorithms: traversal algorithms and analytical iterative algorithms. A description of both follows:

• •

Traversal algorithms involve iterating through vertices of the graph in a graph dependent ordering. Vertices can be traversed one time or multiple times. This class includes search algorithms (such as breadth-first search or depth-first search), single source shortest paths, minimum spanning tree algorithm, connectivity algorithms and so on and so forth.

• •

Analytically iterative algorithms involve iterating over the entire graph multiple times until a convergence condition is reached. This class of algorithms can be efficiently implemented using the Bulk Synchronous Parallel (BSP) model. Algorithms in this class include page rank, community detection, triangle counting and so on.

We will focus in a parallel implementation of the traversal algorithm Breadth-first Search (BFS). This algorithms requires some form of data partition, concurrency control, and thread-level optimizations for optimal performance.

• •

Degree of a graph is the number of edges connected to a vertex.

• •

Distance between to vertices on a graph is the value of the shortest path between to vertices.

2 Graphs partitioning

Parallel SpMV

1D Partitioning-based BFS (Vertex)

In 1D data partitions each processor gets a row of vertices of the adjacency matrix. This kind of partitioning is the default one and has network communication complexity of $\theta(n\times m)$, where n is the number of rows, and m is the number of columns.

[TABLE]

2D Partitioning-based BFS (Edge)

The 2D-partition based BFS is known to perform better than the 1D in graph with low degree. Unfortunately, for high degree graphs the situation becomes the opposite. First of all, as a pre-condition for doing this kind of partition that adjacency matrix must be symmetric in relation to the number of processors, i.e with 36 processors we need to be able to divide the matrix into $6\times 6$ blocks.

[TABLE]

3 Graph 500 challenge

Some concept attaining the Graph 500 challenge are described below.

Performance

the Traverse Edges per Second (TEPS) is used as the main performance metric of the Graph 500 applications. TEPS is proposed by the graph500 benchmark to report the throughput of supercomputers on graph processing.

Energy efficiency

is also a metric which measures the efficiency in terms of energy expressed as traversed edges per joule. Since platforms have heterogeneous hardware resources and different power consumptions, their performance (TEPS) is divided by their power to calculate the energy

Provided reference implementations

The Graph 500 provides five reference implementations with different characteristics. This implementations have the following characteristics

Sequential This reference implementation uses a sequential BFS algorithm and introduces no parallelism but the data vectorizations introduced by the compiler. It could fit for example a UMA architecture, using one only processor on its CPU, such as a mobile phone or a tablet. 2. 2.

OpenMP The benefit of this implementation comes from the hand of the introduced thread parallelism which may be achieved with OpenMP and a multicore architecture. An example architecture which would fit in this category would be a laptop or a Personal computer. 3. 3.

Cray XMT This implementation is specific for Cray™ MPPs. Takes advantage on the specific interconnectors of this HPC supercomputers. 4. 4.

MPI On this implementation the parallelism is based on Flynn’s MIMD category (Table 1 and Figure 2). This is the problem is partitioned and divided into several processors through Message Passing Interface (Section 5) 5. 5.

2D data partitioned MPI.This implementation has been added with posteriority and is based on the 2D problem partitioning purposed in [12]. As the previous one distribute the load among processors using MPI. This reference implementation is the more structurally similar to our one. .

All these implementations are based on CPU processors.

Problem sizes

According to the Graph 500 consortium problem sizes enter into different categories depending on their size. The size is expressed in logarithmic scale and represents the number of total vertices inside the graph. An overview of the categories may be seen in table Table 2.

As it will be described in section 7 and it is cited briefly here, The graph 500 challenge makes use of a Kronecker generator [44] to create the graph. The generated graph (as required in this challenge) contains 16 edges per vertex (Edge factor) . This latter is the reason why the Graph 500 applications works over sparse graphs.

In the table 2 Real world graphs whit a higher number of edges per vertex (degree of the graph) have an scale greater than 30. As it was previously described in this same section a node with logarithmic scale of 30 would have $2^{30}$ vertices and $2^{34}$ Edges.

One note which will also be referenced next in the document is that the meaning of using an Sparse vector with elements of (SpMV) of 64-bit means that a Row (or column for symmetric 2d partitions) of the CSR matrix representing the whole graph, also has 64 elements. This would mean that the maximum represented graph size could have $2^{64}$ Vertices. If we think about this in terms of compression. Is our compressor algorithm works over 32-bit integers (integer is other way of calling the bitmap SpMV) would mean that the maximum archivable graph size (in Table 2) would be small size.

Structure and Kernels

A graph 500 application is formed by 4 different steps, some of which are timed and some of which not. Also, has two differentiated core parts called kernels, in the first one improvements over the structure containing the graph are allowed and timed. The second kernel contains the code of the BFS algorithm. This latter step is repeated 64 times and on each of them the results are verified.

As a result of a benchmark are required some statistical parameters, such as the mean execution time or the Traverse Edges per Second metric (TEPS) calculated as an harmonic mean of all the 64 required iterations. The algorithm of the benchmark may be seen as pseudocode in table 1.

The Green Graph 500 Challenge

The Green Graph 500888http://green.graph500.org/ challenge is an Energy-Aware Benchmark launched in the year 2012. It is based on the similar technical principles than the Graph 500. Some aspects like the metrics (TEPS/Watt), the kernels structure, or the submission rules [27] differ from the previous benchmark.

As it was discussed in the previous section Looking for energy efficient technologies (Section 5), energy is becoming a concern between the HPC community due to the increasingly maintenance costs.

Due to this fact new improved algorithms and implementations are starting to be focused also on a energy-efficient implementation.

7 Input data

1 Synthetic data

In the graph 500 challenge’s algorithm of table 1, the step on line 1 generates a list of edges, that in next steps will be converted into other structures.

The graph 500 consortium sets some specifications for the generated graph and for the used generator. The generator [44] builds open scale graphs which have similar statistic distributions to the graphs observed in the real world. the specifications on the Graph 500 challenge require, apart from the usage of a Kronecker generator, a number of 16 edges per vertex.

This graphs generated synthetically this way have been target of many statistical research. In the figure 8 it may be seen a typical graph 500 graph, following a Power law distribution. One reason for showing how a graph looks like statistically, is to see the comparison with that data that it is distributed by the processors and will be compressed later.

Kroneker generated graphs meet the following relations:

[TABLE]

2 Real world graphs (datasets)

A real world graph usually is greater than a scale 30 kronecker generated graph. Other differences with kronecker are the number of edges per vertex. real world graphs have a much larger edge-factor, and thus a higher diameter (maximum distance between two vertex in the graph).

Chapter 2 Related work

1 Optimizations

2D decomposition This optimization is described in the section Graph computation. 2. 2.

Vertex sorting Vertex sorting is a technique that helps to alleviate one of the major problems of the graph computations. Big size graphs are represented with structures with a very low spatial locality. This means that for example to traverse a graph with a distributed algorithm, the required data to compute the computation may be allocated in other processor. Unfortunately there is no known way to prevent this. The Vertex sorting technique is combined with the previous one (2D decomposition) in such a way that with the decomposition each processor would have a subset of the Adjacency matrix and by adding an extra phase just at the begging to relabel the vertex, we add extra locality to the edges allocated in a same processor. 3. 3.

Direction optimization . Beamer et al. [7, 6] proposed two new techniques to deal with the unnecessary vertex computations when performing the distributed BFS algorithm. As an introduction, the usual way to transverse a graph is called Top-down: the adjacency matrix is used in such direction that provides us with the neighbour nodes to the ones in our frontier queue. One characteristic of the Top-bottom algorithm is that performs bad with high degree graphs (graphs with many neighbours in one vertex). to avoid unnecessary computations. the Adjacency matrix is duplicated to be able to use the opposite direction to traverse. This way when the number of neighbours is greater than a $\delta$, the number of nodes in vertex in the opposite direction is computed and the minimum among them is chosen ase direction to follow. 4. 4.

Data compression Since the latency in the communication between the processors limits the maximum computational speed in the nodes by stalling them without new data, this method intend to alleviate this effect. This is the main goal of this work. 5. 5.

Sparse vector with pop counting. The pop instruction (__popc in CUDA C) are instructions to deal with bitmaps with hardware support. Some new implementations of the graph500 use this optimization. 6. 6.

Adaptive data representation 7. 7.

Overlapped communication The technique of overlapping communication (also called software pipelining) allows to reduce the stall time due to communication. (Described in further detail in section 5). 8. 8.

Shared memory 9. 9.

GPGPU General purpose GPUs are reviewed in section 2. 10. 10.

NUMA aware With NUMA aware specific libraries it is possible to “pin” a full data structure to an specific (core, processor, bank of memory) in systems with the possible penalty of accesses from one processor to a bank of memory belonging to other processor.

2 Other implementations

Hybrid (CPU and GPGPU)

• •

“B40C” (Merrill et al.) [33]

• •

“LoneStar-Merrill” 111http://iss.ices.utexas.edu/?p=projects/galois/lonestargpu [9]

• •

“BFS-Dynamic Parallelism” (Zhang et al.) [65]

• •

“MapGraph” [21]

• •

“GunRock” (Wang et al.) [55]

• •

“Enterprise” (Liu et al.) [30]

• •

SC11 (Ueno et al.)

• •

ISC12 (Ueno et al.) [52]

• •

SC12 (Ueno et al.) [51] 2. 2.

CPU-based

• •

ISC14 (Fujisawa et al.) [60]

3 Compression

In this subsection we present a list of other work on compression applied to Breadth-first Search implementations. We discuss their algorithm and compare this to ours.

“Compression and Sieve: Reducing Communication in Parallel Breadth First Search on Distributed Memory Systems” (Huiwei et al.) [31]

In this work a 1D partition-Based BFS using SpMV multiplications is analyzed and with and without compression on its communication. The approach followed to compress the sparse vectors is to state that they have low cardinality as they are bitmaps. The compression codecs they analyze enter into the category of Bitmap Indexes. The selected codec is WAH (Word Aligned Hybrid) a patented codec with an medium to low compression ratio but on the other hand an acceptable speed. The result of the study shows a reduction in the BFS time of over the 50% thanks to the compression. 2. 2.

“Highly Scalable Graph Search for the Graph500 Benchmark” (Ueno et al.) [51].

The implementation of Ueno et al. had a good position in the graph 500 challenge during June 2012. It was it that year when they introduced compression. In this case the compression was implemented over a 2D Partition-based BFS executing on a GPFPU device. In contrast with Huiwei et al. Uneno treated the data as a sequence of integers instead of a sequence of bits. Also they chose a Varint-based encoder (VLQ) with neither a good compression ratio nor a good compression / decompression speed. The main goal of this algorithm was its implementation simplicity. Once selected they implemented it on both the GPU device and the CPU. Also they reimplemented two more versions, a 32-bit based and a 64-bit.

A not about GPGPU devices is their current lack of native arithmetic support for 64-bit integers. The choose of a full 64-bit compression on the device would have make the algorithm to incur on a performance penalty. For this reason and also to avoid performing compression for low amount of data they implement (as we have also done) a threshold mechanism.

Chapter 3 Problem Analysis

As aforementioned in the Introduction the optimizations presented in this work have been implemented on a previously developed code [25].

In this section, we will analyze the previous implementation and decide (based on the results of the analysis of the code and the instrumentation experiments) what changes will be made to increase the overall performance.

Also, we will review the optimizations already applied to our code based on other state of-the-art research (discussed in Section 2 - Relared work). Next, we will explain how we have measured the data using internal (coded timers) and external instrumentation (profilers). In the case of the external instrumentation we will compare the available candidates and review the chosen one according to our needs.

In addition, we will also include the measured data of our instrumentations. Regarding the communicated data, it will be highlighted what is its type, as it will be useful to understand the effectiveness of the chosen compression scheme.

*Finally, and based on the data of the instrumentation results, we will analyse and decide what approaches are considered to be feasible to alleviate the detected issues and optimize the overall performance. *

1 Our initial implementation

A list of the optimizations implemented on the Baseline version of this work is listed in the table 1.

1 Optimizations table in “Baseline” implementation

The list of optimizations may be seen on table 1.

2 Our “Baseline” algorithm

In our initial algorithm, where the data is partitioned in 2D fashion, each of the processors operates over is own assigned block of the global Matrix. The

As is was discussed in the section General Purpose Graphic Processor Units (Section 4), some requirements are needed to build an optimal algorithm over a GPGPU device.

In this section it will be listed the pseudo-code of the used parallel BFS algorithm implemented over CUDA. This used algorithm is based on the implementation of Merrill et Al. [33, 25] and includes modifications to adapt the code to a 2D-partitioning.

3 General communication algorithm

The code presented in this work makes use of a 2-Dimensional partitioning of the CSR initial Graph Matrix.

The 2D partitioning were proposed as an alternative to the 1D to reduce the amount of data transmitted through the network [63].

Being $G(V,E)$ a Graph with $\lvert V\rvert$ Vertexes and $\lvert E\rvert$ Edges, on an original 1D strategy the distribution unit across the network would be the Vertexes performing a complexity of $O(P)$ where P is the total number of Processors.

On the other hand, on a 2D strategy the transferred unit are the Edges. This leads to a complexity reduction of $O(\sqrt{P})$.

Described briefly, on a 2D partitioning (the used one) there are two phases where each processo operates and transmit part of their Matrixes.

On the first sub-phase the BFS operates over rows and as result performs an “AllReduce” Operation that transmits the Column. This may be viewed on 2 and 3 of Algorithm 3.

The next sub-phase Multicasts rows [25] Line 4 of the algorithm.

4 Data and communications

This section details the types of data that we will use in our communications. These are based on two structures, the first one consists of a sequence of long integers used to distribute vertexes between the nodes 111Update: In our implementation, for the used graph generator and the maximum reachable Scale Factor , these integers are in the range $[1..2^{24}]$ and are sorted in ascending order. Being the smallest, the first one in the sequence. [25].

A second structure is used to prevent nodes of re-processing already visited vertexes.

The Frontier Queue

The frontier Queue is the main data to be compressed. Each of the elements on its sequence is a sparse vector of bits (SpMV). The analysis that we have performed to multiple big buffers of this data shows an uniform distribution, slightly skewed (Section 4).

The Predecessor list bitmap

The bitmap is used in the last call of the BFS iteration is a 32-bit sequence of integers that is used to convert the frontier queue in the predecessor list.

2 Instrumentation

In this section It will be covered the selection of the Performance analysis tool. It will be discussed the used criteria, the needed requirements in our application and some of the available the options.

As mayor requirements it will be set a low execution overhead, and the ability to wrap specific code zones for measurement.

Minor requirements will be the availability of data display through console. This feature will enable the data to be tabbed more easily into diagrams.

• •

TAU 222https://www.cs.uoregon.edu/research/tau/home.php Originally implemented in the University of Oregon. Offers console interface. Narrow code wrapping is not implemented.

• •

Vampir 333https://www.vampir.eu/ Offers graphical interface. We have not tested this option.

• •

Score-P 444http://www.vi-hps.org/projects/score-p/ Offers both graphical (using Cube 555http://www.scalasca.org/) and console interface. Allow direct code zone wrapping.

1 Instrumented zones

Listed below are the instrumented zones corresponding to the second Kernel (timed), which involves the communication, the BFS iteration.

• •

BFSRUN_region_vertexBroadcast initial distribution of the Start vertex. Performed only once.

• •

BFSRUN_region_localExpansion Execution of the BFS algorithm on the GPU device.

• •

BFSRUN_region_testSomethingHasBeenDone Checks if the queue is empty (completion). To do so performs an MPI_Allreduce(isFinished, Logic-AND) between all nodes.

• •

BFSRUN_region_columnCommunication. Performs a communication of the Frontier Queue vectors between ranks in the same column. This steps contains the next one.

• •

BFSRUN_region_rowCommunication. Performs a communication of the Frontier Queue vectors between ranks in the same row.

• •

BFSRUN_region_allReduceBitmap. Merges the predecessors in all nodes and check performs the validation step of the results. Run only once at completion of the BFS iteration.

3 Analysis

In this section it will be analyzed the previous instrumentation and instruction complexity data.

As a result, a solution will be given for each of the issues found. The implementation, results and final conclusions, follow in the next sections of this document.

1 Communication overhead

As is has been mentioned at the beginning of this chapter, our Baseline implementation is mainly based upon a GPGPU BFS kernel, based in the implementation of Merrill et al. [25], with a 2D data partition, which enables the implementation to perform within a multi-GPU platform and in a multiple node environment through Message Passing Interface (MPI).

In the results of the previous work, it was noticed a loss of performance in comparison to the original Merrill et al. [33] implementation. This was despite the application was running on multiple distributed multiGPU nodes within a fast cluster.

As a result of analyzing the reasons of the performance loss, it was decided to mitigate latency generated in the communications through compressing the transfers.

That new improvement is the work presented in this document.

2 Instruction overhead

As we were going to add a compression system which would be enabled at compile time, we still had space to add and test new improvements. The analysis of the code showed for example that the execution using thread parallelism (OpenMP) was slower than the one where the application was running on a single thread. Even though the experiment platforms has capabilities for multiple thread execution.

In a more on detail analysis we saw that some of the code entering in the timed parts of the application performed matrix intensive operations and allowed parallelism by the compilers. We also saw that many scalar improvements could be applied without much complexity.

In the last part we ran the compiler in inspection mode to retrieve the successfully applied optimizations in our Baseline implementation. Also, and to set a reference point with other distributed BFS implementation, we also performed this test on those. The other implementations are: (i) a Graph500 2D partition-based with MPI reference implementation. (ii) and the latest state-of-the-art Graph 500 implementation using GPGPU [51]. The comparison may be seen in table .

Chapter 4 Optimizing data movement by compression

By compressing the data movements in the graph500 execution, we intend to solve the issues detected in the previous section (Section 1).

*In this section we will first start describing some compression-related concepts for a better understanding of the optimizations. We will describe the available compression techniques. For each of them, we will expose how they match the Graph 500 data types. We will also discuss the related and most state-of-the-art compression algorithms (also known as codecs or schemes). After describing the algorithms we will focus on the available open-sourced library choices. Last, we will focus on the integration and discuss the decisions taken here. *

1 Concepts

We list below related compression terms that are used in this work and its related literature [68, 28, 11, 46, 29, 58, 3].

Information Theory

The Information Theory was the result of applying statistical analysis to measure the maximum compression achievable in a communication channel. It was created by Claude Shannon in 1948, in his paper A Mathematical Theory of Communication [45].

In the successive years, the Information theory has applied statistics and optimal codes to assign the shorter codes to the most probable symbols on sequence to increase that way the compression ratio (defined in this section). This ratio depends on some parameters of the data, like its Entropy (also defined in this section).

Information Retrieval Systems (IRSs)

An Information Retrieval Systems (IRSs) is a system which searches and retrieves data in environments of large data-volumes (for example Database or the World Wide Web).

To achieve this, modern techniques like Indexing are used. These techniques make use of specific compression schemes which have required to be adapted to numerical data (indexes). This benefits our Graph 500 application in two ways:

Since our application also uses numerical data - on each communication it transmits a “Frontier queue” (which is a sequence of integers), we are able to use an integer-specific encoding for our purpose. Also, because of the current importance of this IRS systems, this algorithms improves quickly. 2. 2.

These techniques, due to the requirements of this IRS systems, have evolved to offer high speeds with low-latencies. This factor is important since one requirement in our application is the performance.

Indexer

An Indexer is a type of the IRS. As defined before they manage the access to the stored indexes. A part of this Indexer systems is usually compressed due to the high amount of data that they need to store.

Common web search engines like Google111https://www.google.com/ are examples of Indexers. In the case of Google, the indexing system is based in the Variable Byte algorithm (an specific variant called varint-GB) [13] which will be discussed in further detail below.

Inverted Indexes

The Inverted indexes are the most commonly used data-structures (among other options) to implement indexer systems [67]. They consist of two parts:

A lexicon: it is a dictionary-like structure that contains the terms that appear in the documents plus a “document frequency”. The latter is somehow similar to a “ranking”, and indicates in how many documents those terms appear. 2. 2.

Inverted list (also called Posting list), is an structure which contains (for each term) the document identifiers, its frequency, and a list of placements within the documents.

This latter part usually becomes very large when the volume of the data is high, so it is often compressed to reduce the issue. The used schemes are required to meet this criteria:

• •

they need to considerably reduce the size of the Inverted list by achieving a high compression ratio.

• •

the need to allow a very fast decompression to not affect the search speed.

Shannon’s Entropy

As aforementioned Shannon’s Information theory settled the foundation of the data compression. One of the main concepts of this theory: the “Entropy” enables us to know the maximum theoretical limit of an optimum compression. The definition of the concept is given below.

Put simply, Shannon’s Entropy can be defined as the amount of unpredictability of an event in terms of information. In Shannon’s Entropy events are binary.

According to Shannon, in a distribution where a binary event happens with probability $p_{i}$, and the event is sampled $N$ times, the total amount of information we receive is given by (2) 222This expression applies to discrete distributions..

[TABLE]

To illustrate the concept we will use the example of tossing a coin with equal result probability. In this case, because we have no information about the possible result (Each side of the coin has $p(x)=0.5$), the Entropy (Unpredictability) will be maximum according to Equation (2) ($H(X)=1$). Visually, this can be seen333To display Entropy as a continuous function it has been approximated to $H(p)=-plog{(p)}-(1-p)\log{(1-p)}$ [10]. in the Figure 1.

Limits of compressibility

The minimum number of bits per symbol will be given by the Entropy of the symbols in the sequence. Expression (1).

[TABLE]

Compression Ratio

According to the Information Theory [45], the maximum compression resulting from a data sequence, depends on its probability density function. Accordingly, the bigger the predictability in the sequence, the greater the achieved compression.

The resulting compression is proportional to the compression ratio (greater is better). This ratio is shown in the expression (4).

This concept will be useful to measure the impact of the optimization of our communications in the final results.

[TABLE]

“Lossless” and “Lossy” compressions

These terms are related with compression and thus defined below. First, the compression of our data needs to meet these requirements:

the size of the encoded data needs to be lower than the size of the original size. 2. 2.

the encoding algorithm must allow recovering the original sequence from the compressed sequence, without any data loss.

When both previous criteria are met we are talking about “Lossless codecs”. These are the ones used in this work.

In other literature, only the first criteria above is required, and some tolerance is allowed in the second one. These are called “lossy codecs”. They are often used to reduce the size of media data (images, sound and video). These are not being discussed in this work.

Compression algorithms, Schemes, Codecs and Encodings

In this work the compression algorithms will also be named schemes, codecs or encodings. These four terms are equivalent.

Dependency

To perform the compression of a full integer set some of the codecs operate over each integer independently. Others act over a group of integers and return a sequence in which each compressed block have dependencies to the others. These dependencies are called oblivious (independent) and list-adaptive (dependant) in the literature [57, 11]. They also appear as Integer encoders and Integer sequence encoders in the literature [23].

Usually because independent (oblivious) codecs compress each value by its own and discard the information about the rest may miss global characteristics which could benefit the global compression. Codecs belonging to this method are reviewed on this section even though the main focus will be on the dependent (list-adaptive) schemes.

Compression techniques

The compression techniques discussed below are the foundation blocks of the algorithms reviewed in this document.

Next to these techniques are listed others (“Novel techniques and optimizations”) that are usually combined with the former.

As a result, very optimized algorithms, operating over a very specific data domain (integers and integer sequences) are built. These encodings are, as of the moment of writing this document, the start-of-the-art integer compression algorithms.

Dictionary Compression This is the first of the main compression technique exposed in this document and encompass multiple algorithms families.

The dictionary technique is the most general compression method, and can be applied to any data domain. Usually this technique is made by using Huffman or Arithmetical optimal codes. The symbols to be compressed (dictionary elements) are assigned a dense code according to their relative probability of appearance in the document. The higher the probability of the symbol the shorter the assigned code. This technique has the downside of requiring a full recoding if new symbols are added. A way to check the efficacy of this method for our data, is by generating the empirical Entropy (defined next, in this section) of our distribution of numbers. 2. 2.

Prefix Suppression PS (or length encoding) This is the second global compression mechanism described in the document. This technique is applied to data know a-priory to be a integer. Prefix suppression (PS) suppresses the 0’s at the beginning of each element in the sequence. As a direct result, we reduce the sequence’s domain. This technique generates a bit-aligned variable-sized sequence. 3. 3.

Frame of Reference (FOR) The next global compression mechanism described this document is Frame-of-reference. This method was first described in [24]. This compression technique is applied to a sequence of integers. The way it works is by splitting the sequence of elements in $n$ blocks, so that for each block $B$

[TABLE]

This process is repeated for all blocks in the sequence. The size of each of the $n$ blocks $B$ is:

[TABLE]

This technique contrasts with PS in that the generated sequence has a fixed size. This resulted sequence has a byte-level alignment.

This family of algorithms will be the main focus of this work (discussed below) as they give a good compression ratio and low compression / decompression times.

New techniques focused on integer/ integer sequences compression

These optimizations are applied to the previous techniques to improve the algorithm. The technique Binary Packing described below is a full technique per se, and has been used as the core compression algorithm in this work (S4-BP128 codec). Any of these techniques below may be combined.

**Delta compression, differential coding or delta-gaps technique ** The concept is firstly described in [37]. The delta compression, may also appear with the names Differential coding, delta-gaps, d-gaps or $\delta$-gaps in the literature. [68, 28, 11, 46].

The technique works by storing the distances between elements instead of the elements themselves (which have lower size than the main integers).

Regarding the research tendency over this technique, it is changing in the last few years:

• •

the initial work was focused on the gap probability distributions to achieve the optimal compression ratio. For example, depending on the frequency of each term we could adjust the selected Huffman or Golomb’s codings [56].

• •

the most recent work emphasizes in a better compression ratio and lower decompression latency [49]. 2. 2.

Binary packing This technique implements both the concepts of Frame-of-Reference (FOR) and differential coding. If in FOR each integer in a range is coded and then represented in reference to that range. For example, for a block of size 128, and integers in the range [2000, 2127] the can be stored using 7 bits as offsets of the minimum value of the range (i.e. 2000). Using the expression (3) the resulting size can be seen in (2)

[TABLE]

With the technique of differential coding, now the minimum value is coded and the rest of elements are represented as increments $\delta$ to that value.

With this technique there are two main factors which give theresulting size of the compressed sequence. These are:

• •

the number of bits b used to encode the blocks (in the worst case block)

• •

the block length B

Initial schemes of this technique, refered in the literature as PackedBinary [4] or AFOR-1 [14], use a fixed block size of 128 or 32 integers. Future variations of the algorithm make use of variable size blocks (AFOR-2, AFOR-3) [16]. 3. 3.

Patched coding and Exceptions Binary coding may lead at some cases to very bad efficiency, resulting sometimes in compressions needing more bits than the original representation. This is for the example the case of compressing sequences like the following:

[TABLE]

In the sequence given by (3) in case of only needing to compress the 7th first numbers the required bits would be

[TABLE]

Unfortunately, if the 8th number would need to be included in the block, the compression would even result in a size bigger than 32-bit per integer.

To alleviate this problem Zukowski-et-al [68] proposed the concept of exceptions, where the the 8th number in this sample sequence would be stored in a separate space. This is the concept of exceptions, and will encompass all values greater than $2^{b}$, being b the value defined in the Binary packing concept.

2 Compression algorithms (codecs or schemes)

To describe the compression algorithms used in this document, we will group them first, according to their resulting alignment in the compressed sequence. These groups and their features are shown in Table 1 . In addition to this, some terms used here are described in the introduction of the current section.

The main compression groups to be described are Bit aligned, Byte aligned, Word aligned and Dictionary based.

1. A.

Bit aligned algorithms

1. (a)

Unary This oblivious (group independent) codec represents an integer as N one bits followed by a zero bit which acts as a terminator. 2. (b)

Gamma This oblivious codec consists of two parts. In the first it represents an integer as $log_{2}{N}$ bits in Unary scheme. In the second part the $floor(log_{2}{N})$ bits of the integer are also represented with Unary. This scheme is usually inefficient as integers become larger. To avoid this sometimes Delta compression is used in the second half of the encoding. 3. (c)

Huffman-int This scheme is also an independent integer encoder. This is a variation of Gamma in which the first half (the first $log_{2}{N}$ bits) are encoded whith Huffman optimal codes instead of Unary. 4. (d)

Rice This group dependent (List-adaptative) integer sequence encoder consists of two parts. In the first part the $floor(\frac{N}{2*k})$ bits are encoded using in Unary. The secon part uses ($N*mod2*k$) bits for a given $k$. This schema is also known as Ricek in literature. 5. (e)

Golomb This group dependent integer sequence encoder is similar to Ricek. In this scheme $k$ is a power of 2. 2. B.

Byte aligned algorithms

1. (a)

Frame Of Reference Family This codecs are group dependent sequence integer encoders. As described above the FOR technique has multiple scheme implementations. Some of the differences of these schemes are explained here.

1. i.

PFOR [68] Original Patched coding Frame-of-Reference algorithm. In this scheme the exceptions where stored without being compressed including only its relative order in the linked list. 2. ii.

PFOR2008 [64] as a compression optimization it was proposed the compression of the exception storage, using 8, 16 or 32 bits. 3. iii.

NewPFOR [59] The exceptions are stored in a dedicated part of the output. This area is divided in two subareas. Simple-16 is used as exception encoder in one of the areas. 4. iv.

OptPFOR [59] This scheme is very similar to the previous NewPFOR but the size of the exceptions subareas is selected accordingly to achieve better compression ration an lower decompression latency. Simple-16 is also used as encoder scheme for the exceptions. 5. v.

ParaPFOR [5] This is a variation of NewPFOR optimized for the execution on GPGPUs. As an improvement, the compression of the exceptions using Single-16 is moved from the exception part to the original segment. This modifications worsens the compression ratio but leads to a much faster decompression. 6. vi.

Simple-PFOR Fast-PFOR SIMD-FastPFOR [29] In this codec, exceptions are placed in 32 different reserved areas of the output. These are encoded using Vectorized (with SIMD instructions) Binary packing. The main differences in Simple and FastPFOR schemes is in how they compress high bits of the exceptions values. In all other aspects Simple and FastPFOR are identical. Regarding the SIMD-PFOR codec, it is identical to FastPFOR except that its packs relies on vectorised bit packing for some of the exception values. 7. vii.

S4-BP128444This is the codec that has been used in this work. This codec was also introduced by Lemire et al. [28]. It uses the Bit packing technique described above with blocks (B) of 128 integers. For each block, the bit-width (b) ranges fron 0 to 32 bits. The “S4” stands for 4-integer SIMD. This algorithm uses the concept of Metablock to refer to 4 blocks of 128 integers (Partitioned that way to perform better the SIMD instructions). Before Bit packing is performed, the metablocks use some bits to contain two b values as a header. VByte compression is added to the last blocks that can not be grouped in a metapackage. This is the fastest Lemire et al. family of codecs [28] and is the one used in this work. 8. viii.

TurboPFOR A fixed size Vectorized Binary packing is also used to encode the exceptions. As novelty, a flag bitmap is included to map the placement of the exceptions555https://github.com/powturbo/TurboPFor/issues/7. 2. (b)

Variable Byte Code

This popular codec is group independent (Oblivious), as it encodes each value on its own. As a difference with the aforementioned Oblivious codecs (Unary, Gamma and Huffman-int), this one uses byte alignment instead of bit. A benefit of this schemes over the previous ones is that even though the compression ratio is usually worse, the decoding speed is faster [43]. The main algorithm is known under many names: v-byte, Variable Byte, variable-byte, var-byte, VByte, varint, VInt, or VB. An algorithm that enter this category is the Group Varint Encoding (varint-GB) [13], used by Google. Some improved codecs over the initial Variable Byte algorithm will be listed below.

1. i.

varint-G8IU [46] This scheme is based on Intel ™ SSE3 SIMD instructions. By the usage of vectorisation Stephanoph-et-al outperformed by a 300% the original VByte. They also showed that by using vectorisation in the CPU it is possible to improve the decompression performance more than 50%. The varint-G8IU is an SIMD optimized variation of the varint-GB used by Google [29]. It was patented by its authors. 2. ii.

maskedvbyte [40] This version also uses a SIMD instruction set based on Intel ™ SSE2 to achieve the performance improvements over the original VByte algorithm. The codec adds a Mask (Bitmap) to be able to perform the decoding faster with the vectorisation. 3. iii.

maskedvbyte Variable-Length Quantity (VLQ). Only provide in this work to reference the compression added to their Graph 500 implementation by [51]. 3. C.

Word aligned algorithms

1. (a)

Simple** family ** This family of algorithms was firstly described in [2]. It is made by storing as many possible integers within a single n-bit sized word. Within the word, some bits are used to describe the organization of the data. Usually this schemes give a good compression ratio but are generally slow.

1. i.

VSimple 666https://github.com/powturbo/TurboPFor This codec is a variant of Simple-8b. It can use word sizes of 32, 40 and 64 bits. As other improvements it uses a *“Run Length Encoding”*777Note from the author inside the source code files.. 2. ii.

Simple-8b [3] Whereas the Variable-byte family of codecs use a fixed length input to produce a variable sized output, these schemes may be seen as the opposite process, with multiple variable integers they produce a fixed size output. In the case of this particular scheme a 64-bit word is produced. It uses a 4-bit fixed size header to organize the data, and the remaining 60 bits of the word for the data itself. It uses a 14-way encoding approach to compress positive integers. The encoded positive integers may be in a range of $[0,2^{28})$. It is considered to be the most competitive codec in the Simple family for 64-bit CPUs [3]. 3. iii.

Simple-9, Simple-16 [2, 58] These two formats (also named S9 and S16 respectively in the literature) use a word size of 32 bits. The former codec uses a 9-way encoding approach to compress positive integers, but in some cases bits are wasted depending on the sequence. To alleviate the issue the later codec uses a 16-way encoding approach. Both of them use a 4-bit fixed header for describing the organization in the data and 28 bits for the data itself. Due to the limitation that this implies: they are restricted to integers in the range $[0,2^{28})$ they will not be reviewed in further detail as do not meet the Graph 500 requirements for a vertex size (48-bit integers). 4. D.

Dictionary based algorithms

1. (a)

Lempel Ziv Family and PDICT Dictionary algorithms are based in the Dictionary technique described above. Some algorithms that use this technique are the Lempel Ziv algorithms family (commonly named using a LZxx prefix) and PDICT [68]. Their performance compressing / decompressing integer sequences is outperformed by the integer domain-specific algorithms aforementioned. For this reason they will stay of the scope of this work.

3 Compression libraries

For the usage of each of described algorithms in our Graph500 implementation we could have either implemented from scratch each of the above algorithms or either use a state-of-the-art library containing some of those algorithms. As time is a limitation in this work and the compression algorithms are beyond of the scope of our work, we have used several open sourced libraries based on the work of other authors. The selection criteria of each of those libraries is described below.

Regarding these libraries, and because we noted that all of them presented different benefits, we opted to implement a modular library integration, so we were able to easily include (or remove) any of them so we could and test and demonstrate the positive impact of each. The library integration has been made with the following design parameters in mind:

• •

The design is simple and modular. This allows new maintainers to add new desired compression libraries with the minimum effort. This has been done by the usage of a software design pattern called “Factory” [22]. In this way we have managed to improve of three of the Graph 500’s evaluation criteria:

Minimal development time: addition of new libraries by new maintainers in an small amount of time. 2. 2.

Maximal maintainability: modification and debugging of the current code. 3. 3.

Maximal extensibility: addition of new libraries by new maintainers with a minimal impact on the already developed code.

• •

Minimal possible impact in the performance of the BFS kernel. We have managed to achieve this by creating the Factory call prior to the initialization of the second Kernel (BFS cycles and validations) and passing a memory reference to the compression / decompression routine.

1 Lemire et al.

This library888https://github.com/lemire/FastPFor/, has been created as proof-of-concept of various academical works [68, 28, 24]. One of the main contributions of this work is the achievement of a very high performance and compresion-ratio. This has very high value in the compression of Inverted indexes (Section 1), which are very used nowadays (as it has been exposed before).

Lemire et al. achieved this high performance by implementing a version of the algorithm using SIMD instruction sets in the CPUs.

The base of this work is the Frame-of-Reference algorithm first purposed by Goldstein-et-al [24], and later improved by Heman-et-al who among other additions to the algorithm, also enhanced the performance by using Superscalar techniques [68].

Some characteristics about the Lemire’s implementation are listed below.

Compression input types and sizes

This library allows only unsigned 32-bit integer input type.

Compression output

The output type is the same as the input: unsigned 32-bit integer.

Included codecs

FastPFOR, SIMD-FastPFOR, S4-BP128(Used), Vbyte, maskedvbyte, Varint, Varint-G8IU.

Pros and Cons

Good compression ratio for the best case scenario. A 32-bit integer is reduced to near 4-bit (value near to the Entropy) ($ratio=\frac{32}{4}=8$). Also, the library has been implemented in a clean way and the source code is easy to understand.

As disadvantages to fit this library in a Graph 500 application, we have found the next:

The implementation restricts the size of the compressed 32-bit integers. This limitations implies:

• •

for each compression-decompression call, a size and type pre-conversion must be carried out carried out to adjust the type to the required by the Graph 500.

• •

to avoid integer overflows in the previous conversions, it is required to perform a pre-check over all the integers in the Graph just after the graph generation.

2 Turbo-PFOR

The Turbo-PFOR999https://github.com/powturbo/TurboPFor library is an open sourced project containing a full set of state-of-the-art integer compression algorithms from various authors. The library main contribution is an improved version of the Lemire et al. [28] implementation, mixed with the ideas and code from the other authors. The Lemire’s code has been fully optimized101010https://github.com/powturbo/TurboPFor/issues/7 and the library itself makes a massive use of the CPU’s SIMD instructions.

The package also includes an statistical Zipf Synthetic Generator which allows to pass the skewness of the distribution as a parameter. The statistical generator allows to include into the compared terms the value of the empirical Entropy of the generated integer sequence.

The package also compares the new optimized implementation with other state-of-the-art from other libraries. As the result of the optimizations, the author claims a better compression ratio and a higher decompression speed for integer sequences more similar to phenomenological Pareto distributions.

Compression input types and sizes

The library allows sizes of 32 and 64 bits integers without sign.

Compression output

The implementation relays on the C Standard Template Library (STL) and the output type is the same than the input.

Included codecs

“Turbo-PFOR”, Opt-PFOR, Lemire et al.’s set of codecs, Simple-8b, Simple-9, Simple-16, LZ4, Vbyte, maskedvbyte, varint-G8IU.

Pros and Cons

Similar latency than the analogous codec in Lemire et al. library, with an 2X compression ratio improvement.

Also, the integer input size matches the “64-bit unsigned integer” required by Graph 500. This prevents our code from the pre-conversion step in lines 2, 4, 6, 8 of algorithm 5.

As a possible downside, the code is difficult to understand since it heavily relays on C macros.

3 Alenka CUDA

The Alenka111111https://github.com/antonmks/Alenka project is aimed to be a fast database implementation with an inbuilt compressed inverted indexed based on CUDA and GPGPU technology. The query implementation is also based in this latter technology.

The main compression algorithm is an state-of-the-art Patched Frame-Of-Reference Para-PFOR (Section 2). Its main features are:

Compression input types and sizes

One of the main benefits of this implementation is that implements a compression for unsigned integers. These size and sign suit the Graph 500 requirements which are “unsigned 64-bit integers”. Therefore it is not needed the pre conversion step in **lines 2, 4, 6, 8 of algorithm 5 **. The possible integer input types are 32 and 64 bits without sign.

Compression output

the algorithm generates a data structure containing information to re-ensemble the data in the other peer.

Included codecs

Para-PFOR implementation over CUDA (Section 2).

Pros and Cons

As an advantage, the code is clean (easy to read) and short. As a possible downside may be pointed that the transformation of the output structure to an equivalent MPI structure to perform the communication may require some extra programming effort. However, this coding is mostly complete.

4 Integration of the compression

As described in further detail in other parts of this document, the integration of these compression libraries (which combine a broad spectrum of HPC technologies) has been designed flexible (“pluglable”) so that it is simple to add a new “library” from other authors. This library addition would involve some steps like:

The 3rd-party compression library need to be copied to the Graph 500’s compression/ folder. 2. 2.

The “Factory object” must become aware of that a new library. This is done by selecting a name for the new class within that object and adjusting the input and output types. For example, it would be correct to copy an already existing one, choosing a new name, and setting the types. 3. 3.

A compression-specific Object for that new library must exist. This Class contains a wrapper of the new library’s compression & decompression calls so they match our application’s calls. An specific file, that can be copied and modified exists for this purpose. In this new created object it is required to adjust the Class’s name and the new types. 4. 4.

Finally, the file containing a mapping between that compression type and its MPI equivalent needs to be updated.

Other aspects regarding the integration are the low impact of the integration of the integration on our code. The code performs only one call to create the compression object, which remains in memory. This is performed outside the Second kernel (BFS + Communication) to minimize its effect on the the overall performance. The compression object is passed as a memory position to the BFS calls.

Other benefit of this implementation is that the compression library may be updated regardless from the Graph 500 implementation.

1 Experimental performance comparison

In order to select the compression scheme that best suits our purposes, a comparison table (Table 4 and Table 5) containing the results of our compression experiments is presented below. In these tables, the parameters used in the headers are described below.

• •

C. ratio % shows the achieved compression-ratio represented as percentile. The lower the better. The formulat can be seen in expression • ‣ 1

[TABLE]

• •

Bits/Integer Resulting bits required to represent the 32-bit integer. This value has as lower limit the value of entropy (Section 1) of the sequence distribution. The lower the better.

• •

C Speed MI/s represents the compression speed measured in Millions ($10^{6}$) of 32-bit integers processed per second. The higher the better.

• •

D Speed MI/s analogous to the previous parameter. For decompression speed. The higher the better.

The aim of generating the data is to be able to gather compression statistics from each codec and see the one that best suits our needs: we require good decompression speed (like an IRSs), but also a good compression speed, to not penalize the compression calls. A good compression-ratio is as important as the two previous aspects. For building the table we will be using a synthetical Zipf integer sequence generator. The generator is specific for generating sequences of indexes (integers) similar to the ones appearing in web pages, search engines indexes and databases indexes (Further discussion about those topics in the next sections of this chapter).

Regarding graphs datasets, they have been study to follow Power law distributions (like Zipf, pareto or lognormal) [54]. So, that could be a clue of that using a Zipf generator would give a real insight about the data we intend to compress. Unfortunately, regarding our data, we do not compress the vertices. Instead, we do compress sparse vectors (SpMV) represented as integers. These integer sequence have some own numerical characteristicswhich are going to help us to choose the correct codec for of Frontier Queue compression.

If we could do a summary of what family of codecs to choose depending of the properties of the data to compress, it would be as follows:

• •

For compressing an unique integer instead of a sequence, our option would be a Prefix Suppression (PS) type codec. Whiting the available codecs in this broad family we would choose the one which better adjust to or latency , speed and compression-ratio requirements.

• •

If we had a sequence of integers where the distance from one number to the other is a-priori known to not be very high, the Frame Of Reference (FOR) type with a possible delta compression on top would be a good option.

• •

If the data is an integer sequence, where the numbers are known to be low (applying this family over high numbers could result counterproductive), then the Varint family could be a good option,

• •

If the cardinality of the data to compress is low, for example bitmaps with ’0’s and ’1’s, a Bitmap Index family would be a good option. Downsides of bitmap indexes are their low speed and their heavy memory footprint.

• •

for any other case where we have no knowledge about properties of the data to be compressed a general Purpose dictionary compression (for example Lempel Ziv family) would be a good choice.

To gather information about the sequence of integers to be compressed (the SpMVs from our Frontier queue) we have manually extracted several big buffers (bigger than 20M 121212M stands for Millions integers) into files and analysed them afterward with an statistical package [15]. The result of the tests showed an slightly left-skewed Uniform distribution (Figure2). In terms of compressibility the term Uniform distribution is similar to randomness. Because random sequences have a high entropy (Section 1), the compression would most likely be bad. Even though the entropy of those sequences was high: 15-bit out of 32-bit of an integer, we have been able to perform a very high compression (over 90% reduction) very close to the entropy value. This has been thanks to the fact of (1) the data was sorted due to our BFS implementation (2) the gaps in the sorted sequence where small. (3) this kind of sequences have raised interest due to the Big data + Indexer phenomena and have the support of efficient algorithms nowadays.

Also, as aforementioned we are including in this work a compression based on GPU (which we believe can provide interesting performance benefits in comparison to the two ones presented in this work - based on CPU’s SIMD instruction-sets). Regarding compression implemented on GPUs, it has to be noted than currently, due to the limitation of a maximum integer size of 32-bit in the GPUs cores’s ALUs, the use of GPU compression on a graph 500 application would limit the maximum achievable scale of the BFS to 32 ($2^{32}Vertices$). The option of using a 64-bit integers with non native hardware support would penalize the compression. [51].

In order to test compression our chosen Frame of Reference codec in GPGPU we are configuring our benchmarking tool to use a data corpus 131313TREC GOV2 corpus without sensible data, and valid licence for this work. http://ir.dcs.gla.ac.uk/test_collections/access_to_data.html used in papers from other authors, where the algorithm for our CUDA compression (Para-PFOR) is tested with the same dataset and on a similar platform. Our used system is a node of our “Creek” cluster (Table 1).

The result of this comparison can be seen in Table 5. As a first result it can be noted that the decompression speeds of the Para-PFOR algorithm on a Fermi architecture are faster than any other SIMD techniques implemented on CPU. Regarding the compression ratio we can observe that it is lower than other more optimized algorithms. The reasoning behind this is that the Para-PFOR scheme sacrifices its exception management to benefit the overall decompression speed [5].

Our results of this theoretical simulation141414Based on a comparison with other existent research [5] on a comparable environment. of the CPU and CUDA compressions are the following.

• •

the performance of the Para-PFOR using GPGPUs performs better than other SIMD based algorithms, but lacks of a good compression ratio. The latter disadvantage may be alleviated by using the exception management of other PFOR codecs.

• •

the usage of the compression directly in the GPU device may result interesting because it may minimize the memory footprint of the used data structures of our BFS (in case of used encoded buffers directly in the card). Unfortunately, even though we would be able to use more data of the graph in the limited GPU devices’ memory, we would still be limited to an scale of 32, if we want to keep the performance of the compression algorithm low (no native 32-bit integer support on GPU devices [51]).

2 Observations about experimental performance

As it can be seen in the figures of the previous results, the described above algorithms are focused for their usage in IRSs and therefore designed to achieve a high compression ratio with high decompression speed (Table 4).

The higher importance of the decompression speed over the compression one is founded in that these systems need to access (decompress) inverted indexes at near real-time with optimally almost zero latency.

In our integration of these integer sequence compression schemes for a graph 500 application we should focus on the schemes offering a good compression ratio and balanced good $\texttt{compression}\longleftrightarrow\texttt{decompression}$ speeds.

In this work three compression libraries have been integrated (two of them partially). Each of them features a main optimised codec. Based on Table 4 and Table 5 (resulted from comparing results from similar test environments), we will give an insight on how each codec fits in the Graph 500 application.

Turbo-PFOR “Turbo-PFOR” This library has only been partially integrated. The included codec Turbo-PFOR offers a near 50% better compression ratio than any other in this work. The decompression speed performs as well as SIMD-FastPFOR (Lemire et al. library). As downside the compression speed is lower. We believe this scheme may offer promising results. 2. 2.

Lemire et al. “S4-BP128” This is the library used and integrated in this work. Its most competitive codec S4-BP128-D4 (tested in this work) offers both a good compression ratio and good (and balanced) compression and decompression speeds. 3. 3.

Alenka CUDA “Para-PFOR” The partially integrated Alenka library, includes as main algorithm Para-PFOR implemented over CUDA. This codec has been described, and also theoretically compared here. The comparison has been made by adapting the data from other works (a dataset151515TREC GOV2 corpus without sensible data, and valid licence for this work. http://ir.dcs.gla.ac.uk/test_collections/access_to_data.html has been executed in a system with similar capabilities than the one in the paper). As a result this codec seems promising because with a close compression ratio to the one used in this work (S4-BP128) it outperforms all codecs in terms of decompression speed 161616This result has been calculated using theoretical data.. The compression speeds are not listed in the reviewed technical document and therefore unknown.

3 Usage of compression thresholds

By compression thresholds we refer to the minimal size of the integer sequence, required to trigger the compression call. The motivation behind this concept is that the compression routines add an extra overhead to the Processor Unit (CPU or GPU, depending of the used library). This overhead is sometimes not justified in terms of benefit from the compression.

This is the case of very short sequence of integers where the time of performing $\texttt{Compression}\longrightarrow\texttt{Transmission}\longrightarrow\texttt{Decompression}$ is greater than the one achieved by just Transmission.

Regarding this optimization, we have not noticed significant benefit in its use. However, it will be discussed in section Future work (Section 8) this technique may help other improvements.

4 Communication before and after compression

In this section it will be described the algorithms of the communication with the integrated compression call. As it was described in the section Graphs 500 Challenge (Section 3), the specifications of the Graph 500 benchmark specify Vertex labels of 48-bit unsigned integers. To be compressed, our integer sequences (defined as 64-bit integer) we have proposed three codecs (i) S4-BP128 (used), (ii) Turbo-PFOR and (iii) Para-PFOR (over CUDA). The first one operates over 32-bit integers. The last two codecs operate over 64-bit integers.

To be able to fit the 64-bit integer sequence, on the used compressor (S4-BP128, 32-bit), we have favourably used the fact that our Graph generator did not overpassed the 32-bit range $[0,2^{32})$, due to the Scale used in the experiments. As a result of this the 64-bit integers could be pre-converted to 32-bit and post-converted again to 64-bit.

The synthesis of the compression-integration is detailed in the Algorithm

Chapter 5 Optimizing instruction overhead

*As a second group of optimizations made to our BFS base implementation we will be testing some optimizations based upon the efficient programming of our algorithm. In these optimizations we will be covering five areas: (i) scalar optimizations (ii) vectorization (iii) thread parallelism (iv) memory access, and (v) communication patterns. For each of these subsections, small code snippets illustrating the optimizations will be given. At the end of the section, it will be included a checkbox table showing which optimizations have been added. The results are listed in the Final Ressults section (Section 6). *

1 Scalar optimizations

Scalar operations / data are those which involve common data types or functions where the data consist of a single value which can not be operated in parallel. Examples of these in C language are an integer, a char or a data struct, and also any operation involving the use of any on them.

In this section, some methods to improve the performance with scalar operations will be discussed. On each technique it will be annotated whether (or not) that technique has been used in the Graph 500 implementation, and where.

1 Strength reduction

This technique has been intensively used in our code. By Strength reduction we cover transformations in the data and functions that avoid either unnecessary value conversions (or floating point operations), and expensive arithmetical functions by the usage of its reciprocal. Some sample transformations are listed below using short code snippets.

Usage of the correct value for its specific data type. In C language this would encompass the usage of the suffixes “L”,“U” or “LU” for Longs, Unsigned values or both the previous.

This will also include the case of the suffix “f” for a float value. By doing this and keeping a coherence in the arithmetic or floating point operations, unnecessary conversions are avoided. 2. 2.

The replacement of arithmetical operations (with big overhead) by a reciprocal (with less overhead) will also improve the results in the case of operations in critical loops. This is the case of replacements of Divisions (slower) by multiplications (faser).

Replacement of expressions with much overhead with cheaper ones

By expensive expressions we mean long operations involving some of the previous division/multiplication or transcendental functions which can be partially pre-computed and transformed manually outside a critical loop, and substituted with an equivalent expression.

Algebraic transformations. Usage of equivalent functions with less overhead

The usage of transcendental functions with hardware support instead of their analogous. replacement of exp(), log() or pow(x, -0.5) with exp2(), pow2() or 1.0/sqrt(). This kind of strength reduction can also accelerate the global result when the operation is performed inside a critical loop.

Complexity reduction in loops by the use of caching

This technique, also called loop hoisting or loop invariant code-motion in the literature can be placed also in the next subsection memory access optimizations.

This technique involves loops consisting of several nested loops. This kind of transformations are usually detected and performed by the compiler. In some cases where the invariant variable is a vector, the compiler may skip it.

In this optimization, a memory position (e.g. a vector) indexed with the outer loop index does not need to be re-accessed in the inner loop. This variable can be copied to a scalar variable (“cached”) outside the inner loop and therefore prevent the positioning within the vector.

Efficient bitwise operations

The optimal bitwise transformations have been used intensively in this implementation. Bit-level operations are faster than their higher level analogous. Also, the usage of optimal methods which prevent branching whilst the bitwise operation is performed makes the result more optimal.

Control of floating point semantics

These optimizations (also called Precision control) are usually activated in the compiler and therefore would also fit in the subsection about compiler optimizations. However, because they enable a balance between precision and speed of floating point scalars are placed in this section.

The control is done at compile time and allows the selection of multiple levels of precision. The lower the level of precision in the Arithmetic Logic Unit, the faster the resulting code will be.

Because the Graph 500 test is not a floating point intensive benchmark, this implementation, as well as others like Titech’s [51] use the lowest possible precision mode to gain performance.

2 Vectorization

As previously defined in the section Scalar optimizations a variable based upon multiple data values allows the possibility of being vectorized. This is for example the case of an array in C.

By saying that the variable can be vectorized we mean that in case some requirements are met, the variable will be processed by a Vector Processing Unit (VPU) within the CPU.

This will be done with the help of specific Instruction-sets (SSE, AVX, etc in x86 architectures or AltiVec in PowerPC ones). These instruction-sets are an example of SIMD designs described in Section 2.

When a vector is processed in SIMD fashion (by the VPU), as opposed to when it is processed as an scalar by the Arithmetic Processing Unit (ALU), the calculations are accelerated by a factor of the number of values fitting the register. As it can be seen in the table 1, an Intel™ AVX instruction (256-bit width register) operating over 32-bit integers will enable an assembly instruction to operate over 8 integers at once.

A way to look at this is, if the code is not using VPUs where it could use them, it is performing only $\frac{1}{8}$th of the possible achievable arithmetic performance.

To manage the usage of the SIMD instruction-set there are two possible approaches.

The usage of in-line assembly or special functions in higher level languages called intrinsics. 2. 2.

Regular code in a way that allows the compiler to use vector instructions for those loops. This approach is more portable between architectures (requires a re-compilation). Also, it is more portable for future hardware with different instruction-sets.

The vectorization made in this work focuses on the second approach: automatic-vectorization by the compiler.

As the concept of aligned buffer is used multiple times in this sections we define this below.

[TABLE]

Where T is the type of the pointer and p contains its memory address.

Optimization of critical loops

This is one of the main focus on instruction overhead reduction in this work. An introduction to this technique has done in the Vectorisation section above.

As it has been described, the usage of this technique over code with many loop iterations and array operations, may get a higher benefit from this. However, some requirements need to be met for the automatic-vectorisation being performed. These will be explained later in this section.

The main benefit of the usage of this technique on our application is in the second Kernel (BFS + Validation), and more specifically on the validation code.

Our validation code is an adaptation of the reference code provided by the Graph 500111http://www.graph500.org/referencecode community, with the addition of a 2D partitioning [51]. The matrixes operations are performed by the CPU and iterated with non optimized loops which miss the chance of any possible vectorization.

These optimizations based on array vectorisations have been proved to be successful on other state-of-the-art implementations like Ueno et al. [51]. Also, in other works not directly related with the Graph 500, which also makes intensive usage of arrays (the compression library [28]) this optimizations have resulted successful.

In order to successfully manage automatic-vectorisation the requirements listed below must be met.

• •

The loop must be a for-loop. Other kinds like while will not be vectorised.

• •

The number of iterations of the loop must be known either in compile or run time.

• •

The loop cannot contain complex memory accesses, i.e. memory accesses must have a regular pattern (ideally by the use a unit stride. e.g. V[n+STRIDE]).

Note that some compilers may manage these issues, however this is not the general rule.

• •

The loop must not contain branching (for example if-then-else statements).

• •

The data to be vectorised must not depend on other index pointing to itself (for example, a compiler will not vectorise data in the form V[i]= V[i-1]. Sometimes even though a dependency may exist in the code, the compiler can be made aware to ignore it with the use of language keywords, like for example restrict in C)

In addition to these restrictions which applies to the loop, the compiler must be aware of this action by the use of a -O flag. Also, the internal structure of the data vector (this will be explained in further detail later in this chapter) must be declared “aligned” in memory. Only when the code is programmed accordingly, the compiler will detect and automatically vectorize the loops. Even though these improvements may have impact in overall results , due to time constrains, they have been implemented partially. This is reflected in the section Future work (Section 8).

In other implementations of Graph 500, [51], the vectorization of array structures executing on the CPU has been performed on detail. For that matter a fixed alignment of 128-bit (would allow Intel ™ SSE* instructions on x86 architectures) is used for the general memory data buffers. Also, other data structures which require temporal locality use a different alignment constants(“CACHE_LINE_ALIGNMENT”, set to 32-bit for 32-bit integers). These tasks are performed on compile-time by the use of wrappers over the functions that allocate the memory.

Regularization of vectorization patterns

In auto-vectorisation, even though the compiler detects loops with incorrect alignment, it will try to vectorize them by adding scalar iterations at the beginning (called “peel”), and also at the end (called “tail”). The rest of the vector will be vectorized normally. To perform a regularization of the vector a padding is introduced through a second container. Also, the compiler is made aware of this not required operation by the usage of a “pragma”.

This optimization could be done in our Graph 500 application to increase the amount of vectorised loops. These transformations have not been implemented and are listed here as reference.

Usage of data structures with stride memory access

When a data structure must be designed to solve a problem, sometimes two approaches are possible. Both are illustrated with code snippets below.

1 Compiler optimizations

To perform the optimizations, one possible strategy is the use of a preprocessor tool that prepares the build based on selected parameters, or detected features in the host machine’s hardware.

The selection of build tool was difficult as currently exist multiple, and very similar options (e.g. cmake, Automake, etc). In our case we selected our decision on the ease of porting our current Built scripts to other tool. For this matter we chose Automake, a mature project with active support developed by GNU.

Apart from being able to select options in our build or automatically detect optimal hardware features, it is possible to select good build parameters based on good practices. Those are exposed below.

Usage of an unique and optimal instruction-set

As it can be extracted from white papers of commercial or open source compilers[42], a good practise to follow when creating the binary code for an application is to specify the target architecture on which that code will run. This will enable the compile to select that known applicable software support for the specified hardware. Also the avoidance of including several instruction-sets which will increase the size of the generated binary code may impact on the overall performance.

As it has been specified in the introduction of this subsection, by the usage of a preprocessor tool, several tests are done before the build to detect the host platform. The assigned architecture will depend on this result.

The Table 2 compares the used compiler parameters for different Graph 500 implementations. The explicit usage of a pre checked target architecture may be seen under the column march in that table.

Inclusion of performance oriented features

Some research has been done on the benefits on applying Inter-procedural optimizations (IPO) to array-intensive applications to improve the general cache hit rate [66]. By contrast with other implementations like [51] this implementation applies this technique.

Other performance oriented parameters added to the compiler are the -O and -pipe flags. In Table 2 can be seen a comparison of the different configurations of common Graph 500 implementations.

Removal of non performance-oriented features

In contrast with the previous compiler parameters, other parameters are known to reduce the overall performance [39]. The removal of the -fPie flag from our previous implementation may have influenced in the performance increasement.

Ensure auto-vectorisation is enabled in the compiler

All the requirements listed in the section Vector optimisations are enabled by making the compiler aware of this. Usually vector optimisations are performed when the -O level is greater then a value. The value depends on the compiler. For the specific case of GNU’s gcc compiler this value equals “2” (-O2).

Ensure auto-vectorisation has been successfully enabled

The Compiler reporter is a reporting mechanism that outputs all the successfully performed optimizations automatically. It is usually activated by a flag.

By checking that the code transformations that have been performed we can confirm the changes in the generated code.

Due to time restrictions this testing was not fully performed in our Graph 500 application, and the part of Vectorisation is only partially verified.

3 Thread parallelism

The previusly described optimizations act over the achievable parallelism in the data layer. In this new layer it is possible to generate serial codes executed concurrently by the creation of “threads”. In HPC, this threads make use of the parallelism that introduces the multiple cores within the processor.

These threads must be used with a parallel Framework. Some example implementations of parallel threads Frameworks are Pthreads, OpenMP, Cilk+, Boost threads.

We will focus this work on OpenMP as it is the most popular and standardised framework, with the most active support. It implements one of the most common parallel patterns in HPC, the parallel-for or Map-reduce model. On it, when a for needs to iterate n items and uses this model, the framework distributes the items by X.

Concurrent paradigms may increase the overall performance. However, they have the downside of an added programming complexity. The writes and reads in shared memory by multiple codes may create “race conditions”.

Put simple a race condition is a non-deterministic state where the multiple processes have own copies of a variable and operate over it. When at some point the variable is copied to a global position, different threads will change the value and it will be the value of the slowest process, the one remaining in memory.

In our Graph 500 implementation we were required to constantly test the application when the OpenMP threads was re-adjusted.

For a better understanding of the section Thread contention below, the concept of synchronisation of threads, and solutions for the race condition issue will be defined next. The problem of race conditions is usually solved “locking” one or more threads. This is called synchronization.

1 Usage of fork-join model when suitable

In contrast with the parallel-for (or map-reduce model), the OpenMP framework allows also the usage of a very mature model based on Parent and Child processes.

At some cases when it is not possible the usage of a Map-reduce model, it is recommended to study the possibility of including this one.

2 Optimal thread scheduling

The threads framework OpenMP, in the Map-reduce model (parallel-for) allows the implicit specification on how the load of tasks per thread will be balanced. This is done by specifying a scheduling policy. The three main policies available are listed below. Also, it is specified the one that we consider the most suitable.

static(n)

this is the default policy. Each thread has a fixed number n of tasks to perform. Even though other threads have already finished their tasks, if one task in one thread takes longer to run, the other threads must wait.

dynamic(n)

this policy assigns dynamically the tasks according to the load in other threads. One downside of this policy is that threads may not have enough temporal locality resulting on bad cache hit rates.

guided(n, min)

this is the policy chose in our implementation. It mixes both previous scheduling policies. It is basically a dynamic policy with a minimum number of fixed tasks in the thread. The temporal and spatial locality of the accessed data may result in a better cache hit rate.

A note about the scheduling policies is that the correct selection is better done when the execution has been analysed with a profiler.

3 Thread contention prevention

The threads contention may have two main causes. The first one is an excess of synchronisation (defined above). The second one is called false sharing and is defined next.

As it was explained before an address is n-aligned when it is multiple of n. In the case of false sharing the retrieved data of some threads (even though is aligned) is partially divided between two cache lines. This can translate into overhead due to complex memory accesses. In the case of a big-sized critical loop, this alignment issue may affect the maximum performance benefit resulting from the use of threads.

The issue is solved by adding a padding which corrects the alignment of the vector.

This optimization would be recommended in the critical loops of the validation code.

4 NUMA control / NUMA aware

In Shared memory systems, platforms with multiple processors have access to different banks of physical memory. The banks which physically belong to other processor are called “remote memory” for a given processor. The remote memory access have a greater latency penalty compare to local memory access (This can be seen in the Figure 3). This same effect happens in multicore systems and each core’s first level cache.

During the initial assignation of threads in this NUMA model, a thread running in one processor may be using the memory of a different one. All of this results in the threads adding overhead to the application. This issue affecting Non Uniform Memory Access (NUMA) architectures can be alleviated with several approaches.

• •

Because memory allocation occurs not during initial buffer declaration, but upon the first write to the buffer (“first touch”), for better performance in NUMA systems, data can be initialized with the same parallel pattern as during data usage. This approach can be seen on Figure 10. By using this technique an initial assignation is done using the threads framework - this prevents memory-access issues in the threads.

• •

Other approach is the one used by the current Top 1 implementation in the Graph 500. They have implemented their own library 222ULIBC - https://bitbucket.org/yuichiro_yasui/ulibc which allows to “pin” any data structure from memory to any specific CPU, memory bank o core.

In our work, after the analysis of the results we have included in the section Future work (Section 8) a revision of our thread scheduling policy, as well as the introduction of this technique to alleviate the bad performance when the thread execution is active.

4 Memory access

The memory access optimizations performed by alignments and more sequential structures, (as it was stated previously in the document), will help in the vectorisation of many loops. In this section a different approach is defined. This technique helps also with critical loops and may be an advantage in some zones of our second kernel (BFS + verification).

1 Locality of data access

The “Rule of Thumb” for memory optimization is the locality of data access in space and in time.

• •

By spatial locality we mean a correlative access to the data in memory. This is managed by a proper election of the used data structures. The use of Struct of Arrays (SoA) over Arrays of Structs (AoS) improves this situation. The proper packing of the data to avoid having to iterate in 1-stride fashion at the beginning and end of the structure (padding the structure) also improves the locality.

• •

By temporal locality we mean that the required data at one point is as close as possible to the one in the next point. We manage this by changing the order of the operations (e.g., loop tiling, loop fusion or loop hoisting).

2 Merge subsequent loops

In this operation in different but subsequent loops (with same iterator distance) may be groped if the flow of the program allows this. This produces that the retrieved memory block to the CPU cache can be reused more close in time. This transformation can be seen in Figure 12.

3 Variable caching in the outer loop

This optimization (loop hoisting) has been described in the previous section Scalar Optimizations. It is listed here as it is also a memory-access based optimization. As previously stated, this optimization has been included in many parts of the code (both kernels of the Graph 500 application). This transformation can be seen in Figure 14.

5 Communication patterns

Previously in this section, the optimizations over the multiple data in one instruction was defined as the first layer of parallelism. Next, in the case of the threads which are generally used in HPC to control the cores inside a processor, it was defined the second layer of parallelism.

In this section it will be defined a third layer. On it, the “parallel processors” are each one of the nodes communicating by message passing (MPI).

1 Data transfer grouping

Small and subsequence data transfers can be prevented by grouping them in a MPI custom datatype. This improves the communication in two ways: (i) a longer usage of the network, and hence a better throughput is achieved (ii) the memory access improves when relaying on the internal MPI packing system (Zero copy).

We have applied these optimizations several times to alleviate the overhead of successive MPI_broadcasts and MPI_allgatherv for the pre distribution of communication metadata (i.e. buffer sizes in the next transfer)

2 Reduce communication overhead by overlapping

The usage of this technique implies the continuous execution of tasks during the delays produced by the communication delay. To perform this optimization an asynchronous communication pattern, and a double (Send and Receive) buffers must be used.

This optimization has been proposed as other possible improvement in the communications and is listed un the the section Future work (Section 8).

6 Summary of implemented optimizations

In this section we have covered many optimizations, from which only a subset have been added into our work. Mainly, the added optimizations have been as follows:

Scalar optimizations. From this group we have added all of the listed optimizations and covered most of out code. 2. 2.

Vectorization. The optimizations covered in this section been the main aim of our instruction overhead reduction. These optimizations have been added with more or less success (this will be discussed in the section Final results, Section 6) 3. 3.

Thread parallelism. As we have mentioned above, to deal with threads we have used the OpenMP framework. When first analyzing our Baseline implementation we noticed a paradoxical performance loss when using threads compared to the 1 single thread case. This problem has been constant during all the development of our new 2 versions. In the last version we have achieved to slightly increase the performance of the OpenMP version over the 1-threaded one. 4. 4.

Memory access. The memory access optimizations usually require the refactoring of large portions of code. For this reason sometimes it becomes more difficult the optimization of an already created code the its new reimplementation. Hence, we have only applied a few of this optimizations to the new structures that have been required to implement. 5. 5.

Communication patterns. Among these optimizations we have only covered the regrouping of individual transfers performed serially into a custom MPI_Datatype which would achieve in one transmission a bigger amount of data and a bigger use of the network resources. This modifications have been made majorly for the broadcasted metadata prior to other communications. Regarding the communication overlapping it has successfully proved to be a successful addition to a Graph 500 application [60, 51] (Section 2) Despite of this, we have not implemented it in our code.

Chapter 6 Final Results

In this section we present the results from our experiments and try to solve the problem detected in the previous section Analysis (Section 3). In this same section we will discuss the results and set the initial thought for the next sections, Conclusions and Future work.

The chapter starts with data related with the systems used for testing. One is allocated in the Georgia Institute of Technology; other, is in the Rupretcht-Karls Universität Heidelberg.

The results will be grouped depending on the type of performed optimization: (i) communication compression or (ii) instruction overhead optimizations.

1 Experiment platforms

The Tables 1 and 3 show the main details of the used architectures. The big cluster (Keeneland) has been used to perform main tests. The development and analysis of the compression has been performed in the smaller one.

A third table 2 shows the specifications of the platform (also allocated in University of Heidelberg) used to test the overhead added by the compression. The experiments performed on this third platform have not been included as they have only been concluded partially. Even though, some conclusions of the executions on this machine are are added in this chapter. The table is therefore added, as a reference if needed.

2 Results

For this work we are going to focus in three main points of our development. This will help us to identify our two main goals: compression and instruction overhead. They are described below.

Baseline was our first BFS implementation. It its studied in the section Analysis. The Baseline implementation does not compress data and is a good scenario to study traditional overhead reduction on a graph500 application 2. 2.

No-compression is our latest version without compression. The later has been disabled at compile time and introduces zero overhead. Measuring the difference in performance between this and the previous one (none of them have compression capabilities) is a good chance to measure the instruction-oriented improvements. 3. 3.

Compression same code the prior one. It has its compression capabilities enabled.

Here follows a description of the diagrams listed below. The diagrams have been generated from the trace files resulting from the executions. The three main types of graphs shown here are (i) strong scaling diagrams (ii) weak scaling (iii) time breakdowns. To know what conclusions we are getting from each type of graph they will be described.

Strong scaling keeps some value fixed and shows evolution over other parameter. To explain this lets say we use in our system two variables. We can increase the (1) processors or the (2) size of our problem. In Strong scaling we set one with a fixed value. The we see the evolution of the other one. This kind of diagrams let us now how the system behaves with a big peak of load. Lets set this analogy: I have a computing device and a program running that is taking more load each time. The strong scaling shows me a “vertical scalability”, i.e. whether my system is going to function without getting blocked. In this kind of diagrams we are looking forward to see a ascending line, which tell us that the system can manage the load. 2. 2.

Weak scaling Analogous to the previous one, we may vary both of the parameters. In such a way that the x-axis stays constant to that increasement: the ratio of the two increased variables stays constant. We can see this as an “horizontal scaling” where we can see how our system behaves against many parallel tasks. We can see this with an easy example. This would be the scalability of a “modern cloud” when thousand of millions of clients access at the same time. The previous scaling would be the behaviour of an individual computer system agains a very high load. When we look at this type of diagrams we expect to see a horizontal line that shows that our system is capable to manage the load and stay even. 3. 3.

In the time breakdown we will chunk into pieces the our application and visualize through bar-plots the times of each region. This will give insights of the bottlenecks of the application.

The next part of the chapter will be an individual exam of each of the figures. The analysis will be done afterwards, so at this point just some clarifications will be given. The figures listed below are (A) 3 Strong scaling analysis (B) two weak scaling figures (C) two time breakdown diagrams.

1. A.

Strong scaling diagrams

1. (a)

Figure 1 (a) shows the behaviour of the application for each one of the three scenarios on setting a high scale factor fixed and observing the evolution at each number of processors. In this diagram we prefer a high stepped ascending line. 2. (b)

Figure 1 (b) is analogous to the previous case. The Baseline is not shown as we currently don’t have data of times at these time (we do have of TEPS). The measured parameter here is the evolution of time at different processor-number scenarios. In this case we would prefer a stepped line again. However, because we are measuring time we want this to be descending. 3. (c)

Figure 1 (c) shows an different picture of what is happening. We set the number of processor fixed (to one of the highest our application can reach). Then we observe the results for different problem sizes. Again, we look forward a stepped ascending line. 2. B.

Weak scaling diagrams

1. (a)

Figure 2 (a) and 2 (b) show both, the variation of the number of processors and the variation of the problem size in the same axis (X). The scaling is used in logarithmic scale. We do not use a constant number of vertices or edges per increased number of processors. This is due to limitations in our experiments and the fact that since we are using a 2D partitioning of the graph, the required number of processors grows exponentially. Even-though we do not keep the ratio constant we approximate this thanks to the logarithmic increasement of the scale. Regarding the diagrams we measure both Time and TEPS. Again, we would prefer a plain horizontal line as result. As it has been previously noted we are not keeping our ratio constant, so son ascending skewness will be normal. 3. C.

Time breakdown diagrams

1. (a)

Figure 3 (a) and 3 (b) will let us see one of the main goals of this work: the effect compression directly over each zone of the application

On the next section we will analyse these results and prepare the Conclusion in the next chapter.

1 Scalability analysis

To get conclusions about scalability we are first choosing the Compression scenario and the No compression one. We will be leaving the Baseline for the next subsection. Regarding the diagrams, we will be looking at the Strong and Weak scalability diagrams (the time breakdown will be used deduct conclusions about compression).

As we explained at the beginning of this chapter in the strong types of diagrams we will prefer stepped ascending lines (for TEPS) and descending (for Times. The more degree in the line, the best. Looking at Figures 1 (a) and 1 (b) we see the stepped lines. The scenarios shows a better capability to perform high in the case of compression when compared to no compression. This, regardless of being the main BFS algorithm and the state-of-the-art improvements included on this, the one that will let grow more stepped in the diagram.

We are going to start noticing an important behaviour of our application, from now own so it is going to be now when it is going to be described. We will notice that some odd scales (and also number of processors. The explanation follows) produce a very stepped low peak (in the case of TEPS) or a ascending peak in the case of Times. The reason for that lays on the management of odd units (processor ranks or matrixes with odd number of elements). We are partitioning our symmetric NxN matrix on a 2D fashion. For doing this some ranks of the initial MPI_COMM_WORLD will be designed to columns and some others to rows. Then, and from now on the communication happens using these rank partitions. The problem arises when we have to divide by 2 the rank and the number is odd. Some nodes the number of ranks on each node will be different. In our application we make use of a custom MPI_Allreduce algorithm [48] which makes use of Point to Point MPI calls for communication. When that implementation is taken to our application and we have to deal with the odd ranks, we have to create the concept of “residuum” to deal with this. Within our application we make use of residuums in two points: (i) in the column communication (ii) in the predecessor list reduction step which occurs at the end of each BFS iteration implementations. The residuum approach has two mayor downsides. In a similar way matrixes where Nx N is an odd number adds complexity when dividing the Matrix in Row and columns (the matrix are sparse-binary and operated through bitwise operations)

Complexity. The residuum is performed through the addition of many extra loops, branches, conditions, MPI difference cases, special cases when dealing with the bitmaps, more code in the Matrix structure class to deal with this problem, and so on and so forth. All this sums up and makes our implementation very complex and prone to crashes. In fact, we know that some scales (which turn out to be odd) will not run without crashes. (scales 17 and 25). 2. 2.

Low performance The two places where these residuums are located are the column and predecessor reduction steps. In the case of the latter, due to the mechanism to deal with the residuum we can not add compression. As a result the transmitted data in this phase, which takes time only once during the BFS iteration, sums more transmitted data (11 times in for $2^{16}$ Edges) than the whole transmitted (compressed) during the rest of the iteration.

We see a good scalability Strong and Weak fashions for the tested data, using compression.

2 Overhead of the added compression

To measure the added overhead by integrating compression 111our BFS implementation is a proof-of-the-concept implementation to test compression using 3rd party compression libraries and codecs., we performed tests on a third system.

Briefly this is a multiGPU platform with a NUMA architecture: two processors. Each of the processors has attached 8x NVIDIA K80 GPUs. For this experiment we expect that the GPU-BFS kernel uses direct card-to-card communication using the network device (RDMA), avoiding that way the cpu-memory path. As our compression is allocated in the CPU execution layer (we even use the CPU vectorization to introduce data-parallelism), we are going to be able to separate the execution of the BFS (GPU+GRAM+PCIe3+Network device+RDMA) with the rest of the computation, buffer management, BFS queues operated with the CPU, and lastly the compression calls. This latter operations will make use of the 2 CPUS, 2xRAM (one allocated to each processor), and a QPI high speed bus as interconnection.

Because the speed of this setup is very high and the latencies of its “network” (will be mostly bus communication: PCIe 3 and QPI) leave not much space to gain speed by compressing data, we are dealing with the perfect scenario to test how much latency are we adding. First of all to understand the introduced overhead, they are listed and explained in further detail below.

the obvious computational cost of the compression (which is pretty low in the chosen compression scheme, see section 4) 2. 2.

the data conversion before and after compression to adapt it to the requirements of the 3rd party compression libraries. This data conversions have a high memory footprint cost which may be greater than the computational cost of the compression. 3. 3.

NUMA architectures incur in an extra cost as their processors may access banks of memory allocated in other processor. The cost of these accesses is expensive [60]. This add up to the expensive memory footprints from the previous buffer copies. 4. 4.

Due to the fact of this being a compression codecs test-bench, the compression libraries have been integrated modularily. Unfortunately, many times “easy to use” equals low performance. The point being that the “modularity” has been implemented in C++ using a kind of hierarchy called (virtual keyword) which allows to choose the inherited class of other class at runtime. This prevents the compiler from knowing what class will be used, and to inline it. Each compression call using the Lemire et al. library incurs in two virtual call (the virtual keyword is also used in Lemire et al.).

As expected the overhead has been detected in the previous implementation. Also, we have noticed performance loss in other situation: In communication between MPI ranks using shared memory. The performance which compressions in these cases is lower to the one between ranks allocated in different machines (rank communicating through RDMA do not enter here)

The problem could be solved easily without many changes with the use of the implemented compression threshold (Section 4). At the beginning of the application it could be implemented a quick MPI_COMM_WORLD ranks test to detect low compression performance between ranks. We give more details aboit this idea in the Future work section (Section 8)

3 Instruction overhead analysis

In this section, we will compare our Baseline implementation with the transformed application (leaving off the compression capabilities).

We have focused the optimizations on Vectorization. We intend to make use of the Vector Processor Unit (VPU) included in the CPU of modern systems. As an example, a x86 system with SSE™ a 128-bit width VPU, would be able to operate with 4x 32-bit integers at a time. The performance boost is 4X.

Other aspects we have focused are the efficient memory access, and thread parallelism.

To see the successfulness of our optimizations we will be comparing also the scaling diagrams. This time we will be focusing on the Baseline and the No-compression versions.

From what it can be seen in the results the No-compress versions outperforms the Baseline in all scenarios (low and high scales , few and many processors). Even thought these results show success in these changes, it is difficult to assert what changes produced the better (or even some) benefit.

4 Compression analysis

Compression is the main aim of this work and for better gathering conclusions from the results we will be using other sources of data, apart from the diagrams generated in this experiments. The diagrams responsible from giving us information about the successfulness (or not) of the compressed data movement are the time breakdowns (Figures 3 (a), 3 (b)). The zones where the compression could be applied are three: (i) row communication, (ii) column communication and (iii) predecessor list reduction.

Zones (i), and (ii) have been made capable of compressing their data. For the case of (iii) several reasons have made it difficult to have it compressed in this work. For that matter, we present an algorithm for performing that step with compressed data transfers.

The efficiency of the compression has also been measured using a more accurate procedure in the cluster Creek, in University of Heidelberg. The downside, perhaps of using this cluster is that its only 16 GPUs do not allow profiling with big scale transfers.

Next follow the tables generated from the instrumentation results using the profiler.

Chapter 7 Conclusions

The conclusions of this work are based on the subsection Analysis in the previous section Results (Sections 1 and 2).

In this work we have reviewed (i) how does a compression based on “Binary Packing” (Section 4) improves a distributed BFS algorithm (ii) The effects of applying the traditional instruction optimizations to the same distributed algorithm.

Regarding whether or not have we successfully met our goal set in chapter 3 of improving the overall performance of the Baseline implementation it could be stated that yes. For the analyzed scenario, and using a relatively high scale and a relatively high number of processors we have significantly increased

the “horizontal” scalability: our distributed BFS distributes better the load among the nodes. This is an indirect effect of improving the communication latency. 2. 2.

the “vertical” scalability: our algorithm can deal with bigger amounts of load on an unique processor.

1. (a)

the direct reason for this are the instruction overhead optimizations (vectorization, memory access and thread parallelism). For the particular case seen on Figure 1 (a), the direct improvement is constant and stays over the 100% (relation between non compressed and baseline). 2. (b)

the indirect reason for this to happen is the previous optimization in the compression: since the bottleneck in the algorithm is the communication, each of the processors stay still a certain amount of time awaiting to receive the Frontier Queues from its neighbours. The Figure 1 (a) in previous chapter shows a 200% improvement of the compression implementation in relation to the baseline.

Chapter 8 Future work

We believe that some good steps to follow in order to develop a good Graph 500 implementation would be:

1. A

In case we want to run our BFS on GPGPU

1. (a)

Use the Graph 500 implementation SC12 (Ueno et al.) [51] for the BFS kernel. This implementations was already top 1 in the list in the year 2013. (Section 2). 2. (b)

The previous implementation already compress its compression. their used scheme named Variable-Length Quantity VLQ belongs to the Varint family (Section 4). Two characteristics of varint codecs are (i) that they do not perform well with big cardinalities, (ii) are slower and have smaller compression ration than a FOR + delta compression scheme. This latter about the FOR compression is also only valid if the data meets some criteria (as the varints codecs). The criteria for FOR codecs is ordered sequences of integers (unordered is also possible but requires an initial sorting) with low distances between them. Since the Frontier Queue in the sparse graphs that we are studying meet both criteria, this codec is ideal. 3. (c)

Our last step would be to add the NUMA-control111ULIBC - https://bitbucket.org/yuichiro_yasui/ulibc library that has given the position #1 to Fujisawa et al. to deal with the memory movements in our application (including the ones in the compression integration) 2. B

In case we want to run our BFS on a NUMA cluster

1. (a)

As, to our knowledge, there is no open state-of-the-art Graph500 implementation, we would need to develop it based on the paper of the Top 1 implementation [60]. 2. (b)

Since the current Top 1 does not perform compressed movement of data, and since we also know that the algorithm includes SpMV, CSR and 2D graph partitioning, it would be a good fit for the compression scheme presented here. 3. (c)

The last step would also be to add the NUMA control library that they feature to the memory structures in our Compression integration

As we partially described in section 6 when dealing with the problem of the compression-added overhead at some cases, a solution would be to learn about the ranks in the topology and know whether or not apply compression. A more optimum approach for the second case described in that section (nodes with shared memory, there the compression will not have any benefit, but only overhead) would be instead of (1) add an initial discovery cycle (2) broadcast results (3) manage structures. we could do the inverse way: since the 2D partitioning allows manually to create the partitioning for which nodes would go in the row column and the same for the row, we could keep the initial steps: (1) initial discovery of the topology (2) broadcast of the data, and now (3) would be: setup the application rank based on the map: more amount of communication $\longrightarrow$ different MPI rank group (rows or columns)

Other proposed work related with compression, but in a different scope, would encompass the use a Lempel Ziv fast compression algorithm to compress the visited neighborhood bitmap. This has been proved to be successful in other previous work (Huiwei-et-al [31]) and results on a great communications transfer reduction. However, as our bitmap transfer is small (as proven in instrumentation results, Section 2) the implementation-cost vs performance-benefit balance must be previously evaluated.

The usage of an Hybrid (MPI + OpenMP) model for parallelizing the collective operations using spare cores, shows an experimental speedup of 5X in other previous works in real world applications [32]. The revision of the dynamic task scheduling of our OpenMP thread model may also add an extra performance.

One feature wich differes with many others is our implementation is a custom MPI_Allreduce function that we is used to reduce the column communication. The implementation is based on the Rabensefner-et-al algothithm [48] which is dated in 2005. As it was described in the section 6, this algorithm uses some point-to-point communication calls which make neccessary the creation of an special case called residuum for some scenarios where |Rows x Columns| is odd (Section 6). The usage of a more new algorithm, with the possibility of being implemented over GPGPU may be a good improvement, as done in the work of Faraji-et-al [19].

Another possible addition to our implementation is the new acceleration framework OpenACC222http://www.openacc.org/. It has been successfully proved to work for other high computing problems [8], and the results may be interesting on Graph500333http://www.graph500.org/ .