1 2 #include <petsc/private/petscimpl.h> 3 #include <petscviewer.h> 4 5 typedef struct { 6 PetscInt id; 7 PetscInt value; 8 } HeapNode; 9 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 }; 16 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 */ 24 25 #define B 1 /* log2(ARITY) */ 26 #define ARITY (1 << B) /* tree branching factor */ 27 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) 35 36 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 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 } 62 63 PetscErrorCode PetscHeapCreate(PetscInt maxsize,PetscHeap *heap) 64 { 65 PetscErrorCode ierr; 66 PetscHeap h; 67 68 PetscFunctionBegin; 69 *heap = NULL; 70 ierr = PetscMalloc1(1,&h);CHKERRQ(ierr); 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 ierr = PetscCalloc1(h->alloc,&h->base);CHKERRQ(ierr); 75 h->base[0].id = -1; 76 h->base[0].value = PETSC_MIN_INT; 77 *heap = h; 78 PetscFunctionReturn(0); 79 } 80 81 PetscErrorCode PetscHeapAdd(PetscHeap h,PetscInt id,PetscInt val) 82 { 83 PetscInt loc,par; 84 85 PetscFunctionBegin; 86 if (1 < h->end && h->end < ARITY) h->end = ARITY; 87 loc = h->end++; 88 PetscCheckFalse(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)); 89 h->base[loc].id = id; 90 h->base[loc].value = val; 91 92 /* move up until heap condition is satisfied */ 93 while ((void)(par = Parent(loc)), Value(h,par) > val) { 94 Swap(h,loc,par); 95 loc = par; 96 } 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode PetscHeapPop(PetscHeap h,PetscInt *id,PetscInt *val) 101 { 102 PetscInt loc,chld; 103 104 PetscFunctionBegin; 105 if (h->end == 1) { 106 *id = h->base[0].id; 107 *val = h->base[0].value; 108 PetscFunctionReturn(0); 109 } 110 111 *id = h->base[1].id; 112 *val = h->base[1].value; 113 114 /* rotate last entry into first position */ 115 loc = --h->end; 116 if (h->end == ARITY) h->end = 2; /* Skip over hole */ 117 h->base[1].id = Id(h,loc); 118 h->base[1].value = Value(h,loc); 119 120 /* move down until min heap condition is satisfied */ 121 loc = 1; 122 while ((chld = MinChild(h,loc)) && Value(h,loc) > Value(h,chld)) { 123 Swap(h,loc,chld); 124 loc = chld; 125 } 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode PetscHeapPeek(PetscHeap h,PetscInt *id,PetscInt *val) 130 { 131 PetscFunctionBegin; 132 if (h->end == 1) { 133 *id = h->base[0].id; 134 *val = h->base[0].value; 135 PetscFunctionReturn(0); 136 } 137 138 *id = h->base[1].id; 139 *val = h->base[1].value; 140 PetscFunctionReturn(0); 141 } 142 143 PetscErrorCode PetscHeapStash(PetscHeap h,PetscInt id,PetscInt val) 144 { 145 PetscInt loc; 146 147 PetscFunctionBegin; 148 loc = --h->stash; 149 h->base[loc].id = id; 150 h->base[loc].value = val; 151 PetscFunctionReturn(0); 152 } 153 154 PetscErrorCode PetscHeapUnstash(PetscHeap h) 155 { 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 while (h->stash < h->alloc) { 160 PetscInt id = Id(h,h->stash),value = Value(h,h->stash); 161 h->stash++; 162 ierr = PetscHeapAdd(h,id,value);CHKERRQ(ierr); 163 } 164 PetscFunctionReturn(0); 165 } 166 167 PetscErrorCode PetscHeapDestroy(PetscHeap *heap) 168 { 169 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 ierr = PetscFree((*heap)->base);CHKERRQ(ierr); 173 ierr = PetscFree(*heap);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176 177 PetscErrorCode PetscHeapView(PetscHeap h,PetscViewer viewer) 178 { 179 PetscErrorCode ierr; 180 PetscBool iascii; 181 182 PetscFunctionBegin; 183 if (!viewer) { 184 ierr = PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);CHKERRQ(ierr); 185 } 186 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 187 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 188 if (iascii) { 189 ierr = PetscViewerASCIIPrintf(viewer,"Heap size %" PetscInt_FMT " with %" PetscInt_FMT " stashed\n",h->end-1,h->alloc-h->stash);CHKERRQ(ierr); 190 ierr = PetscViewerASCIIPrintf(viewer,"Heap in (id,value) pairs\n");CHKERRQ(ierr); 191 ierr = PetscIntView(2*(h->end-1),(const PetscInt*)(h->base+1),viewer);CHKERRQ(ierr); 192 ierr = PetscViewerASCIIPrintf(viewer,"Stash in (id,value) pairs\n");CHKERRQ(ierr); 193 ierr = PetscIntView(2*(h->alloc-h->stash),(const PetscInt*)(h->base+h->stash),viewer);CHKERRQ(ierr); 194 } 195 PetscFunctionReturn(0); 196 } 197