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