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 PetscInt p = loc >> B; 27 if (p < ARITY) return (PetscInt)(loc != 1); /* Parent(1) is 0, otherwise fix entries ending up in the hole */ 28 return p; 29 } 30 #define Value(h,loc) ((h)->base[loc].value) 31 #define Id(h,loc) ((h)->base[loc].id) 32 33 PETSC_STATIC_INLINE void Swap(PetscHeap h,PetscInt loc,PetscInt loc2) { 34 PetscInt id,val; 35 id = Id(h,loc); 36 val = Value(h,loc); 37 h->base[loc].id = Id(h,loc2); 38 h->base[loc].value = Value(h,loc2); 39 h->base[loc2].id = id; 40 h->base[loc2].value = val; 41 } 42 PETSC_STATIC_INLINE PetscInt MinChild(PetscHeap h,PetscInt loc) { 43 PetscInt min,chld,left,right; 44 left = loc << B; 45 right = PetscMin(left+ARITY-1,h->end-1); 46 chld = 0; 47 min = PETSC_MAX_INT; 48 for ( ; left <= right; left++) { 49 PetscInt val = Value(h,left); 50 if (val < min) { 51 min = val; 52 chld = left; 53 } 54 } 55 return chld; 56 } 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "PetscHeapCreate" 60 PetscErrorCode PetscHeapCreate(PetscInt maxsize,PetscHeap *heap) 61 { 62 PetscErrorCode ierr; 63 PetscHeap h; 64 65 PetscFunctionBegin; 66 *heap = PETSC_NULL; 67 ierr = PetscMalloc(sizeof *h,&h);CHKERRQ(ierr); 68 h->end = 1; 69 h->alloc = maxsize+ARITY; /* We waste all but one slot (loc=1) in the first ARITY slots */ 70 h->stash = h->alloc; 71 ierr = PetscMalloc(h->alloc*sizeof(HeapNode),&h->base);CHKERRQ(ierr); 72 ierr = PetscMemzero(h->base,h->alloc*sizeof(HeapNode));CHKERRQ(ierr); 73 h->base[0].id = -1; 74 h->base[0].value = PETSC_MIN_INT; 75 *heap = h; 76 PetscFunctionReturn(0); 77 } 78 79 #undef __FUNCT__ 80 #define __FUNCT__ "PetscHeapAdd" 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 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)); 89 h->base[loc].id = id; 90 h->base[loc].value = val; 91 92 /* move up until heap condition is satisfied */ 93 while (par = Parent(loc), Value(h,par) > val) { 94 Swap(h,loc,par); 95 loc = par; 96 } 97 PetscFunctionReturn(0); 98 } 99 100 #undef __FUNCT__ 101 #define __FUNCT__ "PetscHeapPop" 102 PetscErrorCode PetscHeapPop(PetscHeap h,PetscInt *id,PetscInt *val) 103 { 104 PetscInt loc,chld; 105 106 PetscFunctionBegin; 107 if (h->end == 1) { 108 *id = h->base[0].id; 109 *val = h->base[0].value; 110 PetscFunctionReturn(0); 111 } 112 113 *id = h->base[1].id; 114 *val = h->base[1].value; 115 116 /* rotate last entry into first position */ 117 loc = --h->end; 118 if (h->end == ARITY) h->end = 2; /* Skip over hole */ 119 h->base[1].id = Id(h,loc); 120 h->base[1].value = Value(h,loc); 121 122 /* move down until min heap condition is satisfied */ 123 loc = 1; 124 while ((chld = MinChild(h,loc)) && Value(h,loc) > Value(h,chld)) { 125 Swap(h,loc,chld); 126 loc = chld; 127 } 128 PetscFunctionReturn(0); 129 } 130 131 #undef __FUNCT__ 132 #define __FUNCT__ "PetscHeapPeek" 133 PetscErrorCode PetscHeapPeek(PetscHeap h,PetscInt *id,PetscInt *val) 134 { 135 PetscFunctionBegin; 136 if (h->end == 1) { 137 *id = h->base[0].id; 138 *val = h->base[0].value; 139 PetscFunctionReturn(0); 140 } 141 142 *id = h->base[1].id; 143 *val = h->base[1].value; 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "PetscHeapStash" 149 PetscErrorCode PetscHeapStash(PetscHeap h,PetscInt id,PetscInt val) 150 { 151 PetscInt loc; 152 153 PetscFunctionBegin; 154 loc = --h->stash; 155 h->base[loc].id = id; 156 h->base[loc].value = val; 157 PetscFunctionReturn(0); 158 } 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "PetscHeapUnstash" 162 PetscErrorCode PetscHeapUnstash(PetscHeap h) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 while (h->stash < h->alloc) { 168 PetscInt id = Id(h,h->stash),value = Value(h,h->stash); 169 h->stash++; 170 ierr = PetscHeapAdd(h,id,value);CHKERRQ(ierr); 171 } 172 PetscFunctionReturn(0); 173 } 174 175 #undef __FUNCT__ 176 #define __FUNCT__ "PetscHeapDestroy" 177 PetscErrorCode PetscHeapDestroy(PetscHeap *heap) 178 { 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 ierr = PetscFree((*heap)->base);CHKERRQ(ierr); 183 ierr = PetscFree(*heap);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "PetscHeapView" 189 PetscErrorCode PetscHeapView(PetscHeap h,PetscViewer viewer) 190 { 191 PetscErrorCode ierr; 192 PetscBool iascii; 193 194 PetscFunctionBegin; 195 if (!viewer) { 196 ierr = PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);CHKERRQ(ierr); 197 } 198 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 199 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 200 if (iascii) { 201 ierr = PetscViewerASCIIPrintf(viewer,"Heap size %D with %D stashed\n",h->end-1,h->alloc-h->stash);CHKERRQ(ierr); 202 ierr = PetscViewerASCIIPrintf(viewer,"Heap in (id,value) pairs\n");CHKERRQ(ierr); 203 ierr = PetscIntView(2*(h->end-1),(const PetscInt*)(h->base+1),viewer);CHKERRQ(ierr); 204 ierr = PetscViewerASCIIPrintf(viewer,"Stash in (id,value) pairs\n");CHKERRQ(ierr); 205 ierr = PetscIntView(2*(h->alloc-h->stash),(const PetscInt*)(h->base+h->stash),viewer);CHKERRQ(ierr); 206 } 207 PetscFunctionReturn(0); 208 } 209