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++) weights[i - s] = i;
 11:   return 0;
 12: }

 14: PetscErrorCode MatColoringCreateRandomWeights(MatColoring mc, PetscReal *weights)
 15: {
 16:   PetscInt    i, s, e;
 17:   PetscRandom rand;
 18:   PetscReal   r;
 19:   Mat         G = mc->mat;

 21:   /* each weight should be the degree plus a random perturbation */
 22:   PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand);
 23:   PetscRandomSetFromOptions(rand);
 24:   MatGetOwnershipRange(G, &s, &e);
 25:   for (i = s; i < e; i++) {
 26:     PetscRandomGetValueReal(rand, &r);
 27:     weights[i - s] = PetscAbsReal(r);
 28:   }
 29:   PetscRandomDestroy(&rand);
 30:   return 0;
 31: }

 33: PetscErrorCode MatColoringGetDegrees(Mat G, PetscInt distance, PetscInt *degrees)
 34: {
 35:   PetscInt        j, i, s, e, n, ln, lm, degree, bidx, idx, dist;
 36:   Mat             lG, *lGs;
 37:   IS              ris;
 38:   PetscInt       *seen;
 39:   const PetscInt *gidx;
 40:   PetscInt       *idxbuf;
 41:   PetscInt       *distbuf;
 42:   PetscInt        ncols;
 43:   const PetscInt *cols;
 44:   PetscBool       isSEQAIJ;
 45:   Mat_SeqAIJ     *aij;
 46:   PetscInt       *Gi, *Gj;

 48:   MatGetOwnershipRange(G, &s, &e);
 49:   n = e - s;
 50:   ISCreateStride(PetscObjectComm((PetscObject)G), n, s, 1, &ris);
 51:   MatIncreaseOverlap(G, 1, &ris, distance);
 52:   ISSort(ris);
 53:   MatCreateSubMatrices(G, 1, &ris, &ris, MAT_INITIAL_MATRIX, &lGs);
 54:   lG = lGs[0];
 55:   PetscObjectBaseTypeCompare((PetscObject)lG, MATSEQAIJ, &isSEQAIJ);
 57:   MatGetSize(lG, &ln, &lm);
 58:   aij = (Mat_SeqAIJ *)lG->data;
 59:   Gi  = aij->i;
 60:   Gj  = aij->j;
 61:   PetscMalloc3(lm, &seen, lm, &idxbuf, lm, &distbuf);
 62:   for (i = 0; i < ln; i++) seen[i] = -1;
 63:   ISGetIndices(ris, &gidx);
 64:   for (i = 0; i < ln; i++) {
 65:     if (gidx[i] >= e || gidx[i] < s) continue;
 66:     bidx   = -1;
 67:     ncols  = Gi[i + 1] - Gi[i];
 68:     cols   = &(Gj[Gi[i]]);
 69:     degree = 0;
 70:     /* place the distance-one neighbors on the queue */
 71:     for (j = 0; j < ncols; j++) {
 72:       bidx++;
 73:       seen[cols[j]] = i;
 74:       distbuf[bidx] = 1;
 75:       idxbuf[bidx]  = cols[j];
 76:     }
 77:     while (bidx >= 0) {
 78:       /* pop */
 79:       idx  = idxbuf[bidx];
 80:       dist = distbuf[bidx];
 81:       bidx--;
 82:       degree++;
 83:       if (dist < distance) {
 84:         ncols = Gi[idx + 1] - Gi[idx];
 85:         cols  = &(Gj[Gi[idx]]);
 86:         for (j = 0; j < ncols; j++) {
 87:           if (seen[cols[j]] != i) {
 88:             bidx++;
 89:             seen[cols[j]] = i;
 90:             idxbuf[bidx]  = cols[j];
 91:             distbuf[bidx] = dist + 1;
 92:           }
 93:         }
 94:       }
 95:     }
 96:     degrees[gidx[i] - s] = degree;
 97:   }
 98:   ISRestoreIndices(ris, &gidx);
 99:   ISDestroy(&ris);
100:   PetscFree3(seen, idxbuf, distbuf);
101:   MatDestroyMatrices(1, &lGs);
102:   return 0;
103: }

