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