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