Actual source code: weights.c
1: #include <petsc/private/matimpl.h>
2: #include <../src/mat/impls/aij/seq/aij.h>
4: static PetscErrorCode MatColoringCreateLexicalWeights(MatColoring mc, PetscReal *weights)
5: {
6: PetscInt i, s, e;
7: Mat G = mc->mat;
9: PetscFunctionBegin;
10: PetscCall(MatGetOwnershipRange(G, &s, &e));
11: for (i = s; i < e; i++) weights[i - s] = i;
12: PetscFunctionReturn(PETSC_SUCCESS);
13: }
15: static PetscErrorCode MatColoringCreateRandomWeights(MatColoring mc, PetscReal *weights)
16: {
17: PetscInt i, s, e;
18: PetscRandom rand;
19: PetscReal r;
20: Mat G = mc->mat;
22: PetscFunctionBegin;
23: /* each weight should be the degree plus a random perturbation */
24: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand));
25: PetscCall(PetscRandomSetFromOptions(rand));
26: PetscCall(MatGetOwnershipRange(G, &s, &e));
27: for (i = s; i < e; i++) {
28: PetscCall(PetscRandomGetValueReal(rand, &r));
29: weights[i - s] = PetscAbsReal(r);
30: }
31: PetscCall(PetscRandomDestroy(&rand));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: /*@
36: MatColoringGetDegrees - Compute the number of vertices reachable in the graph within a given distance for each locally owned row.
38: Collective
40: Input Parameters:
41: + G - the graph matrix
42: - distance - the distance (in graph edges) used to define the neighborhood
44: Output Parameter:
45: . degrees - array (of local size) of neighborhood sizes for each row
47: Level: developer
49: .seealso: `MatColoring`, `MatColoringCreateWeights()`, `MatColoringSetWeights()`
50: @*/
51: PetscErrorCode MatColoringGetDegrees(Mat G, PetscInt distance, PetscInt *degrees)
52: {
53: PetscInt j, i, s, e, n, ln, lm, degree, bidx, idx, dist;
54: Mat lG, *lGs;
55: IS ris;
56: PetscInt *seen;
57: const PetscInt *gidx;
58: PetscInt *idxbuf;
59: PetscInt *distbuf;
60: PetscInt ncols;
61: const PetscInt *cols;
62: PetscBool isSEQAIJ;
63: Mat_SeqAIJ *aij;
64: PetscInt *Gi, *Gj;
66: PetscFunctionBegin;
67: PetscCall(MatGetOwnershipRange(G, &s, &e));
68: n = e - s;
69: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)G), n, s, 1, &ris));
70: PetscCall(MatIncreaseOverlap(G, 1, &ris, distance));
71: PetscCall(ISSort(ris));
72: PetscCall(MatCreateSubMatrices(G, 1, &ris, &ris, MAT_INITIAL_MATRIX, &lGs));
73: lG = lGs[0];
74: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lG, MATSEQAIJ, &isSEQAIJ));
75: PetscCheck(isSEQAIJ, PetscObjectComm((PetscObject)G), PETSC_ERR_SUP, "Requires an MPI/SEQAIJ Matrix");
76: PetscCall(MatGetSize(lG, &ln, &lm));
77: aij = (Mat_SeqAIJ *)lG->data;
78: Gi = aij->i;
79: Gj = aij->j;
80: PetscCall(PetscMalloc3(lm, &seen, lm, &idxbuf, lm, &distbuf));
81: for (i = 0; i < ln; i++) seen[i] = -1;
82: PetscCall(ISGetIndices(ris, &gidx));
83: for (i = 0; i < ln; i++) {
84: if (gidx[i] >= e || gidx[i] < s) continue;
85: bidx = -1;
86: ncols = Gi[i + 1] - Gi[i];
87: cols = &(Gj[Gi[i]]);
88: degree = 0;
89: /* place the distance-one neighbors on the queue */
90: for (j = 0; j < ncols; j++) {
91: bidx++;
92: seen[cols[j]] = i;
93: distbuf[bidx] = 1;
94: idxbuf[bidx] = cols[j];
95: }
96: while (bidx >= 0) {
97: /* pop */
98: idx = idxbuf[bidx];
99: dist = distbuf[bidx];
100: bidx--;
101: degree++;
102: if (dist < distance) {
103: ncols = Gi[idx + 1] - Gi[idx];
104: cols = &(Gj[Gi[idx]]);
105: for (j = 0; j < ncols; j++) {
106: if (seen[cols[j]] != i) {
107: bidx++;
108: seen[cols[j]] = i;
109: idxbuf[bidx] = cols[j];
110: distbuf[bidx] = dist + 1;
111: }
112: }
113: }
114: }
115: degrees[gidx[i] - s] = degree;
116: }
117: PetscCall(ISRestoreIndices(ris, &gidx));
118: PetscCall(ISDestroy(&ris));
119: PetscCall(PetscFree3(seen, idxbuf, distbuf));
120: PetscCall(MatDestroyMatrices(1, &lGs));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: static PetscErrorCode MatColoringCreateLargestFirstWeights(MatColoring mc, PetscReal *weights)
125: {
126: PetscInt i, s, e, n, ncols;
127: PetscRandom rand;
128: PetscReal r;
129: PetscInt *degrees;
130: Mat G = mc->mat;
132: PetscFunctionBegin;
133: /* each weight should be the degree plus a random perturbation */
134: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand));
135: PetscCall(PetscRandomSetFromOptions(rand));
136: PetscCall(MatGetOwnershipRange(G, &s, &e));
137: n = e - s;
138: PetscCall(PetscMalloc1(n, °rees));
139: PetscCall(MatColoringGetDegrees(G, mc->dist, degrees));
140: for (i = s; i < e; i++) {
141: PetscCall(MatGetRow(G, i, &ncols, NULL, NULL));
142: PetscCall(PetscRandomGetValueReal(rand, &r));
143: weights[i - s] = ncols + PetscAbsReal(r);
144: PetscCall(MatRestoreRow(G, i, &ncols, NULL, NULL));
145: }
146: PetscCall(PetscRandomDestroy(&rand));
147: PetscCall(PetscFree(degrees));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: static PetscErrorCode MatColoringCreateSmallestLastWeights(MatColoring mc, PetscReal *weights)
152: {
153: PetscInt *degrees, *degb, *llprev, *llnext;
154: PetscInt j, i, s, e, n, ln, lm, degree, maxdegree = 0, bidx, idx, dist, distance = mc->dist;
155: Mat lG, *lGs;
156: IS ris;
157: PetscInt *seen;
158: const PetscInt *gidx;
159: PetscInt *idxbuf;
160: PetscInt *distbuf;
161: PetscInt ncols, nxt, prv, cur;
162: const PetscInt *cols;
163: PetscBool isSEQAIJ;
164: Mat_SeqAIJ *aij;
165: PetscInt *Gi, *Gj, *rperm;
166: Mat G = mc->mat;
167: PetscReal *lweights, r;
168: PetscRandom rand;
170: PetscFunctionBegin;
171: PetscCall(MatGetOwnershipRange(G, &s, &e));
172: n = e - s;
173: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)G), n, s, 1, &ris));
174: PetscCall(MatIncreaseOverlap(G, 1, &ris, distance + 1));
175: PetscCall(ISSort(ris));
176: PetscCall(MatCreateSubMatrices(G, 1, &ris, &ris, MAT_INITIAL_MATRIX, &lGs));
177: lG = lGs[0];
178: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lG, MATSEQAIJ, &isSEQAIJ));
179: PetscCheck(isSEQAIJ, PetscObjectComm((PetscObject)G), PETSC_ERR_ARG_WRONGSTATE, "Requires an MPI/SEQAIJ Matrix");
180: PetscCall(MatGetSize(lG, &ln, &lm));
181: aij = (Mat_SeqAIJ *)lG->data;
182: Gi = aij->i;
183: Gj = aij->j;
184: PetscCall(PetscMalloc3(lm, &seen, lm, &idxbuf, lm, &distbuf));
185: PetscCall(PetscMalloc1(lm, °rees));
186: PetscCall(PetscMalloc1(lm, &lweights));
187: for (i = 0; i < ln; i++) {
188: seen[i] = -1;
189: lweights[i] = 1.;
190: }
191: PetscCall(ISGetIndices(ris, &gidx));
192: for (i = 0; i < ln; i++) {
193: bidx = -1;
194: ncols = Gi[i + 1] - Gi[i];
195: cols = &(Gj[Gi[i]]);
196: degree = 0;
197: /* place the distance-one neighbors on the queue */
198: for (j = 0; j < ncols; j++) {
199: bidx++;
200: seen[cols[j]] = i;
201: distbuf[bidx] = 1;
202: idxbuf[bidx] = cols[j];
203: }
204: while (bidx >= 0) {
205: /* pop */
206: idx = idxbuf[bidx];
207: dist = distbuf[bidx];
208: bidx--;
209: degree++;
210: if (dist < distance) {
211: ncols = Gi[idx + 1] - Gi[idx];
212: cols = &(Gj[Gi[idx]]);
213: for (j = 0; j < ncols; j++) {
214: if (seen[cols[j]] != i) {
215: bidx++;
216: seen[cols[j]] = i;
217: idxbuf[bidx] = cols[j];
218: distbuf[bidx] = dist + 1;
219: }
220: }
221: }
222: }
223: degrees[i] = degree;
224: if (degree > maxdegree) maxdegree = degree;
225: }
226: /* bucket by degree by some random permutation */
227: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand));
228: PetscCall(PetscRandomSetFromOptions(rand));
229: PetscCall(PetscMalloc1(ln, &rperm));
230: for (i = 0; i < ln; i++) {
231: PetscCall(PetscRandomGetValueReal(rand, &r));
232: lweights[i] = r;
233: rperm[i] = i;
234: }
235: PetscCall(PetscSortRealWithPermutation(lm, lweights, rperm));
236: PetscCall(PetscMalloc1(maxdegree + 1, °b));
237: PetscCall(PetscMalloc2(ln, &llnext, ln, &llprev));
238: for (i = 0; i < maxdegree + 1; i++) degb[i] = -1;
239: for (i = 0; i < ln; i++) {
240: llnext[i] = -1;
241: llprev[i] = -1;
242: seen[i] = -1;
243: }
244: for (i = 0; i < ln; i++) {
245: idx = rperm[i];
246: llnext[idx] = degb[degrees[idx]];
247: if (degb[degrees[idx]] > 0) llprev[degb[degrees[idx]]] = idx;
248: degb[degrees[idx]] = idx;
249: }
250: PetscCall(PetscFree(rperm));
251: /* remove the lowest degree one */
252: i = 0;
253: while (i != maxdegree + 1) {
254: for (i = 1; i < maxdegree + 1; i++) {
255: if (degb[i] > 0) {
256: cur = degb[i];
257: degrees[cur] = 0;
258: degb[i] = llnext[cur];
259: bidx = -1;
260: ncols = Gi[cur + 1] - Gi[cur];
261: cols = &(Gj[Gi[cur]]);
262: /* place the distance-one neighbors on the queue */
263: for (j = 0; j < ncols; j++) {
264: if (cols[j] != cur) {
265: bidx++;
266: seen[cols[j]] = i;
267: distbuf[bidx] = 1;
268: idxbuf[bidx] = cols[j];
269: }
270: }
271: while (bidx >= 0) {
272: /* pop */
273: idx = idxbuf[bidx];
274: dist = distbuf[bidx];
275: bidx--;
276: nxt = llnext[idx];
277: prv = llprev[idx];
278: if (degrees[idx] > 0) {
279: /* change up the degree of the neighbors still in the graph */
280: if (lweights[idx] <= lweights[cur]) lweights[idx] = lweights[cur] + 1;
281: if (nxt > 0) llprev[nxt] = prv;
282: if (prv > 0) {
283: llnext[prv] = nxt;
284: } else {
285: degb[degrees[idx]] = nxt;
286: }
287: degrees[idx]--;
288: llnext[idx] = degb[degrees[idx]];
289: llprev[idx] = -1;
290: if (degb[degrees[idx]] >= 0) llprev[degb[degrees[idx]]] = idx;
291: degb[degrees[idx]] = idx;
292: if (dist < distance) {
293: ncols = Gi[idx + 1] - Gi[idx];
294: cols = &(Gj[Gi[idx]]);
295: for (j = 0; j < ncols; j++) {
296: if (seen[cols[j]] != i) {
297: bidx++;
298: seen[cols[j]] = i;
299: idxbuf[bidx] = cols[j];
300: distbuf[bidx] = dist + 1;
301: }
302: }
303: }
304: }
305: }
306: break;
307: }
308: }
309: }
310: for (i = 0; i < lm; i++) {
311: if (gidx[i] >= s && gidx[i] < e) weights[gidx[i] - s] = lweights[i];
312: }
313: PetscCall(PetscRandomDestroy(&rand));
314: PetscCall(PetscFree(degb));
315: PetscCall(PetscFree2(llnext, llprev));
316: PetscCall(PetscFree(degrees));
317: PetscCall(PetscFree(lweights));
318: PetscCall(ISRestoreIndices(ris, &gidx));
319: PetscCall(ISDestroy(&ris));
320: PetscCall(PetscFree3(seen, idxbuf, distbuf));
321: PetscCall(MatDestroyMatrices(1, &lGs));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: MatColoringCreateWeights - Create per-row weights and, optionally, a decreasing-order permutation of them
327: for a `MatColoring`, using the weight scheme set on the `MatColoring`.
329: Collective
331: Input Parameter:
332: . mc - the `MatColoring`
334: Output Parameters:
335: + weights - the array of per-row weights, or `NULL` if not needed
336: - lperm - the permutation that sorts `weights` in decreasing order, or `NULL` if not needed
338: Level: developer
340: Note:
341: Both output arrays are allocated by this routine and must be freed by the caller with `PetscFree()`.
343: .seealso: `MatColoring`, `MatColoringWeightType`, `MatColoringSetWeights()`, `MatColoringGetDegrees()`
344: @*/
345: PetscErrorCode MatColoringCreateWeights(MatColoring mc, PetscReal **weights, PetscInt **lperm)
346: {
347: PetscInt i, s, e, n;
348: PetscReal *wts;
350: PetscFunctionBegin;
351: /* create weights of the specified type */
352: PetscCall(MatGetOwnershipRange(mc->mat, &s, &e));
353: n = e - s;
354: PetscCall(PetscMalloc1(n, &wts));
355: switch (mc->weight_type) {
356: case MAT_COLORING_WEIGHT_RANDOM:
357: PetscCall(MatColoringCreateRandomWeights(mc, wts));
358: break;
359: case MAT_COLORING_WEIGHT_LEXICAL:
360: PetscCall(MatColoringCreateLexicalWeights(mc, wts));
361: break;
362: case MAT_COLORING_WEIGHT_LF:
363: PetscCall(MatColoringCreateLargestFirstWeights(mc, wts));
364: break;
365: case MAT_COLORING_WEIGHT_SL:
366: PetscCall(MatColoringCreateSmallestLastWeights(mc, wts));
367: break;
368: }
369: if (lperm) {
370: PetscCall(PetscMalloc1(n, lperm));
371: for (i = 0; i < n; i++) (*lperm)[i] = i;
372: PetscCall(PetscSortRealWithPermutation(n, wts, *lperm));
373: for (i = 0; i < n / 2; i++) {
374: PetscInt swp;
375: swp = (*lperm)[i];
376: (*lperm)[i] = (*lperm)[n - 1 - i];
377: (*lperm)[n - 1 - i] = swp;
378: }
379: }
380: if (weights) *weights = wts;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@
385: MatColoringSetWeights - Provide user weights (and optionally a permutation ordering) for a `MatColoring`.
387: Collective
389: Input Parameters:
390: + mc - the `MatColoring`
391: . weights - array of per-row weights, or `NULL` to clear any previously set weights
392: - lperm - the permutation sorting `weights` in decreasing order, or `NULL` to have it computed
394: Level: developer
396: Note:
397: The arrays are copied into the `MatColoring`; the caller retains ownership of `weights` and `lperm`.
399: .seealso: `MatColoring`, `MatColoringCreateWeights()`, `MatColoringGetDegrees()`
400: @*/
401: PetscErrorCode MatColoringSetWeights(MatColoring mc, PetscReal *weights, PetscInt *lperm)
402: {
403: PetscInt i, s, e, n;
405: PetscFunctionBegin;
406: PetscCall(MatGetOwnershipRange(mc->mat, &s, &e));
407: n = e - s;
408: if (weights) {
409: PetscCall(PetscMalloc2(n, &mc->user_weights, n, &mc->user_lperm));
410: for (i = 0; i < n; i++) mc->user_weights[i] = weights[i];
411: if (!lperm) {
412: for (i = 0; i < n; i++) mc->user_lperm[i] = i;
413: PetscCall(PetscSortRealWithPermutation(n, mc->user_weights, mc->user_lperm));
414: for (i = 0; i < n / 2; i++) {
415: PetscInt swp;
416: swp = mc->user_lperm[i];
417: mc->user_lperm[i] = mc->user_lperm[n - 1 - i];
418: mc->user_lperm[n - 1 - i] = swp;
419: }
420: } else {
421: for (i = 0; i < n; i++) mc->user_lperm[i] = lperm[i];
422: }
423: } else {
424: mc->user_weights = NULL;
425: mc->user_lperm = NULL;
426: }
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }