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

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

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

 17: /*
 18:  * The arity of the heap can be changed via the parameter B below. Consider the B=2 (arity=4 case below)
 19:  *
 20:  * [00 (sentinel); 01 (min node); 10 (unused); 11 (unused); 0100 (first child); 0101; 0110; 0111; ...]
 21:  *
 22:  * Slots 10 and 11 are referred to as the "hole" below in the implementation.
 23:  */

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

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

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

 69:   *heap            = NULL;
 70:   PetscMalloc1(1,&h);
 71:   h->end           = 1;
 72:   h->alloc         = maxsize+ARITY; /* We waste all but one slot (loc=1) in the first ARITY slots */
 73:   h->stash         = h->alloc;
 74:   PetscMalloc1(h->alloc,&h->base);
 75:   PetscMemzero(h->base,h->alloc*sizeof(HeapNode));
 76:   h->base[0].id    = -1;
 77:   h->base[0].value = PETSC_MIN_INT;
 78:   *heap            = h;
 79:   return(0);
 80: }

 82: PetscErrorCode PetscHeapAdd(PetscHeap h,PetscInt id,PetscInt val)
 83: {
 84:   PetscInt loc,par;

 87:   if (1 < h->end && h->end < ARITY) h->end = ARITY;
 88:   loc = h->end++;
 89:   if (h->end > h->stash) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Addition would exceed allocation %D (%D stashed)",h->alloc,(h->alloc-h->stash));
 90:   h->base[loc].id    = id;
 91:   h->base[loc].value = val;

 93:   /* move up until heap condition is satisfied */
 94:   while ((void)(par = Parent(loc)), Value(h,par) > val) {
 95:     Swap(h,loc,par);
 96:     loc = par;
 97:   }
 98:   return(0);
 99: }

101: PetscErrorCode PetscHeapPop(PetscHeap h,PetscInt *id,PetscInt *val)
102: {
103:   PetscInt loc,chld;

106:   if (h->end == 1) {
107:     *id  = h->base[0].id;
108:     *val = h->base[0].value;
109:     return(0);
110:   }

112:   *id  = h->base[1].id;
113:   *val = h->base[1].value;

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

121:   /* move down until min heap condition is satisfied */
122:   loc = 1;
123:   while ((chld = MinChild(h,loc)) && Value(h,loc) > Value(h,chld)) {
124:     Swap(h,loc,chld);
125:     loc = chld;
126:   }
127:   return(0);
128: }

130: PetscErrorCode PetscHeapPeek(PetscHeap h,PetscInt *id,PetscInt *val)
131: {
133:   if (h->end == 1) {
134:     *id  = h->base[0].id;
135:     *val = h->base[0].value;
136:     return(0);
137:   }

139:   *id  = h->base[1].id;
140:   *val = h->base[1].value;
141:   return(0);
142: }

144: PetscErrorCode PetscHeapStash(PetscHeap h,PetscInt id,PetscInt val)
145: {
146:   PetscInt loc;

149:   loc                = --h->stash;
150:   h->base[loc].id    = id;
151:   h->base[loc].value = val;
152:   return(0);
153: }

155: PetscErrorCode PetscHeapUnstash(PetscHeap h)
156: {

160:   while (h->stash < h->alloc) {
161:     PetscInt id = Id(h,h->stash),value = Value(h,h->stash);
162:     h->stash++;
163:     PetscHeapAdd(h,id,value);
164:   }
165:   return(0);
166: }

168: PetscErrorCode PetscHeapDestroy(PetscHeap *heap)
169: {

173:   PetscFree((*heap)->base);
174:   PetscFree(*heap);
175:   return(0);
176: }

178: PetscErrorCode PetscHeapView(PetscHeap h,PetscViewer viewer)
179: {
181:   PetscBool      iascii;

184:   if (!viewer) {
185:     PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);
186:   }
188:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
189:   if (iascii) {
190:     PetscViewerASCIIPrintf(viewer,"Heap size %D with %D stashed\n",h->end-1,h->alloc-h->stash);
191:     PetscViewerASCIIPrintf(viewer,"Heap in (id,value) pairs\n");
192:     PetscIntView(2*(h->end-1),(const PetscInt*)(h->base+1),viewer);
193:     PetscViewerASCIIPrintf(viewer,"Stash in (id,value) pairs\n");
194:     PetscIntView(2*(h->alloc-h->stash),(const PetscInt*)(h->base+h->stash),viewer);
195:   }
196:   return(0);
197: }