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