1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2 #include <../src/mat/impls/aij/seq/aij.h>
3
MatColoringCreateLexicalWeights(MatColoring mc,PetscReal * weights)4 static PetscErrorCode MatColoringCreateLexicalWeights(MatColoring mc, PetscReal *weights)
5 {
6 PetscInt i, s, e;
7 Mat G = mc->mat;
8
9 PetscFunctionBegin;
10 PetscCall(MatGetOwnershipRange(G, &s, &e));
11 for (i = s; i < e; i++) weights[i - s] = i;
12 PetscFunctionReturn(PETSC_SUCCESS);
13 }
14
MatColoringCreateRandomWeights(MatColoring mc,PetscReal * weights)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;
21
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 }
34
MatColoringGetDegrees(Mat G,PetscInt distance,PetscInt * degrees)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;
49
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 }
107
MatColoringCreateLargestFirstWeights(MatColoring mc,PetscReal * weights)108 static 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;
115
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 }
134
MatColoringCreateSmallestLastWeights(MatColoring mc,PetscReal * weights)135 static 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;
153
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 }
308
MatColoringCreateWeights(MatColoring mc,PetscReal ** weights,PetscInt ** lperm)309 PetscErrorCode MatColoringCreateWeights(MatColoring mc, PetscReal **weights, PetscInt **lperm)
310 {
311 PetscInt i, s, e, n;
312 PetscReal *wts;
313
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 }
347
MatColoringSetWeights(MatColoring mc,PetscReal * weights,PetscInt * lperm)348 PetscErrorCode MatColoringSetWeights(MatColoring mc, PetscReal *weights, PetscInt *lperm)
349 {
350 PetscInt i, s, e, n;
351
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 }
376