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