105: PetscErrorCode MatColoringCreateLargestFirstWeights(MatColoring mc, PetscReal *weights)
106: {
107:   PetscInt    i, s, e, n, ncols;
108:   PetscRandom rand;
109:   PetscReal   r;
110:   PetscInt   *degrees;
111:   Mat         G = mc->mat;

113:   /* each weight should be the degree plus a random perturbation */
114:   PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand);
115:   PetscRandomSetFromOptions(rand);
116:   MatGetOwnershipRange(G, &s, &e);
117:   n = e - s;
118:   PetscMalloc1(n, &degrees);
119:   MatColoringGetDegrees(G, mc->dist, degrees);
120:   for (i = s; i < e; i++) {
121:     MatGetRow(G, i, &ncols, NULL, NULL);
122:     PetscRandomGetValueReal(rand, &r);
123:     weights[i - s] = ncols + PetscAbsReal(r);
124:     MatRestoreRow(G, i, &ncols, NULL, NULL);
125:   }
126:   PetscRandomDestroy(&rand);
127:   PetscFree(degrees);
128:   return 0;
129: }

131: PetscErrorCode MatColoringCreateSmallestLastWeights(MatColoring mc, PetscReal *weights)
132: {
133:   PetscInt       *degrees, *degb, *llprev, *llnext;
134:   PetscInt        j, i, s, e, n, ln, lm, degree, maxdegree = 0, bidx, idx, dist, distance = mc->dist;
135:   Mat             lG, *lGs;
136:   IS              ris;
137:   PetscInt       *seen;
138:   const PetscInt *gidx;
139:   PetscInt       *idxbuf;
140:   PetscInt       *distbuf;
141:   PetscInt        ncols, nxt, prv, cur;
142:   const PetscInt *cols;
143:   PetscBool       isSEQAIJ;
144:   Mat_SeqAIJ     *aij;
145:   PetscInt       *Gi, *Gj, *rperm;
146:   Mat             G = mc->mat;
147:   PetscReal      *lweights, r;
148:   PetscRandom     rand;

150:   MatGetOwnershipRange(G, &s, &e);
151:   n = e - s;
152:   ISCreateStride(PetscObjectComm((PetscObject)G), n, s, 1, &ris);
153:   MatIncreaseOverlap(G, 1, &ris, distance + 1);
154:   ISSort(ris);
155:   MatCreateSubMatrices(G, 1, &ris, &ris, MAT_INITIAL_MATRIX, &lGs);
156:   lG = lGs[0];
157:   PetscObjectBaseTypeCompare((PetscObject)lG, MATSEQAIJ, &isSEQAIJ);
159:   MatGetSize(lG, &ln, &lm);
160:   aij = (Mat_SeqAIJ *)lG->data;
161:   Gi  = aij->i;
162:   Gj  = aij->j;
163:   PetscMalloc3(lm, &seen, lm, &idxbuf, lm, &distbuf);
164:   PetscMalloc1(lm, &degrees);
165:   PetscMalloc1(lm, &lweights);
166:   for (i = 0; i < ln; i++) {
167:     seen[i]     = -1;
168:     lweights[i] = 1.;
169:   }
170:   ISGetIndices(ris, &gidx);
171:   for (i = 0; i < ln; i++) {
172:     bidx   = -1;
173:     ncols  = Gi[i + 1] - Gi[i];
174:     cols   = &(Gj[Gi[i]]);
175:     degree = 0;
176:     /* place the distance-one neighbors on the queue */
177:     for (j = 0; j < ncols; j++) {
178:       bidx++;
179:       seen[cols[j]] = i;
180:       distbuf[bidx] = 1;
181:       idxbuf[bidx]  = cols[j];
182:     }
183:     while (bidx >= 0) {
184:       /* pop */
185:       idx  = idxbuf[bidx];
186:       dist = distbuf[bidx];
187:       bidx--;
188:       degree++;
189:       if (dist < distance) {
190:         ncols = Gi[idx + 1] - Gi[idx];
191:         cols  = &(Gj[Gi[idx]]);
192:         for (j = 0; j < ncols; j++) {
193:           if (seen[cols[j]] != i) {
194:             bidx++;
195:             seen[cols[j]] = i;
196:             idxbuf[bidx]  = cols[j];
197:             distbuf[bidx] = dist + 1;
198:           }
199:         }
200:       }
201:     }
202:     degrees[i] = degree;
203:     if (degree > maxdegree) maxdegree = degree;
204:   }
205:   /* bucket by degree by some random permutation */
206:   PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand);
207:   PetscRandomSetFromOptions(rand);
208:   PetscMalloc1(ln, &rperm);
209:   for (i = 0; i < ln; i++) {
210:     PetscRandomGetValueReal(rand, &r);
211:     lweights[i] = r;
212:     rperm[i]    = i;
213:   }
214:   PetscSortRealWithPermutation(lm, lweights, rperm);
215:   PetscMalloc1(maxdegree + 1, &degb);
216:   PetscMalloc2(ln, &llnext, ln, &llprev);
217:   for (i = 0; i < maxdegree + 1; i++) degb[i] = -1;
218:   for (i = 0; i < ln; i++) {
219:     llnext[i] = -1;
220:     llprev[i] = -1;
221:     seen[i]   = -1;
222:   }
223:   for (i = 0; i < ln; i++) {
224:     idx         = rperm[i];
225:     llnext[idx] = degb[degrees[idx]];
226:     if (degb[degrees[idx]] > 0) llprev[degb[degrees[idx]]] = idx;
227:     degb[degrees[idx]] = idx;
228:   }
229:   PetscFree(rperm);
230:   /* remove the lowest degree one */
231:   i = 0;
232:   while (i != maxdegree + 1) {
233:     for (i = 1; i < maxdegree + 1; i++) {
234:       if (degb[i] > 0) {
235:         cur          = degb[i];
236:         degrees[cur] = 0;
237:         degb[i]      = llnext[cur];
238:         bidx         = -1;
239:         ncols        = Gi[cur + 1] - Gi[cur];
240:         cols         = &(Gj[Gi[cur]]);
241:         /* place the distance-one neighbors on the queue */
242:         for (j = 0; j < ncols; j++) {
243:           if (cols[j] != cur) {
244:             bidx++;
245:             seen[cols[j]] = i;
246:             distbuf[bidx] = 1;
247:             idxbuf[bidx]  = cols[j];
248:           }
249:         }
250:         while (bidx >= 0) {
251:           /* pop */
252:           idx  = idxbuf[bidx];
253:           dist = distbuf[bidx];
254:           bidx--;
255:           nxt = llnext[idx];
256:           prv = llprev[idx];
257:           if (degrees[idx] > 0) {
258:             /* change up the degree of the neighbors still in the graph */
259:             if (lweights[idx] <= lweights[cur]) lweights[idx] = lweights[cur] + 1;
260:             if (nxt > 0) llprev[nxt] = prv;
261:             if (prv > 0) {
262:               llnext[prv] = nxt;
263:             } else {
264:               degb[degrees[idx]] = nxt;
265:             }
266:             degrees[idx]--;
267:             llnext[idx] = degb[degrees[idx]];
268:             llprev[idx] = -1;
269:             if (degb[degrees[idx]] >= 0) llprev[degb[degrees[idx]]] = idx;
270:             degb[degrees[idx]] = idx;
271:             if (dist < distance) {
272:               ncols = Gi[idx + 1] - Gi[idx];
273:               cols  = &(Gj[Gi[idx]]);
274:               for (j = 0; j < ncols; j++) {
275:                 if (seen[cols[j]] != i) {
276:                   bidx++;
277:                   seen[cols[j]] = i;
278:                   idxbuf[bidx]  = cols[j];
279:                   distbuf[bidx] = dist + 1;
280:                 }
281:               }
282:             }
283:           }
284:         }
285:         break;
286:       }
287:     }
288:   }
289:   for (i = 0; i < lm; i++) {
290:     if (gidx[i] >= s && gidx[i] < e) weights[gidx[i] - s] = lweights[i];
291:   }
292:   PetscRandomDestroy(&rand);
293:   PetscFree(degb);
294:   PetscFree2(llnext, llprev);
295:   PetscFree(degrees);
296:   PetscFree(lweights);
297:   ISRestoreIndices(ris, &gidx);
298:   ISDestroy(&ris);
299:   PetscFree3(seen, idxbuf, distbuf);
300:   MatDestroyMatrices(1, &lGs);
301:   return 0;
302: }

