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: {
      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: }