Actual source code: weights.c
1: #include <petsc/private/matimpl.h>
2: #include <../src/mat/impls/aij/seq/aij.h>
4: 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: 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: PetscErrorCode MatColoringGetDegrees(Mat G, PetscInt distance, PetscInt *degrees)
36: {
37: PetscInt j, i, s, e, n, ln, lm, degree, bidx, idx, dist;
38: Mat lG, *lGs;
39: IS ris;
40: PetscInt *seen;
41: const PetscInt *gidx;
42: PetscInt *idxbuf;
43: PetscInt *distbuf;
44: PetscInt ncols;
45: const PetscInt *cols;
46: PetscBool isSEQAIJ;
47: Mat_SeqAIJ *aij;
48: PetscInt *Gi, *Gj;
50: PetscFunctionBegin;
51: PetscCall(MatGetOwnershipRange(G, &s, &e));
52: n = e - s;
53: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)G), n, s, 1, &ris));
54: PetscCall(MatIncreaseOverlap(G, 1, &ris, distance));
55: PetscCall(ISSort(ris));
56: PetscCall(MatCreateSubMatrices(G, 1, &ris, &ris, MAT_INITIAL_MATRIX, &lGs));
57: lG = lGs[0];
58: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lG, MATSEQAIJ, &isSEQAIJ));
59: PetscCheck(isSEQAIJ, PetscObjectComm((PetscObject)G), PETSC_ERR_SUP, "Requires an MPI/SEQAIJ Matrix");
60: PetscCall(MatGetSize(lG, &ln, &lm));
61: aij = (Mat_SeqAIJ *)lG->data;
62: Gi = aij->i;
63: Gj = aij->j;
64: PetscCall(PetscMalloc3(lm, &seen, lm, &idxbuf, lm, &distbuf));
65: for (i = 0; i < ln; i++) seen[i] = -1;
66: PetscCall(ISGetIndices(ris, &gidx));
67: for (i = 0; i < ln; i++) {
68: if (gidx[i] >= e || gidx[i] < s) continue;
69: bidx = -1;
70: ncols = Gi[i + 1] - Gi[i];
71: cols = &(Gj[Gi[i]]);
72: degree = 0;
73: /* place the distance-one neighbors on the queue */
74: for (j = 0; j < ncols; j++) {
75: bidx++;
76: seen[cols[j]] = i;
77: distbuf[bidx] = 1;
78: idxbuf[bidx] = cols[j];
79: }
80: while (bidx >= 0) {
81: /* pop */
82: idx = idxbuf[bidx];
83: dist = distbuf[bidx];
84: bidx--;
85: degree++;
86: if (dist < distance) {
87: ncols = Gi[idx + 1] - Gi[idx];
88: cols = &(Gj[Gi[idx]]);
89: for (j = 0; j < ncols; j++) {
90: if (seen[cols[j]] != i) {
91: bidx++;
92: seen[cols[j]] = i;
93: idxbuf[bidx] = cols[j];
94: distbuf[bidx] = dist + 1;
95: }
96: }
97: }
98: }
99: degrees[gidx[i] - s] = degree;
100: }
101: PetscCall(ISRestoreIndices(ris, &gidx));
102: PetscCall(ISDestroy(&ris));
103: PetscCall(PetscFree3(seen, idxbuf, distbuf));
104: PetscCall(MatDestroyMatrices(1, &lGs));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: PetscErrorCode MatColoringCreateLargestFirstWeights(MatColoring mc, PetscReal *weights)
109: {
110: PetscInt i, s, e, n, ncols;
111: PetscRandom rand;
112: PetscReal r;
113: PetscInt *degrees;
114: Mat G = mc->mat;
116: PetscFunctionBegin;
117: /* each weight should be the degree plus a random perturbation */
118: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand));
119: PetscCall(PetscRandomSetFromOptions(rand));
120: PetscCall(MatGetOwnershipRange(G, &s, &e));
121: n = e - s;
122: PetscCall(PetscMalloc1(n, °rees));
123: PetscCall(MatColoringGetDegrees(G, mc->dist, degrees));
124: for (i = s; i < e; i++) {
125: PetscCall(MatGetRow(G, i, &ncols, NULL, NULL));
126: PetscCall(PetscRandomGetValueReal(rand, &r));
127: weights[i - s] = ncols + PetscAbsReal(r);
128: PetscCall(MatRestoreRow(G, i, &ncols, NULL, NULL));
129: }
130: PetscCall(PetscRandomDestroy(&rand));
131: PetscCall(PetscFree(degrees));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: PetscErrorCode MatColoringCreateSmallestLastWeights(MatColoring mc, PetscReal *weights)
136: {
137: PetscInt *degrees, *degb, *llprev, *llnext;
138: PetscInt j, i, s, e, n, ln, lm, degree, maxdegree = 0, bidx, idx, dist, distance = mc->dist;
139: Mat lG, *lGs;
140: IS ris;
141: PetscInt *seen;
142: const PetscInt *gidx;
143: PetscInt *idxbuf;
144: PetscInt *distbuf;
145: PetscInt ncols, nxt, prv, cur;
146: const PetscInt *cols;
147: PetscBool isSEQAIJ;
148: Mat_SeqAIJ *aij;
149: PetscInt *Gi, *Gj, *rperm;
150: Mat G = mc->mat;
151: PetscReal *lweights, r;
152: PetscRandom rand;
154: PetscFunctionBegin;
155: PetscCall(MatGetOwnershipRange(G, &s, &e));
156: n = e - s;
157: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)G), n, s, 1, &ris));
158: PetscCall(MatIncreaseOverlap(G, 1, &ris, distance + 1));
159: PetscCall(ISSort(ris));
160: PetscCall(MatCreateSubMatrices(G, 1, &ris, &ris, MAT_INITIAL_MATRIX, &lGs));
161: lG = lGs[0];
162: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lG, MATSEQAIJ, &isSEQAIJ));
163: PetscCheck(isSEQAIJ, PetscObjectComm((PetscObject)G), PETSC_ERR_ARG_WRONGSTATE, "Requires an MPI/SEQAIJ Matrix");
164: PetscCall(MatGetSize(lG, &ln, &lm));
165: aij = (Mat_SeqAIJ *)lG->data;
166: Gi = aij->i;
167: Gj = aij->j;
168: PetscCall(PetscMalloc3(lm, &seen, lm, &idxbuf, lm, &distbuf));
169: PetscCall(PetscMalloc1(lm, °rees));
170: PetscCall(PetscMalloc1(lm, &lweights));
171: for (i = 0; i < ln; i++) {
172: seen[i] = -1;
173: lweights[i] = 1.;
174: }
175: PetscCall(ISGetIndices(ris, &gidx));
176: for (i = 0; i < ln; i++) {
177: bidx = -1;
178: ncols = Gi[i + 1] - Gi[i];
179: cols = &(Gj[Gi[i]]);
180: degree = 0;
181: /* place the distance-one neighbors on the queue */
182: for (j = 0; j < ncols; j++) {
183: bidx++;
184: seen[cols[j]] = i;
185: distbuf[bidx] = 1;
186: idxbuf[bidx] = cols[j];
187: }
188: while (bidx >= 0) {
189: /* pop */
190: idx = idxbuf[bidx];
191: dist = distbuf[bidx];
192: bidx--;
193: degree++;
194: if (dist < distance) {
195: ncols = Gi[idx + 1] - Gi[idx];
196: cols = &(Gj[Gi[idx]]);
197: for (j = 0; j < ncols; j++) {
198: if (seen[cols[j]] != i) {
199: bidx++;
200: seen[cols[j]] = i;
201: idxbuf[bidx] = cols[j];
202: distbuf[bidx] = dist + 1;
203: }
204: }
205: }
206: }
207: degrees[i] = degree;
208: if (degree > maxdegree) maxdegree = degree;
209: }
210: /* bucket by degree by some random permutation */
211: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand));
212: PetscCall(PetscRandomSetFromOptions(rand));
213: PetscCall(PetscMalloc1(ln, &rperm));
214: for (i = 0; i < ln; i++) {
215: PetscCall(PetscRandomGetValueReal(rand, &r));
216: lweights[i] = r;
217: rperm[i] = i;
218: }
219: PetscCall(PetscSortRealWithPermutation(lm, lweights, rperm));
220: PetscCall(PetscMalloc1(maxdegree + 1, °b));
221: PetscCall(PetscMalloc2(ln, &llnext, ln, &llprev));
222: for (i = 0; i < maxdegree + 1; i++) degb[i] = -1;
223: for (i = 0; i < ln; i++) {
224: llnext[i] = -1;
225: llprev[i] = -1;
226: seen[i] = -1;
227: }
228: for (i = 0; i < ln; i++) {
229: idx = rperm[i];
230: llnext[idx] = degb[degrees[idx]];
231: if (degb[degrees[idx]] > 0) llprev[degb[degrees[idx]]] = idx;
232: degb[degrees[idx]] = idx;
233: }
234: PetscCall(PetscFree(rperm));
235: /* remove the lowest degree one */
236: i = 0;
237: while (i != maxdegree + 1) {
238: for (i = 1; i < maxdegree + 1; i++) {
239: if (degb[i] > 0) {
240: cur = degb[i];
241: degrees[cur] = 0;
242: degb[i] = llnext[cur];
243: bidx = -1;
244: ncols = Gi[cur + 1] - Gi[cur];
245: cols = &(Gj[Gi[cur]]);
246: /* place the distance-one neighbors on the queue */
247: for (j = 0; j < ncols; j++) {
248: if (cols[j] != cur) {
249: bidx++;
250: seen[cols[j]] = i;
251: distbuf[bidx] = 1;
252: idxbuf[bidx] = cols[j];
253: }
254: }
255: while (bidx >= 0) {
256: /* pop */
257: idx = idxbuf[bidx];
258: dist = distbuf[bidx];
259: bidx--;
260: nxt = llnext[idx];
261: prv = llprev[idx];
262: if (degrees[idx] > 0) {
263: /* change up the degree of the neighbors still in the graph */
264: if (lweights[idx] <= lweights[cur]) lweights[idx] = lweights[cur] + 1;
265: if (nxt > 0) llprev[nxt] = prv;
266: if (prv > 0) {
267: llnext[prv] = nxt;
268: } else {
269: degb[degrees[idx]] = nxt;
270: }
271: degrees[idx]--;
272: llnext[idx] = degb[degrees[idx]];
273: llprev[idx] = -1;
274: if (degb[degrees[idx]] >= 0) llprev[degb[degrees[idx]]] = idx;
275: degb[degrees[idx]] = idx;
276: if (dist < distance) {
277: ncols = Gi[idx + 1] - Gi[idx];
278: cols = &(Gj[Gi[idx]]);
279: for (j = 0; j < ncols; j++) {
280: if (seen[cols[j]] != i) {
281: bidx++;
282: seen[cols[j]] = i;
283: idxbuf[bidx] = cols[j];
284: distbuf[bidx] = dist + 1;
285: }
286: }
287: }
288: }
289: }
290: break;
291: }
292: }
293: }
294: for (i = 0; i < lm; i++) {
295: if (gidx[i] >= s && gidx[i] < e) weights[gidx[i] - s] = lweights[i];
296: }
297: PetscCall(PetscRandomDestroy(&rand));
298: PetscCall(PetscFree(degb));
299: PetscCall(PetscFree2(llnext, llprev));
300: PetscCall(PetscFree(degrees));
301: PetscCall(PetscFree(lweights));
302: PetscCall(ISRestoreIndices(ris, &gidx));
303: PetscCall(ISDestroy(&ris));
304: PetscCall(PetscFree3(seen, idxbuf, distbuf));
305: PetscCall(MatDestroyMatrices(1, &lGs));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: PetscErrorCode MatColoringCreateWeights(MatColoring mc, PetscReal **weights, PetscInt **lperm)
310: {
311: PetscInt i, s, e, n;
312: PetscReal *wts;
314: PetscFunctionBegin;
315: /* create weights of the specified type */
316: PetscCall(MatGetOwnershipRange(mc->mat, &s, &e));
317: n = e - s;
318: PetscCall(PetscMalloc1(n, &wts));
319: switch (mc->weight_type) {
320: case MAT_COLORING_WEIGHT_RANDOM:
321: PetscCall(MatColoringCreateRandomWeights(mc, wts));
322: break;
323: case MAT_COLORING_WEIGHT_LEXICAL:
324: PetscCall(MatColoringCreateLexicalWeights(mc, wts));
325: break;
326: case MAT_COLORING_WEIGHT_LF:
327: PetscCall(MatColoringCreateLargestFirstWeights(mc, wts));
328: break;
329: case MAT_COLORING_WEIGHT_SL:
330: PetscCall(MatColoringCreateSmallestLastWeights(mc, wts));
331: break;
332: }
333: if (lperm) {
334: PetscCall(PetscMalloc1(n, lperm));
335: for (i = 0; i < n; i++) (*lperm)[i] = i;
336: PetscCall(PetscSortRealWithPermutation(n, wts, *lperm));
337: for (i = 0; i < n / 2; i++) {
338: PetscInt swp;
339: swp = (*lperm)[i];
340: (*lperm)[i] = (*lperm)[n - 1 - i];
341: (*lperm)[n - 1 - i] = swp;
342: }
343: }
344: if (weights) *weights = wts;
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: PetscErrorCode MatColoringSetWeights(MatColoring mc, PetscReal *weights, PetscInt *lperm)
349: {
350: PetscInt i, s, e, n;
352: PetscFunctionBegin;
353: PetscCall(MatGetOwnershipRange(mc->mat, &s, &e));
354: n = e - s;
355: if (weights) {
356: PetscCall(PetscMalloc2(n, &mc->user_weights, n, &mc->user_lperm));
357: for (i = 0; i < n; i++) mc->user_weights[i] = weights[i];
358: if (!lperm) {
359: for (i = 0; i < n; i++) mc->user_lperm[i] = i;
360: PetscCall(PetscSortRealWithPermutation(n, mc->user_weights, mc->user_lperm));
361: for (i = 0; i < n / 2; i++) {
362: PetscInt swp;
363: swp = mc->user_lperm[i];
364: mc->user_lperm[i] = mc->user_lperm[n - 1 - i];
365: mc->user_lperm[n - 1 - i] = swp;
366: }
367: } else {
368: for (i = 0; i < n; i++) mc->user_lperm[i] = lperm[i];
369: }
370: } else {
371: mc->user_weights = NULL;
372: mc->user_lperm = NULL;
373: }
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }