Actual source code: pheap.c

      1: #include <petsc/private/petscimpl.h>
      2: #include <petscviewer.h>

      4: typedef struct {
      5:   PetscInt id;
      6:   PetscInt value;
      7: } HeapNode;

      9: struct _PetscHeap {
     10:   PetscInt  end;   /* one past the last item */
     11:   PetscInt  alloc; /* length of array */
     12:   PetscInt  stash; /* stash grows down, this points to last item */
     13:   HeapNode *base;
     14: };

     16: /*
     17:   The arity of the heap can be changed via the parameter B below. Consider the B=2 (arity=4 case below)

     19:   [00 (sentinel); 01 (min node); 10 (unused); 11 (unused); 0100 (first child); 0101; 0110; 0111; ...]

     21:   Slots 10 and 11 are referred to as the "hole" below in the implementation.
     22: */

     24: #define B     1        /* log2(ARITY) */
     25: #define ARITY (1 << B) /* tree branching factor */
     26: static inline PetscInt Parent(PetscInt loc)
     27: {
     28:   PetscInt p = loc >> B;
     29:   if (p < ARITY) return (PetscInt)(loc != 1); /* Parent(1) is 0, otherwise fix entries ending up in the hole */
     30:   return p;
     31: }
     32: #define Value(h, loc) ((h)->base[loc].value)
     33: #define Id(h, loc)    ((h)->base[loc].id)

     35: static inline void Swap(PetscHeap h, PetscInt loc, PetscInt loc2)
     36: {
     37:   PetscInt id, val;
     38:   id                  = Id(h, loc);
     39:   val                 = Value(h, loc);
     40:   h->base[loc].id     = Id(h, loc2);
     41:   h->base[loc].value  = Value(h, loc2);
     42:   h->base[loc2].id    = id;
     43:   h->base[loc2].value = val;
     44: }
     45: static inline PetscInt MinChild(PetscHeap h, PetscInt loc)
     46: {
     47:   PetscInt min, chld, left, right;
     48:   left  = loc << B;
     49:   right = PetscMin(left + ARITY - 1, h->end - 1);
     50:   chld  = 0;
     51:   min   = PETSC_INT_MAX;
     52:   for (; left <= right; left++) {
     53:     PetscInt val = Value(h, left);
     54:     if (val < min) {
     55:       min  = val;
     56:       chld = left;
     57:     }
     58:   }
     59:   return chld;
     60: }

     62: PetscErrorCode PetscHeapCreate(PetscInt maxsize, PetscHeap *heap)
     63: {
     64:   PetscHeap h;

     66:   PetscFunctionBegin;
     67:   *heap = NULL;
     68:   PetscCall(PetscMalloc1(1, &h));
     69:   h->end   = 1;
     70:   h->alloc = maxsize + ARITY; /* We waste all but one slot (loc=1) in the first ARITY slots */
     71:   h->stash = h->alloc;
     72:   PetscCall(PetscCalloc1(h->alloc, &h->base));
     73:   h->base[0].id    = -1;
     74:   h->base[0].value = PETSC_INT_MIN;
     75:   *heap            = h;
     76:   PetscFunctionReturn(PETSC_SUCCESS);
     77: }

     79: PetscErrorCode PetscHeapAdd(PetscHeap h, PetscInt id, PetscInt val)
     80: {
     81:   PetscInt loc, par;

     83:   PetscFunctionBegin;
     84:   if (1 < h->end && h->end < ARITY) h->end = ARITY;
     85:   loc = h->end++;
     86:   PetscCheck(h->end <= h->stash, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Addition would exceed allocation %" PetscInt_FMT " (%" PetscInt_FMT " stashed)", h->alloc, h->alloc - h->stash);
     87:   h->base[loc].id    = id;
     88:   h->base[loc].value = val;

     90:   /* move up until heap condition is satisfied */
     91:   while ((void)(par = Parent(loc)), Value(h, par) > val) {
     92:     Swap(h, loc, par);
     93:     loc = par;
     94:   }
     95:   PetscFunctionReturn(PETSC_SUCCESS);
     96: }

     98: PetscErrorCode PetscHeapPop(PetscHeap h, PetscInt *id, PetscInt *val)
     99: {
    100:   PetscInt loc, chld;

    102:   PetscFunctionBegin;
    103:   if (h->end == 1) {
    104:     *id  = h->base[0].id;
    105:     *val = h->base[0].value;
    106:     PetscFunctionReturn(PETSC_SUCCESS);
    107:   }

    109:   *id  = h->base[1].id;
    110:   *val = h->base[1].value;

    112:   /* rotate last entry into first position */
    113:   loc = --h->end;
    114:   if (h->end == ARITY) h->end = 2; /* Skip over hole */
    115:   h->base[1].id    = Id(h, loc);
    116:   h->base[1].value = Value(h, loc);

    118:   /* move down until min heap condition is satisfied */
    119:   loc = 1;
    120:   while ((chld = MinChild(h, loc)) && Value(h, loc) > Value(h, chld)) {
    121:     Swap(h, loc, chld);
    122:     loc = chld;
    123:   }
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: PetscErrorCode PetscHeapPeek(PetscHeap h, PetscInt *id, PetscInt *val)
    128: {
    129:   PetscFunctionBegin;
    130:   if (h->end == 1) {
    131:     *id  = h->base[0].id;
    132:     *val = h->base[0].value;
    133:     PetscFunctionReturn(PETSC_SUCCESS);
    134:   }

    136:   *id  = h->base[1].id;
    137:   *val = h->base[1].value;
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: PetscErrorCode PetscHeapStash(PetscHeap h, PetscInt id, PetscInt val)
    142: {
    143:   PetscInt loc;

    145:   PetscFunctionBegin;
    146:   loc                = --h->stash;
    147:   h->base[loc].id    = id;
    148:   h->base[loc].value = val;
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: PetscErrorCode PetscHeapUnstash(PetscHeap h)
    153: {
    154:   PetscFunctionBegin;
    155:   while (h->stash < h->alloc) {
    156:     PetscInt id = Id(h, h->stash), value = Value(h, h->stash);
    157:     h->stash++;
    158:     PetscCall(PetscHeapAdd(h, id, value));
    159:   }
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    163: PetscErrorCode PetscHeapDestroy(PetscHeap *heap)
    164: {
    165:   PetscFunctionBegin;
    166:   PetscCall(PetscFree((*heap)->base));
    167:   PetscCall(PetscFree(*heap));
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: PetscErrorCode PetscHeapView(PetscHeap h, PetscViewer viewer)
    172: {
    173:   PetscBool iascii;

    175:   PetscFunctionBegin;
    176:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer));
    178:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    179:   if (iascii) {
    180:     PetscCall(PetscViewerASCIIPrintf(viewer, "Heap size %" PetscInt_FMT " with %" PetscInt_FMT " stashed\n", h->end - 1, h->alloc - h->stash));
    181:     PetscCall(PetscViewerASCIIPrintf(viewer, "Heap in (id,value) pairs\n"));
    182:     PetscCall(PetscIntView(2 * (h->end - 1), (const PetscInt *)(h->base + 1), viewer));
    183:     PetscCall(PetscViewerASCIIPrintf(viewer, "Stash in (id,value) pairs\n"));
    184:     PetscCall(PetscIntView(2 * (h->alloc - h->stash), (const PetscInt *)(h->base + h->stash), viewer));
    185:   }
    186:   PetscFunctionReturn(PETSC_SUCCESS);
    187: }