# Implementing a Heap

November 23, 2013

In this article, we implement a [Heap data structure](http://en.wikipedia.org/wiki/Heap(datastructure)). Since this is a fundamental data structure, the conceptual details of Heaps are discussed in many data structure textbooks. I shall, therefore, not repeat them here. Nonetheless, I shall discuss the Heap implementation with an example usage in mind, so that it is more concrete and provides a way to adapt the implementation for other purposes.

## Usage in mind

I wrote this module while I was a research fellow at Cardiff University. We were trying to find ways in which we could decrease particle simulation time using GPGPU, with the hope that it could be applied to the reduction of radiotherapy dose calculation times.

Within the system that was developed, a Heap data structure named ParticleRepository maintains the particles that are being simulated. Each of these particles are encapsulated by the Particle data structure. To parallelise a simulation using GPU threads, the simulation space represented by a Cuboid is first partitioned orthogonally into multiple subcuboids of equal dimensions. Every particle inside the simulation space therefore belongs to one of these subcuboids, and this is recorded inside a Particle by the unique index subcuboid. The subcuboid index is updated every time a particle leaves the current subcuboid.

Within the particles repository heap, all of the particles are ordered using the subcuboid index. This allows the simulation driver to dispatch particle groups that are inside the same subcuboid to architecturally related group of threads (e.g., CUDA thread block), thus maximising sharing of simulation environment that are common to all particles inside a subcuboid (e.g., CSG models of machine components). The particles repository heap itself is stored in memory as a fixed array of particles, named head, which is an array implementation of a complete binary tree. In the next post, we shall consider the case where we use paging to allow the particles repository heap to grow when the particles simulated no longer fits inside the fixed array.

## Add particle to particle repository

Function heap_insert_fast(pr,t) inserts a new particle t into the particles repository pr. We bubble up the new node to its rightful place using the containing subcuboid as the comparison key. This bubbling is done by climbing up the tree starting at the array position dictated by the array implementation of a complete binary tree, and moving the node from child to parent as we climb up. The insertion is complete when the node is finally placed at its rightful place based on the outcome of the key comparison.

void heap_insert_fast(ParticleRepository *pr, Particle *t)
{
uint32_t p, n; /* indices of parent and current node */
Particle *fp;
fp = pr->head; /* complete binary tree: array of particles */
n = ++pr->count; /* start bubbling from the last node */
p = n >> 1;
while (p > 0) { /* bubble up the tree */
if (fp[p].subcuboid < t->subcuboid)
fp[n] = fp[p];  /* move from parent to child */
else break; /* rightful place found */
n = p; /* climb up the tree */
p = n >> 1;
}
if (0 == t->id)
t->id = pr->pid++; /* give unique id to particle */
fp[n] = *t; /* place node */
}

## Remove particle from particle repository

Function heap_remove_fast(pr,t) removes the particle at the top of the heap representing the particles repository pr, and stores the particle into t. We remove the first node, which is always the top of the heap. We then fill this void by promoting the last node of the complete binary tree to become the root. Finally, we re-balance the heap by bubbling down the root until we find its rightful place, or we find that it has become a leaf. While bubbling down, we move either the left or right subtree root to the current node position as we climb down the tree. This is again based on the outcome of a key comparison.

void heap_remove_fast(ParticleRepository *pr, Particle *t)
{
uint32_t n, l, r, q; /* indices to node, left, right and last */
Particle temp, *fp;
fp = pr->head; /* complete binary tree: array of particles */
*t = fp[1]; /* particle to return */
q = pr->count--;
if (q) {
temp = fp[q]; /* the last node to be promoted */
n = 1; /* start at the root */
l = 2;
/* bubble down the tree choosing left or right subtree */
while (l < q) {
r = l + 1;
if (fp[q].subcuboid < fp[l].subcuboid) {
if (r < q && fp[l].subcuboid < fp[r].subcuboid) {
fp[n] = fp[r];
n = r;
} else {
fp[n] = fp[l];
n = l;
}
} else {
if (r < q && fp[q].subcuboid < fp[r].subcuboid) {
fp[n] = fp[r];
n = r;
} else break; /* rightful place found */
}
l = n << 1; /* go deeper until we have passed a leaf */
}
fp[n] = temp; /* place node */
}
}