Actual source code: weights.c

  1: #include <petsc/private/matimpl.h>
  2: #include <../src/mat/impls/aij/seq/aij.h>

  4: static PetscErrorCode MatColoringCreateLexicalWeights(MatColoring mc, PetscReal *weights)
  5: {
  6:   PetscInt i, s, e;
  7:   Mat      G = mc->mat;

  9:   PetscFunctionBegin;
 10:   PetscCall(MatGetOwnershipRange(G, &s, &e));
 11:   for (i = s; i < e; i++) weights[i - s] = i;
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

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

 22:   PetscFunctionBegin;
 23:   /* each weight should be the degree plus a random perturbation */
 24:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)mc), &rand));
 25:   PetscCall(PetscRandomSetFromOptions(rand));
 26:   PetscCall(MatGetOwnershipRange(G, &s, &e));
 27:   for (i = s; i < e; i++) {
 28:     PetscCall(PetscRandomGetValueReal(rand, &r));
 29:     weights[i - s] = PetscAbsReal(r);
 30:   }
 31:   PetscCall(PetscRandomDestroy(&rand));
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBegin;
 51:   PetscCall(MatGetOwnershipRange(G, &s, &e));
 52:   n = e - s;
 53:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)G), n, s, 1, &ris));
 54:   PetscCall(MatIncreaseOverlap(G, 1, &ris, distance));
 55:   PetscCall(ISSort(ris));
 56:   PetscCall(MatCreateSubMatrices(G, 1, &ris, &ris, MAT_INITIAL_MATRIX, &lGs));
 57:   lG = lGs[0];
 58:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)lG, MATSEQAIJ, &isSEQAIJ));
 59:   PetscCheck(isSEQAIJ, PetscObjectComm((PetscObject)G), PETSC_ERR_SUP, "Requires an MPI/SEQAIJ Matrix");
 60:   PetscCall(MatGetSize(lG, &ln, &lm));
 61:   aij = (Mat_SeqAIJ *)lG->data;
 62:   Gi  = aij->i;
 63:   Gj  = aij->j;
 64:   PetscCall(PetscMalloc3(lm, &seen, lm, &idxbuf, lm, &distbuf));
 65:   for (i = 0; i < ln; i++) seen[i] = -1;
 66:   PetscCall(ISGetIndices(ris, &gidx));
 67:   for (i = 0; i < ln; i++) {
 68:     if (gidx[i] >= e || gidx[i] < s) continue;
 69:     bidx   = -1;
 70:     ncols  = Gi[i + 1] - Gi[i];
 71:     cols   = &(Gj[Gi[i]]);
 72:     degree = 0;
 73:     /* place the distance-one neighbors on the queue */
 74:     for (j = 0; j < ncols; j++) {
 75:       bidx++;
 76:       seen[cols[j]] = i;
 77:       distbuf[bidx] = 1;
 78:       idxbuf[bidx]  = cols[j];
 79:     }
 80:     while (bidx >= 0) {
 81:       /* pop */
 82:       idx  = idxbuf[bidx];
 83:       dist = distbuf[bidx];
 84:       bidx--;
 85:       degree++;
 86:       if (dist < distance) {
 87:         ncols = Gi[idx + 1] - Gi[idx];
 88:         cols  = &(Gj[Gi[idx]]);
 89:         for (j = 0; j < ncols; j++) {
 90:           if (seen[cols[j]] != i) {
 91:             bidx++;
 92:             seen[cols[j]] = i;
 93:             idxbuf[bidx]  = cols[j];
 94:             distbuf[bidx] = dist + 1;
 95:           }
 96:         }
 97:       }
 98:     }
 99:     degrees[gidx[i] - s] = degree;
100:   }
101:   PetscCall(ISRestoreIndices(ris, &gidx));
102:   PetscCall(ISDestroy(&ris));
103:   PetscCall(PetscFree3(seen, idxbuf, distbuf));
104:   PetscCall(MatDestroyMatrices(1, &lGs));
105:   PetscFunctionReturn(PETSC_SUCCESS);
106: }

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

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

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

309: PetscErrorCode MatColoringCreateWeights(MatColoring mc, PetscReal **weights, PetscInt **lperm)
310: {
311:   PetscInt   i, s, e, n;
312:   PetscReal *wts;

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

348: PetscErrorCode MatColoringSetWeights(MatColoring mc, PetscReal *weights, PetscInt *lperm)
349: {
350:   PetscInt i, s, e, n;

352:   PetscFunctionBegin;
353:   PetscCall(MatGetOwnershipRange(mc->mat, &s, &e));
354:   n = e - s;
355:   if (weights) {
356:     PetscCall(PetscMalloc2(n, &mc->user_weights, n, &mc->user_lperm));
357:     for (i = 0; i < n; i++) mc->user_weights[i] = weights[i];
358:     if (!lperm) {
359:       for (i = 0; i < n; i++) mc->user_lperm[i] = i;
360:       PetscCall(PetscSortRealWithPermutation(n, mc->user_weights, mc->user_lperm));
361:       for (i = 0; i < n / 2; i++) {
362:         PetscInt swp;
363:         swp                       = mc->user_lperm[i];
364:         mc->user_lperm[i]         = mc->user_lperm[n - 1 - i];
365:         mc->user_lperm[n - 1 - i] = swp;
366:       }
367:     } else {
368:       for (i = 0; i < n; i++) mc->user_lperm[i] = lperm[i];
369:     }
370:   } else {
371:     mc->user_weights = NULL;
372:     mc->user_lperm   = NULL;
373:   }
374:   PetscFunctionReturn(PETSC_SUCCESS);
375: }