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 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) 35 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 } 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 = PetscMalloc1(h->alloc,&h->base);CHKERRQ(ierr); 75 ierr = PetscMemzero(h->base,h->alloc*sizeof(HeapNode));CHKERRQ(ierr); 76 h->base[0].id = -1; 77 h->base[0].value = PETSC_MIN_INT; 78 *heap = h; 79 PetscFunctionReturn(0); 80 } 81 82 PetscErrorCode PetscHeapAdd(PetscHeap h,PetscInt id,PetscInt val) 83 { 84 PetscInt loc,par; 85 86 PetscFunctionBegin; 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; 92 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 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode PetscHeapPop(PetscHeap h,PetscInt *id,PetscInt *val) 102 { 103 PetscInt loc,chld; 104 105 PetscFunctionBegin; 106 if (h->end == 1) { 107 *id = h->base[0].id; 108 *val = h->base[0].value; 109 PetscFunctionReturn(0); 110 } 111 112 *id = h->base[1].id; 113 *val = h->base[1].value; 114 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); 120 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 PetscFunctionReturn(0); 128 } 129 130 PetscErrorCode PetscHeapPeek(PetscHeap h,PetscInt *id,PetscInt *val) 131 { 132 PetscFunctionBegin; 133 if (h->end == 1) { 134 *id = h->base[0].id; 135 *val = h->base[0].value; 136 PetscFunctionReturn(0); 137 } 138 139 *id = h->base[1].id; 140 *val = h->base[1].value; 141 PetscFunctionReturn(0); 142 } 143 144 PetscErrorCode PetscHeapStash(PetscHeap h,PetscInt id,PetscInt val) 145 { 146 PetscInt loc; 147 148 PetscFunctionBegin; 149 loc = --h->stash; 150 h->base[loc].id = id; 151 h->base[loc].value = val; 152 PetscFunctionReturn(0); 153 } 154 155 PetscErrorCode PetscHeapUnstash(PetscHeap h) 156 { 157 PetscErrorCode ierr; 158 159 PetscFunctionBegin; 160 while (h->stash < h->alloc) { 161 PetscInt id = Id(h,h->stash),value = Value(h,h->stash); 162 h->stash++; 163 ierr = PetscHeapAdd(h,id,value);CHKERRQ(ierr); 164 } 165 PetscFunctionReturn(0); 166 } 167 168 PetscErrorCode PetscHeapDestroy(PetscHeap *heap) 169 { 170 PetscErrorCode ierr; 171 172 PetscFunctionBegin; 173 ierr = PetscFree((*heap)->base);CHKERRQ(ierr); 174 ierr = PetscFree(*heap);CHKERRQ(ierr); 175 PetscFunctionReturn(0); 176 } 177 178 PetscErrorCode PetscHeapView(PetscHeap h,PetscViewer viewer) 179 { 180 PetscErrorCode ierr; 181 PetscBool iascii; 182 183 PetscFunctionBegin; 184 if (!viewer) { 185 ierr = PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);CHKERRQ(ierr); 186 } 187 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 188 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 189 if (iascii) { 190 ierr = PetscViewerASCIIPrintf(viewer,"Heap size %D with %D stashed\n",h->end-1,h->alloc-h->stash);CHKERRQ(ierr); 191 ierr = PetscViewerASCIIPrintf(viewer,"Heap in (id,value) pairs\n");CHKERRQ(ierr); 192 ierr = PetscIntView(2*(h->end-1),(const PetscInt*)(h->base+1),viewer);CHKERRQ(ierr); 193 ierr = PetscViewerASCIIPrintf(viewer,"Stash in (id,value) pairs\n");CHKERRQ(ierr); 194 ierr = PetscIntView(2*(h->alloc-h->stash),(const PetscInt*)(h->base+h->stash),viewer);CHKERRQ(ierr); 195 } 196 PetscFunctionReturn(0); 197 } 198