Actual source code: pheap.c
1: #include <petsc/private/petscimpl.h>
3: typedef struct {
4: PetscInt id;
5: PetscInt value;
6: } HeapNode;
8: struct _n_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: };
15: /*
16: The arity of the heap can be changed via the parameter B below. Consider the B=2 (arity=4 case below)
18: [00 (sentinel); 01 (min node); 10 (unused); 11 (unused); 0100 (first child); 0101; 0110; 0111; ...]
20: Slots 10 and 11 are referred to as the "hole" below in the implementation.
21: */
23: #define B 1 /* log2(ARITY) */
24: #define ARITY (1 << B) /* tree branching factor */
25: 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)
34: 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: 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_INT_MAX;
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: }
61: /*@
62: PetscHeapCreate - Creates a `PetscHeap` object, a simple min-heap for `(id, value)` pairs.
64: Not Collective
66: Input Parameter:
67: . maxsize - the maximum number of items the heap can hold at once
69: Output Parameter:
70: . heap - the newly created `PetscHeap` object
72: Level: developer
74: Note:
75: The heap is ordered by `value`; items with equal values may be returned in any order.
77: .seealso: `PetscHeap`, `PetscHeapAdd()`, `PetscHeapPop()`, `PetscHeapPeek()`, `PetscHeapStash()`, `PetscHeapUnstash()`, `PetscHeapView()`, `PetscHeapDestroy()`
78: @*/
79: PetscErrorCode PetscHeapCreate(PetscInt maxsize, PetscHeap *heap)
80: {
81: PetscHeap h;
83: PetscFunctionBegin;
84: *heap = NULL;
85: PetscCall(PetscMalloc1(1, &h));
86: h->end = 1;
87: h->alloc = maxsize + ARITY; /* We waste all but one slot (loc=1) in the first ARITY slots */
88: h->stash = h->alloc;
89: PetscCall(PetscCalloc1(h->alloc, &h->base));
90: h->base[0].id = -1;
91: h->base[0].value = PETSC_INT_MIN;
92: *heap = h;
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: /*@
97: PetscHeapAdd - Insert an item into a `PetscHeap`.
99: Not Collective
101: Input Parameters:
102: + h - the `PetscHeap`
103: . id - the item identifier
104: - val - the value used for heap ordering
106: Level: developer
108: .seealso: `PetscHeap`, `PetscHeapCreate()`, `PetscHeapPop()`, `PetscHeapPeek()`, `PetscHeapStash()`, `PetscHeapUnstash()`, `PetscHeapDestroy()`
109: @*/
110: PetscErrorCode PetscHeapAdd(PetscHeap h, PetscInt id, PetscInt val)
111: {
112: PetscInt loc, par;
114: PetscFunctionBegin;
115: if (1 < h->end && h->end < ARITY) h->end = ARITY;
116: loc = h->end++;
117: PetscCheck(h->end <= h->stash, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Addition would exceed allocation %" PetscInt_FMT " (%" PetscInt_FMT " stashed)", h->alloc, h->alloc - h->stash);
118: h->base[loc].id = id;
119: h->base[loc].value = val;
121: /* move up until heap condition is satisfied */
122: while ((void)(par = Parent(loc)), Value(h, par) > val) {
123: Swap(h, loc, par);
124: loc = par;
125: }
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*@
130: PetscHeapPop - Remove and return the minimum item from a `PetscHeap`.
132: Not Collective
134: Input Parameter:
135: . h - the `PetscHeap`
137: Output Parameters:
138: + id - identifier of the popped item, or `-1` if the heap is empty
139: - val - value of the popped item, or `PETSC_INT_MIN` if the heap is empty
141: Level: developer
143: .seealso: `PetscHeap`, `PetscHeapCreate()`, `PetscHeapAdd()`, `PetscHeapPeek()`, `PetscHeapStash()`, `PetscHeapUnstash()`, `PetscHeapDestroy()`
144: @*/
145: PetscErrorCode PetscHeapPop(PetscHeap h, PetscInt *id, PetscInt *val)
146: {
147: PetscInt loc, chld;
149: PetscFunctionBegin;
150: if (h->end == 1) {
151: *id = h->base[0].id;
152: *val = h->base[0].value;
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: *id = h->base[1].id;
157: *val = h->base[1].value;
159: /* rotate last entry into first position */
160: loc = --h->end;
161: if (h->end == ARITY) h->end = 2; /* Skip over hole */
162: h->base[1].id = Id(h, loc);
163: h->base[1].value = Value(h, loc);
165: /* move down until min heap condition is satisfied */
166: loc = 1;
167: while ((chld = MinChild(h, loc)) && Value(h, loc) > Value(h, chld)) {
168: Swap(h, loc, chld);
169: loc = chld;
170: }
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@
175: PetscHeapPeek - Return the minimum item of a `PetscHeap` without removing it.
177: Not Collective
179: Input Parameter:
180: . h - the `PetscHeap`
182: Output Parameters:
183: + id - identifier of the minimum item, or `-1` if the heap is empty
184: - val - value of the minimum item, or `PETSC_INT_MIN` if the heap is empty
186: Level: developer
188: .seealso: `PetscHeap`, `PetscHeapCreate()`, `PetscHeapAdd()`, `PetscHeapPop()`, `PetscHeapDestroy()`
189: @*/
190: PetscErrorCode PetscHeapPeek(PetscHeap h, PetscInt *id, PetscInt *val)
191: {
192: PetscFunctionBegin;
193: if (h->end == 1) {
194: *id = h->base[0].id;
195: *val = h->base[0].value;
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: *id = h->base[1].id;
200: *val = h->base[1].value;
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@
205: PetscHeapStash - Set aside an item in a `PetscHeap` for later insertion via `PetscHeapUnstash()`.
207: Not Collective
209: Input Parameters:
210: + h - the `PetscHeap`
211: . id - the item identifier
212: - val - the value used for heap ordering
214: Level: developer
216: Note:
217: Stashed items are held in the trailing portion of the heap's storage and do not participate in
218: heap ordering until `PetscHeapUnstash()` reinserts them.
220: .seealso: `PetscHeap`, `PetscHeapCreate()`, `PetscHeapAdd()`, `PetscHeapUnstash()`, `PetscHeapDestroy()`
221: @*/
222: PetscErrorCode PetscHeapStash(PetscHeap h, PetscInt id, PetscInt val)
223: {
224: PetscInt loc;
226: PetscFunctionBegin;
227: loc = --h->stash;
228: h->base[loc].id = id;
229: h->base[loc].value = val;
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: /*@
234: PetscHeapUnstash - Reinsert all items previously stashed with `PetscHeapStash()` into the heap.
236: Not Collective
238: Input Parameter:
239: . h - the `PetscHeap`
241: Level: developer
243: .seealso: `PetscHeap`, `PetscHeapCreate()`, `PetscHeapAdd()`, `PetscHeapStash()`, `PetscHeapDestroy()`
244: @*/
245: PetscErrorCode PetscHeapUnstash(PetscHeap h)
246: {
247: PetscFunctionBegin;
248: while (h->stash < h->alloc) {
249: PetscInt id = Id(h, h->stash), value = Value(h, h->stash);
250: h->stash++;
251: PetscCall(PetscHeapAdd(h, id, value));
252: }
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: PetscHeapDestroy - Destroys a `PetscHeap` created with `PetscHeapCreate()`.
259: Not Collective
261: Input Parameter:
262: . heap - the `PetscHeap` to destroy; set to `NULL` on return
264: Level: developer
266: .seealso: `PetscHeap`, `PetscHeapCreate()`
267: @*/
268: PetscErrorCode PetscHeapDestroy(PetscHeap *heap)
269: {
270: PetscFunctionBegin;
271: PetscCall(PetscFree((*heap)->base));
272: PetscCall(PetscFree(*heap));
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: /*@
277: PetscHeapView - View the contents of a `PetscHeap`, including any stashed items.
279: Not Collective
281: Input Parameters:
282: + h - the `PetscHeap`
283: - viewer - a `PetscViewer`, or `NULL` to use `PETSC_VIEWER_STDOUT_SELF`
285: Level: developer
287: .seealso: `PetscHeap`, `PetscHeapCreate()`, `PetscHeapAdd()`, `PetscHeapPop()`
288: @*/
289: PetscErrorCode PetscHeapView(PetscHeap h, PetscViewer viewer)
290: {
291: PetscBool isascii;
293: PetscFunctionBegin;
294: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer));
296: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
297: if (isascii) {
298: PetscCall(PetscViewerASCIIPrintf(viewer, "Heap size %" PetscInt_FMT " with %" PetscInt_FMT " stashed\n", h->end - 1, h->alloc - h->stash));
299: PetscCall(PetscViewerASCIIPrintf(viewer, "Heap in (id,value) pairs\n"));
300: PetscCall(PetscIntView(2 * (h->end - 1), (const PetscInt *)(h->base + 1), viewer));
301: PetscCall(PetscViewerASCIIPrintf(viewer, "Stash in (id,value) pairs\n"));
302: PetscCall(PetscIntView(2 * (h->alloc - h->stash), (const PetscInt *)(h->base + h->stash), viewer));
303: }
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }