1 #include <../src/mat/utils/petscheap.h> 2 3 typedef struct { 4 PetscInt id; 5 PetscInt value; 6 } HeapNode; 7 8 struct _PetscHeap { 9 PetscInt end; /* one past the last item */ 10 PetscInt alloc; /* length of array */ 11 PetscInt stash; /* stash grows down, this points to last item */ 12 HeapNode *base; 13 }; 14 15 /* 16 * The arity of the heap can be changed via the parameter B below. Consider the B=2 (arity=4 case below) 17 * 18 * [00 (sentinel); 01 (min node); 10 (unused); 11 (unused); 0100 (first child); 0101; 0110; 0111; ...] 19 * 20 * Slots 10 and 11 are referred to as the "hole" below in the implementation. 21 */ 22 23 #define B 1 /* log2(ARITY) */ 24 #define ARITY (1 << B) /* tree branching factor */ 25 PETSC_STATIC_INLINE PetscInt Parent(PetscInt loc) 26 { 27 PetscInt p = loc >> B; 28 if (p < ARITY) return (PetscInt)(loc != 1); /* Parent(1) is 0, otherwise fix entries ending up in the hole */ 29 return p; 30 } 31 #define Value(h,loc) ((h)->base[loc].value) 32 #define Id(h,loc) ((h)->base[loc].id) 33 34 PETSC_STATIC_INLINE void Swap(PetscHeap h,PetscInt loc,PetscInt loc2) 35 { 36 PetscInt id,val; 37 id = Id(h,loc); 38 val = Value(h,loc); 39 h->base[loc].id = Id(h,loc2); 40 h->base[loc].value = Value(h,loc2); 41 h->base[loc2].id = id; 42 h->base[loc2].value = val; 43 } 44 PETSC_STATIC_INLINE PetscInt MinChild(PetscHeap h,PetscInt loc) 45 { 46 PetscInt min,chld,left,right; 47 left = loc << B; 48 right = PetscMin(left+ARITY-1,h->end-1); 49 chld = 0; 50 min = PETSC_MAX_INT; 51 for (; left <= right; left++) { 52 PetscInt val = Value(h,left); 53 if (val < min) { 54 min = val; 55 chld = left; 56 } 57 } 58 return chld; 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "PetscHeapCreate" 63 PetscErrorCode PetscHeapCreate(PetscInt maxsize,PetscHeap *heap) 64 { 65 PetscErrorCode ierr; 66 PetscHeap h; 67 68 PetscFunctionBegin; 69 *heap = PETSC_NULL; 70 ierr = PetscMalloc(sizeof(*h),&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 = PetscMalloc(h->alloc*sizeof(HeapNode),&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 #undef __FUNCT__ 83 #define __FUNCT__ "PetscHeapAdd" 84 PetscErrorCode PetscHeapAdd(PetscHeap h,PetscInt id,PetscInt val) 85 { 86 PetscInt loc,par; 87 88 PetscFunctionBegin; 89 if (1 < h->end && h->end < ARITY) h->end = ARITY; 90 loc = h->end++; 91 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)); 92 h->base[loc].id = id; 93 h->base[loc].value = val; 94 95 /* move up until heap condition is satisfied */ 96 while (par = Parent(loc), Value(h,par) > val) { 97 Swap(h,loc,par); 98 loc = par; 99 } 100 PetscFunctionReturn(0); 101 } 102 103 #undef __FUNCT__ 104 #define __FUNCT__ "PetscHeapPop" 105 PetscErrorCode PetscHeapPop(PetscHeap h,PetscInt *id,PetscInt *val) 106 { 107 PetscInt loc,chld; 108 109 PetscFunctionBegin; 110 if (h->end == 1) { 111 *id = h->base[0].id; 112 *val = h->base[0].value; 113 PetscFunctionReturn(0); 114 } 115 116 *id = h->base[1].id; 117 *val = h->base[1].value; 118 119 /* rotate last entry into first position */ 120 loc = --h->end; 121 if (h->end == ARITY) h->end = 2; /* Skip over hole */ 122 h->base[1].id = Id(h,loc); 123 h->base[1].value = Value(h,loc); 124 125 /* move down until min heap condition is satisfied */ 126 loc = 1; 127 while ((chld = MinChild(h,loc)) && Value(h,loc) > Value(h,chld)) { 128 Swap(h,loc,chld); 129 loc = chld; 130 } 131 PetscFunctionReturn(0); 132 } 133 134 #undef __FUNCT__ 135 #define __FUNCT__ "PetscHeapPeek" 136 PetscErrorCode PetscHeapPeek(PetscHeap h,PetscInt *id,PetscInt *val) 137 { 138 PetscFunctionBegin; 139 if (h->end == 1) { 140 *id = h->base[0].id; 141 *val = h->base[0].value; 142 PetscFunctionReturn(0); 143 } 144 145 *id = h->base[1].id; 146 *val = h->base[1].value; 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "PetscHeapStash" 152 PetscErrorCode PetscHeapStash(PetscHeap h,PetscInt id,PetscInt val) 153 { 154 PetscInt loc; 155 156 PetscFunctionBegin; 157 loc = --h->stash; 158 h->base[loc].id = id; 159 h->base[loc].value = val; 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "PetscHeapUnstash" 165 PetscErrorCode PetscHeapUnstash(PetscHeap h) 166 { 167 PetscErrorCode ierr; 168 169 PetscFunctionBegin; 170 while (h->stash < h->alloc) { 171 PetscInt id = Id(h,h->stash),value = Value(h,h->stash); 172 h->stash++; 173 ierr = PetscHeapAdd(h,id,value);CHKERRQ(ierr); 174 } 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "PetscHeapDestroy" 180 PetscErrorCode PetscHeapDestroy(PetscHeap *heap) 181 { 182 PetscErrorCode ierr; 183 184 PetscFunctionBegin; 185 ierr = PetscFree((*heap)->base);CHKERRQ(ierr); 186 ierr = PetscFree(*heap);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 #undef __FUNCT__ 191 #define __FUNCT__ "PetscHeapView" 192 PetscErrorCode PetscHeapView(PetscHeap h,PetscViewer viewer) 193 { 194 PetscErrorCode ierr; 195 PetscBool iascii; 196 197 PetscFunctionBegin; 198 if (!viewer) { 199 ierr = PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);CHKERRQ(ierr); 200 } 201 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 202 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 203 if (iascii) { 204 ierr = PetscViewerASCIIPrintf(viewer,"Heap size %D with %D stashed\n",h->end-1,h->alloc-h->stash);CHKERRQ(ierr); 205 ierr = PetscViewerASCIIPrintf(viewer,"Heap in (id,value) pairs\n");CHKERRQ(ierr); 206 ierr = PetscIntView(2*(h->end-1),(const PetscInt*)(h->base+1),viewer);CHKERRQ(ierr); 207 ierr = PetscViewerASCIIPrintf(viewer,"Stash in (id,value) pairs\n");CHKERRQ(ierr); 208 ierr = PetscIntView(2*(h->alloc-h->stash),(const PetscInt*)(h->base+h->stash),viewer);CHKERRQ(ierr); 209 } 210 PetscFunctionReturn(0); 211 } 212