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