Actual source code: weights.c

petsc-3.8.4 2018-03-24
Report Typos and Errors
  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,&degrees);
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,&degrees);
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,&degb);
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: }