304: PetscErrorCode MatColoringCreateWeights(MatColoring mc, PetscReal **weights, PetscInt **lperm)
305: {
306:   PetscInt   i, s, e, n;
307:   PetscReal *wts;

309:   /* create weights of the specified type */
310:   MatGetOwnershipRange(mc->mat, &s, &e);
311:   n = e - s;
312:   PetscMalloc1(n, &wts);
313:   switch (mc->weight_type) {
314:   case MAT_COLORING_WEIGHT_RANDOM:
315:     MatColoringCreateRandomWeights(mc, wts);
316:     break;
317:   case MAT_COLORING_WEIGHT_LEXICAL:
318:     MatColoringCreateLexicalWeights(mc, wts);
319:     break;
320:   case MAT_COLORING_WEIGHT_LF:
321:     MatColoringCreateLargestFirstWeights(mc, wts);
322:     break;
323:   case MAT_COLORING_WEIGHT_SL:
324:     MatColoringCreateSmallestLastWeights(mc, wts);
325:     break;
326:   }
327:   if (lperm) {
328:     PetscMalloc1(n, lperm);
329:     for (i = 0; i < n; i++) (*lperm)[i] = i;
330:     PetscSortRealWithPermutation(n, wts, *lperm);
331:     for (i = 0; i < n / 2; i++) {
332:       PetscInt swp;
333:       swp                 = (*lperm)[i];
334:       (*lperm)[i]         = (*lperm)[n - 1 - i];
335:       (*lperm)[n - 1 - i] = swp;
336:     }
337:   }
338:   if (weights) *weights = wts;
339:   return 0;
340: }

342: PetscErrorCode MatColoringSetWeights(MatColoring mc, PetscReal *weights, PetscInt *lperm)
343: {
344:   PetscInt i, s, e, n;

346:   MatGetOwnershipRange(mc->mat, &s, &e);
347:   n = e - s;
348:   if (weights) {
349:     PetscMalloc2(n, &mc->user_weights, n, &mc->user_lperm);
350:     for (i = 0; i < n; i++) mc->user_weights[i] = weights[i];
351:     if (!lperm) {
352:       for (i = 0; i < n; i++) mc->user_lperm[i] = i;
353:       PetscSortRealWithPermutation(n, mc->user_weights, mc->user_lperm);
354:       for (i = 0; i < n / 2; i++) {
355:         PetscInt swp;
356:         swp                       = mc->user_lperm[i];
357:         mc->user_lperm[i]         = mc->user_lperm[n - 1 - i];
358:         mc->user_lperm[n - 1 - i] = swp;
359:       }
360:     } else {
361:       for (i = 0; i < n; i++) mc->user_lperm[i] = lperm[i];
362:     }
363:   } else {
364:     mc->user_weights = NULL;
365:     mc->user_lperm   = NULL;
366:   }
367:   return 0;
368: }