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,&degrees);
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,&degrees);
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,&degb);
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: }