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 */
Parent(PetscInt loc)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
Swap(PetscHeap h,PetscInt loc,PetscInt loc2)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 }
MinChild(PetscHeap h,PetscInt loc)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
PetscHeapCreate(PetscInt maxsize,PetscHeap * heap)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
PetscHeapAdd(PetscHeap h,PetscInt id,PetscInt val)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 ((void)(par = Parent(loc)), Value(h,par) > val) {
94 Swap(h,loc,par);
95 loc = par;
96 }
97 PetscFunctionReturn(0);
98 }
99
PetscHeapPop(PetscHeap h,PetscInt * id,PetscInt * val)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
PetscHeapPeek(PetscHeap h,PetscInt * id,PetscInt * val)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
PetscHeapStash(PetscHeap h,PetscInt id,PetscInt val)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
PetscHeapUnstash(PetscHeap h)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
PetscHeapDestroy(PetscHeap * heap)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
PetscHeapView(PetscHeap h,PetscViewer viewer)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 %D with %D 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