In this article, we implement a [Heap data
structure](http://en.wikipedia.org/wiki/Heap*(data*structure)). 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 */
}
}
```