Actual source code: hem.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  2: #include <petsc/private/matimpl.h>    /*I "petscmat.h" I*/
  3: #include <../src/mat/impls/aij/seq/aij.h>
  4: #include <../src/mat/impls/aij/mpi/mpiaij.h>

  6: /* linked list methods
  7:  *
  8:  *  PetscCDCreate
  9:  */
 12: PetscErrorCode PetscCDCreate(PetscInt a_size, PetscCoarsenData **a_out)
 13: {
 14:   PetscErrorCode   ierr;
 15:   PetscCoarsenData *ail;
 16:   PetscInt         ii;

 19:   /* alocate pool, partially */
 20:   PetscNew(&ail);
 21:   *a_out               = ail;
 22:   ail->pool_list.next  = NULL;
 23:   ail->pool_list.array = NULL;
 24:   ail->chk_sz          = 0;
 25:   /* allocate array */
 26:   ail->size = a_size;
 27:   PetscMalloc1(a_size, &ail->array);
 28:   for (ii=0; ii<a_size; ii++) ail->array[ii] = NULL;
 29:   ail->extra_nodes = NULL;
 30:   ail->mat         = NULL;
 31:   /* ail->removedIS = NULL; */
 32:   return(0);
 33: }

 35: /* NPDestroy
 36:  */
 39: PetscErrorCode PetscCDDestroy(PetscCoarsenData *ail)
 40: {
 42:   PetscCDArrNd   *n = &ail->pool_list;

 45:   n = n->next;
 46:   while (n) {
 47:     PetscCDArrNd *lstn = n;
 48:     n    = n->next;
 49:     PetscFree(lstn);
 50:   }
 51:   if (ail->pool_list.array) {
 52:     PetscFree(ail->pool_list.array);
 53:   }
 54:   /* if (ail->removedIS) { */
 55:   /*   ISDestroy(&ail->removedIS); */
 56:   /* } */
 57:   /* delete this (+array) */
 58:   PetscFree(ail->array);
 59:   /* delete this (+agg+pool array) */
 60:   PetscFree(ail);
 61:   return(0);
 62: }
 63: /* PetscCDSetChuckSize
 64:  */
 67: PetscErrorCode PetscCDSetChuckSize(PetscCoarsenData *ail, PetscInt a_sz)
 68: {
 70:   ail->chk_sz = a_sz;
 71:   return(0);
 72: }
 73: /*  PetscCDGetNewNode
 74:  */
 77: PetscErrorCode PetscCDGetNewNode(PetscCoarsenData *ail, PetscCDIntNd **a_out, PetscInt a_id)
 78: {

 82:   if (ail->extra_nodes) {
 83:     PetscCDIntNd *node = ail->extra_nodes;
 84:     ail->extra_nodes = node->next;
 85:     node->gid        = a_id;
 86:     node->next       = NULL;
 87:     *a_out           = node;
 88:   } else {
 89:     if (!ail->pool_list.array) {
 90:       if (!ail->chk_sz) ail->chk_sz = 10; /* use a chuck size of ail->size? */
 91:       PetscMalloc1(ail->chk_sz, &ail->pool_list.array);
 92:       ail->new_node       = ail->pool_list.array;
 93:       ail->new_left       = ail->chk_sz;
 94:       ail->new_node->next = NULL;
 95:     } else if (!ail->new_left) {
 96:       PetscCDArrNd *node;
 97:       PetscMalloc(ail->chk_sz*sizeof(PetscCDIntNd) + sizeof(PetscCDArrNd), &node);
 98:       node->array         = (PetscCDIntNd*)(node + 1);
 99:       node->next          = ail->pool_list.next;
100:       ail->pool_list.next = node;
101:       ail->new_left       = ail->chk_sz;
102:       ail->new_node       = node->array;
103:     }
104:     ail->new_node->gid  = a_id;
105:     ail->new_node->next = NULL;
106:     *a_out              = ail->new_node++; ail->new_left--;
107:   }
108:   return(0);
109: }

111: /* PetscLLNSetID
112:  */
115: PetscErrorCode PetscLLNSetID(PetscCDIntNd *a_this, PetscInt a_id)
116: {
118:   a_this->gid = a_id;
119:   return(0);
120: }
121: /* PetscLLNGetID
122:  */
125: PetscErrorCode PetscLLNGetID(const PetscCDIntNd *a_this, PetscInt *a_gid)
126: {
128:   *a_gid = a_this->gid;
129:   return(0);
130: }
131: /* PetscCDGetHeadPos
132:  */
135: PetscErrorCode PetscCDGetHeadPos(const PetscCoarsenData *ail, PetscInt a_idx, PetscCDPos *pos)
136: {
138:   if (a_idx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"a_idx >= ail->size: a_idx=%d.",a_idx);
139:   *pos = ail->array[a_idx];
140:   return(0);
141: }
142: /* PetscCDGetNextPos
143:  */
146: PetscErrorCode PetscCDGetNextPos(const PetscCoarsenData *ail, PetscInt l_idx, PetscCDPos *pos)
147: {
149:   if (!(*pos)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"NULL input position.");
150:   *pos = (*pos)->next;
151:   return(0);
152: }

154: /* PetscCDAppendID
155:  */
158: PetscErrorCode PetscCDAppendID(PetscCoarsenData *ail, PetscInt a_idx, PetscInt a_id)
159: {
161:   PetscCDIntNd   *n,*n2;

164:   PetscCDGetNewNode(ail, &n, a_id);
165:   if (a_idx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_idx);
166:   if (!(n2=ail->array[a_idx])) ail->array[a_idx] = n;
167:   else {
168:     do {
169:       if (!n2->next) {
170:         n2->next = n;
171:         if (n->next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"n should not have a next");
172:         break;
173:       }
174:       n2 = n2->next;
175:     } while (n2);
176:     if (!n2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"n2 should be non-null");
177:   }
178:   return(0);
179: }
180: /* PetscCDAppendNode
181:  */
184: PetscErrorCode PetscCDAppendNode(PetscCoarsenData *ail, PetscInt a_idx,  PetscCDIntNd *a_n)
185: {
186:   PetscCDIntNd *n2;

189:   if (a_idx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_idx);
190:   if (!(n2=ail->array[a_idx])) ail->array[a_idx] = a_n;
191:   else {
192:     do {
193:       if (!n2->next) {
194:         n2->next  = a_n;
195:         a_n->next = NULL;
196:         break;
197:       }
198:       n2 = n2->next;
199:     } while (n2);
200:     if (!n2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"n2 should be non-null");
201:   }
202:   return(0);
203: }

205: /* PetscCDRemoveNextNode: a_last->next, this exposes single linked list structure to API
206:  */
209: PetscErrorCode PetscCDRemoveNextNode(PetscCoarsenData *ail, PetscInt a_idx,  PetscCDIntNd *a_last)
210: {
211:   PetscCDIntNd *del;

214:   if (a_idx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_idx);
215:   if (!a_last->next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"a_last should have a next");
216:   del          = a_last->next;
217:   a_last->next = del->next;
218:   /* del->next = NULL; -- this still used in a iterator so keep it intact -- need to fix this with a double linked list */
219:   /* could reuse n2 but PetscCDAppendNode sometimes uses it */
220:   return(0);
221: }

223: /* PetscCDPrint
224:  */
227: PetscErrorCode PetscCDPrint(const PetscCoarsenData *ail, MPI_Comm comm)
228: {
230:   PetscCDIntNd   *n;
231:   PetscInt       ii,kk;
232:   PetscMPIInt    rank;

235:   MPI_Comm_rank(comm, &rank);
236:   for (ii=0; ii<ail->size; ii++) {
237:     kk = 0;
238:     n  = ail->array[ii];
239:     if (n) PetscPrintf(comm,"[%d]%s list %d:\n",rank,__FUNCT__,ii);
240:     while (n) {
241:       PetscPrintf(comm,"\t[%d] %d) id %d\n",rank,++kk,n->gid);
242:       n = n->next;
243:     }
244:   }
245:   return(0);
246: }
247: /* PetscCDAppendRemove
248:  */
251: PetscErrorCode PetscCDAppendRemove(PetscCoarsenData *ail, PetscInt a_destidx, PetscInt a_srcidx)
252: {
253:   PetscCDIntNd *n;

256:   if (a_srcidx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_srcidx);
257:   if (a_destidx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_destidx);
258:   if (a_destidx==a_srcidx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"a_destidx==a_srcidx %d.",a_destidx);
259:   n = ail->array[a_destidx];
260:   if (!n) ail->array[a_destidx] = ail->array[a_srcidx];
261:   else {
262:     do {
263:       if (!n->next) {
264:         n->next = ail->array[a_srcidx];
265:         break;
266:       }
267:       n = n->next;
268:     } while (1);
269:   }
270:   ail->array[a_srcidx] = NULL;
271:   return(0);
272: }

274: /* PetscCDRemoveAll
275:  */
278: PetscErrorCode PetscCDRemoveAll(PetscCoarsenData *ail, PetscInt a_idx)
279: {
280:   PetscCDIntNd *rem,*n1;

283:   if (a_idx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_idx);
284:   rem               = ail->array[a_idx];
285:   ail->array[a_idx] = NULL;
286:   if (!(n1=ail->extra_nodes)) ail->extra_nodes = rem;
287:   else {
288:     while (n1->next) n1 = n1->next;
289:     n1->next = rem;
290:   }
291:   return(0);
292: }

294: /* PetscCDSizeAt
295:  */
298: PetscErrorCode PetscCDSizeAt(const PetscCoarsenData *ail, PetscInt a_idx, PetscInt *a_sz)
299: {
300:   PetscCDIntNd *n1;
301:   PetscInt     sz = 0;

304:   if (a_idx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_idx);
305:   n1 = ail->array[a_idx];
306:   while (n1) {
307:     n1 = n1->next;
308:     sz++;
309:   }
310:   *a_sz = sz;
311:   return(0);
312: }

314: /* PetscCDEmptyAt
315:  */
318: PetscErrorCode PetscCDEmptyAt(const PetscCoarsenData *ail, PetscInt a_idx, PetscBool *a_e)
319: {
321:   if (a_idx>=ail->size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Index %d out of range.",a_idx);
322:   *a_e = (PetscBool)(ail->array[a_idx]==NULL);
323:   return(0);
324: }

326: /* PetscCDGetMIS
327:  */
330: PetscErrorCode PetscCDGetMIS(PetscCoarsenData *ail, IS *a_mis)
331: {
333:   PetscCDIntNd   *n;
334:   PetscInt       ii,kk;
335:   PetscInt       *permute;

338:   for (ii=kk=0; ii<ail->size; ii++) {
339:     n = ail->array[ii];
340:     if (n) kk++;
341:   }
342:   PetscMalloc1(kk, &permute);
343:   for (ii=kk=0; ii<ail->size; ii++) {
344:     n = ail->array[ii];
345:     if (n) permute[kk++] = ii;
346:   }
347:   ISCreateGeneral(PETSC_COMM_SELF, kk, permute, PETSC_OWN_POINTER, a_mis);
348:   return(0);
349: }
350: /* PetscCDGetMat
351:  */
354: PetscErrorCode PetscCDGetMat(const PetscCoarsenData *ail, Mat *a_mat)
355: {
357:   *a_mat = ail->mat;
358:   return(0);
359: }

361: /* PetscCDSetMat
362:  */
365: PetscErrorCode PetscCDSetMat(PetscCoarsenData *ail, Mat a_mat)
366: {
368:   ail->mat = a_mat;
369:   return(0);
370: }


373: /* PetscCDGetASMBlocks
374:  */
377: PetscErrorCode PetscCDGetASMBlocks(const PetscCoarsenData *ail, const PetscInt a_bs, PetscInt *a_sz, IS **a_local_is)
378: {
380:   PetscCDIntNd   *n;
381:   PetscInt       lsz,ii,kk,*idxs,jj;
382:   IS             *is_loc;

385:   for (ii=kk=0; ii<ail->size; ii++) {
386:     if (ail->array[ii]) kk++;
387:   }
388:   *a_sz = kk; /* out */

390:   PetscMalloc1(kk, &is_loc);
391:   *a_local_is = is_loc; /* out */

393:   for (ii=kk=0; ii<ail->size; ii++) {
394:     for (lsz=0, n=ail->array[ii]; n; lsz++, n=n->next) /* void */;
395:     if (lsz) {
396:       PetscMalloc1(a_bs*lsz, &idxs);
397:       for (lsz = 0, n=ail->array[ii]; n; n = n->next) {
398:         PetscInt gid;
399:         PetscLLNGetID(n, &gid);
400:         for (jj=0; jj<a_bs; lsz++,jj++) idxs[lsz] = a_bs*gid + jj;
401:       }
402:       ISCreateGeneral(PETSC_COMM_SELF, lsz, idxs, PETSC_OWN_POINTER, &is_loc[kk++]);
403:     }
404:   }
405:   if (*a_sz != kk) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"*a_sz %D != kk %D",*a_sz,kk);
406:   return(0);
407: }


410: /* PetscCDSetRemovedIS
411:  */
412: /* PetscErrorCode PetscCDSetRemovedIS(PetscCoarsenData *ail, MPI_Comm comm, const PetscInt a_sz, PetscInt a_ids[]) */
413: /* { */

416: /*   if (ail->removedIS) { */
417: /*     ISDestroy(&ail->removedIS); */
418: /*   } */
419: /*   ISCreateGeneral(comm, a_sz, a_ids, PETSC_COPY_VALUES, &ail->removedIS); */
420: /*   return(0); */
421: /* } */

423: /* PetscCDGetRemovedIS
424:  */
425: /* PetscErrorCode PetscCDGetRemovedIS(PetscCoarsenData *ail, IS *a_is) */
426: /* { */
427: /*   *a_is = ail->removedIS; */
428: /*   ail->removedIS = NULL; /\* hack to relinquish ownership *\/ */
429: /*   return(0); */
430: /* } */

432: /* ********************************************************************** */
433: /* edge for priority queue */
434: typedef struct edge_tag {
435:   PetscReal weight;
436:   PetscInt  lid0,gid1,cpid1;
437: } Edge;

439: static int gamg_hem_compare(const void *a, const void *b)
440: {
441:   PetscReal va = ((Edge*)a)->weight, vb = ((Edge*)b)->weight;
442:   return (va < vb) ? 1 : (va == vb) ? 0 : -1; /* 0 for equal */
443: }

445: /* -------------------------------------------------------------------------- */
446: /*
447:    heavyEdgeMatchAgg - parallel heavy edge matching (HEM). MatAIJ specific!!!

449:    Input Parameter:
450:    . perm - permutation
451:    . a_Gmat - glabal matrix of graph (data not defined)

453:    Output Parameter:
454:    . a_locals_llist - array of list of local nodes rooted at local node
455: */
458: static PetscErrorCode heavyEdgeMatchAgg(IS perm,Mat a_Gmat,PetscCoarsenData **a_locals_llist)
459: {
460:   PetscErrorCode   ierr;
461:   PetscBool        isMPI;
462:   MPI_Comm         comm;
463:   PetscInt         sub_it,kk,n,ix,*idx,*ii,iter,Iend,my0;
464:   PetscMPIInt      rank,size;
465:   const PetscInt   nloc = a_Gmat->rmap->n,n_iter=6; /* need to figure out how to stop this */
466:   PetscInt         *lid_cprowID,*lid_gid;
467:   PetscBool        *lid_matched;
468:   Mat_SeqAIJ       *matA, *matB=0;
469:   Mat_MPIAIJ       *mpimat     =0;
470:   PetscScalar      one         =1.;
471:   PetscCoarsenData *agg_llists = NULL,*deleted_list = NULL;
472:   Mat              cMat,tMat,P;
473:   MatScalar        *ap;
474:   PetscMPIInt      tag1,tag2;

477:   PetscObjectGetComm((PetscObject)a_Gmat,&comm);
478:   MPI_Comm_rank(comm, &rank);
479:   MPI_Comm_size(comm, &size);
480:   MatGetOwnershipRange(a_Gmat, &my0, &Iend);
481:   PetscCommGetNewTag(comm, &tag1);
482:   PetscCommGetNewTag(comm, &tag2);

484:   PetscMalloc1(nloc, &lid_gid); /* explicit array needed */
485:   PetscMalloc1(nloc, &lid_cprowID);
486:   PetscMalloc1(nloc, &lid_matched);

488:   PetscCDCreate(nloc, &agg_llists);
489:   /* PetscCDSetChuckSize(agg_llists, nloc+1); */
490:   *a_locals_llist = agg_llists;
491:   PetscCDCreate(size, &deleted_list);
492:   PetscCDSetChuckSize(deleted_list, 100);
493:   /* setup 'lid_gid' for scatters and add self to all lists */
494:   for (kk=0; kk<nloc; kk++) {
495:     lid_gid[kk] = kk + my0;
496:     PetscCDAppendID(agg_llists, kk, my0+kk);
497:   }

499:   /* make a copy of the graph, this gets destroyed in iterates */
500:   MatDuplicate(a_Gmat,MAT_COPY_VALUES,&cMat);
501:   PetscObjectTypeCompare((PetscObject)a_Gmat, MATMPIAIJ, &isMPI);
502:   iter = 0;
503:   while (iter++ < n_iter) {
504:     PetscScalar    *cpcol_gid,*cpcol_max_ew,*cpcol_max_pe,*lid_max_ew;
505:     PetscBool      *cpcol_matched;
506:     PetscMPIInt    *cpcol_pe,proc;
507:     Vec            locMaxEdge,locMaxPE,ghostMaxEdge,ghostMaxPE;
508:     PetscInt       nEdges,n_nz_row,jj;
509:     Edge           *Edges;
510:     PetscInt       gid;
511:     const PetscInt *perm_ix, n_sub_its = 120;

513:     /* get submatrices of cMat */
514:     if (isMPI) {
515:       mpimat = (Mat_MPIAIJ*)cMat->data;
516:       matA   = (Mat_SeqAIJ*)mpimat->A->data;
517:       matB   = (Mat_SeqAIJ*)mpimat->B->data;
518:       if (!matB->compressedrow.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"matB must have compressed row usage");
519:     } else {
520:       matA = (Mat_SeqAIJ*)cMat->data;
521:     }

523:     /* set max edge on nodes */
524:     MatCreateVecs(cMat, &locMaxEdge, 0);
525:     MatCreateVecs(cMat, &locMaxPE, 0);

527:     /* get 'cpcol_pe' & 'cpcol_gid' & init. 'cpcol_matched' using 'mpimat->lvec' */
528:     if (mpimat) {
529:       Vec         vec;
530:       PetscScalar vval;

532:       MatCreateVecs(cMat, &vec, 0);
533:       /* cpcol_pe */
534:       vval = (PetscScalar)(rank);
535:       for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
536:         VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
537:       }
538:       VecAssemblyBegin(vec);
539:       VecAssemblyEnd(vec);
540:       VecScatterBegin(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
541:       VecScatterEnd(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
542:       VecGetArray(mpimat->lvec, &cpcol_gid); /* get proc ID in 'cpcol_gid' */
543:       VecGetLocalSize(mpimat->lvec, &n);
544:       PetscMalloc1(n, &cpcol_pe);
545:       for (kk=0; kk<n; kk++) cpcol_pe[kk] = (PetscMPIInt)PetscRealPart(cpcol_gid[kk]);
546:       VecRestoreArray(mpimat->lvec, &cpcol_gid);

548:       /* cpcol_gid */
549:       for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
550:         vval = (PetscScalar)(gid);
551:         VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
552:       }
553:       VecAssemblyBegin(vec);
554:       VecAssemblyEnd(vec);
555:       VecScatterBegin(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
556:       VecScatterEnd(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
557:       VecDestroy(&vec);
558:       VecGetArray(mpimat->lvec, &cpcol_gid); /* get proc ID in 'cpcol_gid' */

560:       /* cpcol_matched */
561:       VecGetLocalSize(mpimat->lvec, &n);
562:       PetscMalloc1(n, &cpcol_matched);
563:       for (kk=0; kk<n; kk++) cpcol_matched[kk] = PETSC_FALSE;
564:     }

566:     /* need an inverse map - locals */
567:     for (kk=0; kk<nloc; kk++) lid_cprowID[kk] = -1;
568:     /* set index into compressed row 'lid_cprowID' */
569:     if (matB) {
570:       ii = matB->compressedrow.i;
571:       for (ix=0; ix<matB->compressedrow.nrows; ix++) {
572:         lid_cprowID[matB->compressedrow.rindex[ix]] = ix;
573:       }
574:     }

576:     /* get removed IS, use '' */
577:     /* if (iter==1) { */
578:     /*   PetscInt *lid_rem,idx; */
579:     /*   PetscMalloc1(nloc, &lid_rem); */
580:     /*   for (kk=idx=0;kk<nloc;kk++) { */
581:     /*     PetscInt nn,lid=kk; */
582:     /*     ii = matA->i; nn = ii[lid+1] - ii[lid]; */
583:     /*     if ((ix=lid_cprowID[lid]) != -1) { /\* if I have any ghost neighbors *\/ */
584:     /*       ii = matB->compressedrow.i; */
585:     /*       nn += ii[ix+1] - ii[ix]; */
586:     /*     } */
587:     /*     if (nn < 2) { */
588:     /*       lid_rem[idx++] = kk + my0; */
589:     /*     } */
590:     /*   } */
591:     /*   PetscCDSetRemovedIS(agg_llists, comm, idx, lid_rem); */
592:     /*   PetscFree(lid_rem); */
593:     /* } */

595:     /* compute 'locMaxEdge' & 'locMaxPE', and create list of edges, count edges' */
596:     for (nEdges=0,kk=0,gid=my0; kk<nloc; kk++,gid++) {
597:       PetscReal   max_e = 0., tt;
598:       PetscScalar vval;
599:       PetscInt    lid   = kk;
600:       PetscMPIInt max_pe=rank,pe;
601:       ii = matA->i; n = ii[lid+1] - ii[lid]; idx = matA->j + ii[lid];
602:       ap = matA->a + ii[lid];
603:       for (jj=0; jj<n; jj++) {
604:         PetscInt lidj = idx[jj];
605:         if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
606:         if (lidj > lid) nEdges++;
607:       }
608:       if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
609:         ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
610:         ap  = matB->a + ii[ix];
611:         idx = matB->j + ii[ix];
612:         for (jj=0; jj<n; jj++) {
613:           if ((tt=PetscRealPart(ap[jj])) > max_e) max_e = tt;
614:           nEdges++;
615:           if ((pe=cpcol_pe[idx[jj]]) > max_pe) max_pe = pe;
616:         }
617:       }
618:       vval = max_e;
619:       VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES);

621:       vval = (PetscScalar)max_pe;
622:       VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES);
623:     }
624:     VecAssemblyBegin(locMaxEdge);
625:     VecAssemblyEnd(locMaxEdge);
626:     VecAssemblyBegin(locMaxPE);
627:     VecAssemblyEnd(locMaxPE);

629:     /* get 'cpcol_max_ew' & 'cpcol_max_pe' */
630:     if (mpimat) {
631:       VecDuplicate(mpimat->lvec, &ghostMaxEdge);
632:       VecScatterBegin(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
633:       VecScatterEnd(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
634:       VecGetArray(ghostMaxEdge, &cpcol_max_ew);

636:       VecDuplicate(mpimat->lvec, &ghostMaxPE);
637:       VecScatterBegin(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
638:         VecScatterEnd(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
639:       VecGetArray(ghostMaxPE, &cpcol_max_pe);
640:     }

642:     /* setup sorted list of edges */
643:     PetscMalloc1(nEdges, &Edges);
644:     ISGetIndices(perm, &perm_ix);
645:     for (nEdges=n_nz_row=kk=0; kk<nloc; kk++) {
646:       PetscInt nn, lid = perm_ix[kk];
647:       ii = matA->i; nn = n = ii[lid+1] - ii[lid]; idx = matA->j + ii[lid];
648:       ap = matA->a + ii[lid];
649:       for (jj=0; jj<n; jj++) {
650:         PetscInt lidj = idx[jj];
651:         if (lidj > lid) {
652:           Edges[nEdges].lid0   = lid;
653:           Edges[nEdges].gid1   = lidj + my0;
654:           Edges[nEdges].cpid1  = -1;
655:           Edges[nEdges].weight = PetscRealPart(ap[jj]);
656:           nEdges++;
657:         }
658:       }
659:       if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
660:         ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
661:         ap  = matB->a + ii[ix];
662:         idx = matB->j + ii[ix];
663:         nn += n;
664:         for (jj=0; jj<n; jj++) {
665:           Edges[nEdges].lid0   = lid;
666:           Edges[nEdges].gid1   = (PetscInt)PetscRealPart(cpcol_gid[idx[jj]]);
667:           Edges[nEdges].cpid1  = idx[jj];
668:           Edges[nEdges].weight = PetscRealPart(ap[jj]);
669:           nEdges++;
670:         }
671:       }
672:       if (nn > 1) n_nz_row++;
673:       else if (iter == 1) {
674:         /* should select this because it is technically in the MIS but lets not */
675:         PetscCDRemoveAll(agg_llists, lid);
676:       }
677:     }
678:     ISRestoreIndices(perm,&perm_ix);

680:     qsort(Edges, nEdges, sizeof(Edge), gamg_hem_compare);

682:     /* projection matrix */
683:     MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 1, 0, 1, 0, &P);

685:     /* clear matched flags */
686:     for (kk=0; kk<nloc; kk++) lid_matched[kk] = PETSC_FALSE;
687:     /* process - communicate - process */
688:     for (sub_it=0; sub_it<n_sub_its; sub_it++) {
689:       PetscInt nactive_edges;

691:       VecGetArray(locMaxEdge, &lid_max_ew);
692:       for (kk=nactive_edges=0; kk<nEdges; kk++) {
693:         /* HEM */
694:         const Edge     *e   = &Edges[kk];
695:         const PetscInt lid0 =e->lid0,gid1=e->gid1,cpid1=e->cpid1,gid0=lid0+my0,lid1=gid1-my0;
696:         PetscBool      isOK = PETSC_TRUE;

698:         /* skip if either (local) vertex is done already */
699:         if (lid_matched[lid0] || (gid1>=my0 && gid1<Iend && lid_matched[gid1-my0])) continue;

701:         /* skip if ghost vertex is done */
702:         if (cpid1 != -1 && cpcol_matched[cpid1]) continue;

704:         nactive_edges++;
705:         /* skip if I have a bigger edge someplace (lid_max_ew gets updated) */
706:         if (PetscRealPart(lid_max_ew[lid0]) > e->weight + 1.e-12) continue;

708:         if (cpid1 == -1) {
709:           if (PetscRealPart(lid_max_ew[lid1]) > e->weight + 1.e-12) continue;
710:         } else {
711:           /* see if edge might get matched on other proc */
712:           PetscReal g_max_e = PetscRealPart(cpcol_max_ew[cpid1]);
713:           if (g_max_e > e->weight + 1.e-12) continue;
714:           else if (e->weight > g_max_e - 1.e-12 && (PetscMPIInt)PetscRealPart(cpcol_max_pe[cpid1]) > rank) {
715:             /* check for max_e == to this edge and larger processor that will deal with this */
716:             continue;
717:           }
718:         }

720:         /* check ghost for v0 */
721:         if (isOK) {
722:           PetscReal max_e,ew;
723:           if ((ix=lid_cprowID[lid0]) != -1) { /* if I have any ghost neighbors */
724:             ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
725:             ap  = matB->a + ii[ix];
726:             idx = matB->j + ii[ix];
727:             for (jj=0; jj<n && isOK; jj++) {
728:               PetscInt lidj = idx[jj];
729:               if (cpcol_matched[lidj]) continue;
730:               ew = PetscRealPart(ap[jj]); max_e = PetscRealPart(cpcol_max_ew[lidj]);
731:               /* check for max_e == to this edge and larger processor that will deal with this */
732:               if (ew > max_e - 1.e-12 && ew > PetscRealPart(lid_max_ew[lid0]) - 1.e-12 && (PetscMPIInt)PetscRealPart(cpcol_max_pe[lidj]) > rank) {
733:                 isOK = PETSC_FALSE;
734:               }
735:             }
736:           }

738:           /* for v1 */
739:           if (cpid1 == -1 && isOK) {
740:             if ((ix=lid_cprowID[lid1]) != -1) { /* if I have any ghost neighbors */
741:               ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
742:               ap  = matB->a + ii[ix];
743:               idx = matB->j + ii[ix];
744:               for (jj=0; jj<n && isOK; jj++) {
745:                 PetscInt lidj = idx[jj];
746:                 if (cpcol_matched[lidj]) continue;
747:                 ew = PetscRealPart(ap[jj]); max_e = PetscRealPart(cpcol_max_ew[lidj]);
748:                 /* check for max_e == to this edge and larger processor that will deal with this */
749:                 if (ew > max_e - 1.e-12 && ew > PetscRealPart(lid_max_ew[lid1]) - 1.e-12 && (PetscMPIInt)PetscRealPart(cpcol_max_pe[lidj]) > rank) {
750:                   isOK = PETSC_FALSE;
751:                 }
752:               }
753:             }
754:           }
755:         }

757:         /* do it */
758:         if (isOK) {
759:           if (cpid1 == -1) {
760:             lid_matched[lid1] = PETSC_TRUE;  /* keep track of what we've done this round */
761:             PetscCDAppendRemove(agg_llists, lid0, lid1);
762:           } else if (sub_it != n_sub_its-1) {
763:             /* add gid1 to list of ghost deleted by me -- I need their children */
764:             proc = cpcol_pe[cpid1];

766:             cpcol_matched[cpid1] = PETSC_TRUE; /* messing with VecGetArray array -- needed??? */

768:             PetscCDAppendID(deleted_list, proc, cpid1); /* cache to send messages */
769:             PetscCDAppendID(deleted_list, proc, lid0);
770:           } else continue;

772:           lid_matched[lid0] = PETSC_TRUE; /* keep track of what we've done this round */
773:           /* set projection */
774:           MatSetValues(P,1,&gid0,1,&gid0,&one,INSERT_VALUES);
775:           MatSetValues(P,1,&gid1,1,&gid0,&one,INSERT_VALUES);
776:         } /* matched */
777:       } /* edge loop */

779:       /* deal with deleted ghost on first pass */
780:       if (size>1 && sub_it != n_sub_its-1) {
781:         PetscCDPos pos;  PetscBool ise = PETSC_FALSE;
782:         PetscInt   nSend1, **sbuffs1,nSend2;
783: #define REQ_BF_SIZE 100
784:         MPI_Request *sreqs2[REQ_BF_SIZE],*rreqs2[REQ_BF_SIZE];
785:         MPI_Status status;

787:         /* send request */
788:         for (proc=0,nSend1=0; proc<size; proc++) {
789:           PetscCDEmptyAt(deleted_list,proc,&ise);
790:           if (!ise) nSend1++;
791:         }
792:         PetscMalloc1(nSend1, &sbuffs1);
793:         /* PetscMalloc4(nSend1, sbuffs1, nSend1, rbuffs1, nSend1, sreqs1, nSend1, rreqs1); */
794:         /* PetscFree4(sbuffs1,rbuffs1,sreqs1,rreqs1); */
795:         for (proc=0,nSend1=0; proc<size; proc++) {
796:           /* count ghosts */
797:           PetscCDSizeAt(deleted_list,proc,&n);
798:           if (n>0) {
799: #define CHUNCK_SIZE 100
800:             PetscInt    *sbuff,*pt;
801:             MPI_Request *request;
802:             n   /= 2;
803:             PetscMalloc1(2 + 2*n + n*CHUNCK_SIZE + 2, &sbuff);
804:             /* PetscMalloc4(2+2*n,sbuffs1[nSend1],n*CHUNCK_SIZE,rbuffs1[nSend1],1,rreqs2[nSend1],1,sreqs2[nSend1]); */
805:             /* save requests */
806:             sbuffs1[nSend1] = sbuff;
807:             request         = (MPI_Request*)sbuff;
808:             sbuff           = pt = (PetscInt*)(request+1);
809:             *pt++           = n; *pt++ = rank;

811:             PetscCDGetHeadPos(deleted_list,proc,&pos);
812:             while (pos) {
813:               PetscInt lid0, cpid, gid;
814:               PetscLLNGetID(pos, &cpid);
815:               gid   = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
816:               PetscCDGetNextPos(deleted_list,proc,&pos);
817:               PetscLLNGetID(pos, &lid0);
818:               PetscCDGetNextPos(deleted_list,proc,&pos);
819:               *pt++ = gid; *pt++ = lid0;
820:             }
821:             /* send request tag1 [n, proc, n*[gid1,lid0] ] */
822:             MPI_Isend(sbuff, 2*n+2, MPIU_INT, proc, tag1, comm, request);
823:             /* post recieve */
824:             request        = (MPI_Request*)pt;
825:             rreqs2[nSend1] = request; /* cache recv request */
826:             pt             = (PetscInt*)(request+1);
827:             MPI_Irecv(pt, n*CHUNCK_SIZE, MPIU_INT, proc, tag2, comm, request);
828:             /* clear list */
829:             PetscCDRemoveAll(deleted_list, proc);
830:             nSend1++;
831:           }
832:         }
833:         /* recieve requests, send response, clear lists */
834:         kk     = nactive_edges;
835:         MPI_Allreduce(&kk,&nactive_edges,1,MPIU_INT,MPI_SUM,comm); /* not correct syncronization and global */
836:         nSend2 = 0;
837:         while (1) {
838: #define BF_SZ 10000
839:           PetscMPIInt flag,count;
840:           PetscInt rbuff[BF_SZ],*pt,*pt2,*pt3,count2,*sbuff,count3;
841:           MPI_Request *request;
842:           MPI_Iprobe(MPI_ANY_SOURCE, tag1, comm, &flag, &status);
843:           if (!flag) break;
844:           MPI_Get_count(&status, MPIU_INT, &count);
845:           if (count > BF_SZ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"buffer too small for recieve: %d",count);
846:           proc = status.MPI_SOURCE;
847:           /* recieve request tag1 [n, proc, n*[gid1,lid0] ] */
848:           MPI_Recv(rbuff, count, MPIU_INT, proc, tag1, comm, &status);
849:           /* count sends */
850:           pt = rbuff; count3 = count2 = 0;
851:           n  = *pt++; kk = *pt++;
852:           while (n--) {
853:             PetscInt gid1=*pt++, lid1=gid1-my0; kk=*pt++;
854:             if (lid_matched[lid1]) {
855:               PetscPrintf(PETSC_COMM_SELF,"\t *** [%d]%s %d) ERROR recieved deleted gid %d, deleted by (lid) %d from proc %d\n",rank,__FUNCT__,sub_it,gid1,kk);
856:               PetscSleep(1);
857:             }
858:             lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
859:             PetscCDSizeAt(agg_llists, lid1, &kk);
860:             count2           += kk + 2;
861:             count3++; /* number of verts requested (n) */
862:           }
863:           if (count2 > count3*CHUNCK_SIZE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Irecv will be too small: %d",count2);
864:           /* send tag2 *[lid0, n, n*[gid] ] */
865:           PetscMalloc(count2*sizeof(PetscInt) + sizeof(MPI_Request), &sbuff);
866:           request          = (MPI_Request*)sbuff;
867:           sreqs2[nSend2++] = request; /* cache request */
868:           if (nSend2==REQ_BF_SIZE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"buffer too small for requests: %d",nSend2);
869:           pt2 = sbuff = (PetscInt*)(request+1);
870:           pt  = rbuff;
871:           n   = *pt++; kk = *pt++;
872:           while (n--) {
873:             /* read [n, proc, n*[gid1,lid0] */
874:             PetscInt gid1=*pt++, lid1=gid1-my0, lid0=*pt++;
875:             /* write [lid0, n, n*[gid] ] */
876:             *pt2++ = lid0;
877:             pt3    = pt2++; /* save pointer for later */
878:             /* for (pos=PetscCDGetHeadPos(agg_llists,lid1) ; pos ; pos=PetscCDGetNextPos(agg_llists,lid1,pos)) { */
879:             PetscCDGetHeadPos(agg_llists,lid1,&pos);
880:             while (pos) {
881:               PetscInt gid;
882:               PetscLLNGetID(pos, &gid);
883:               PetscCDGetNextPos(agg_llists,lid1,&pos);
884:               *pt2++ = gid;
885:             }
886:             *pt3 = (pt2-pt3)-1;
887:             /* clear list */
888:             PetscCDRemoveAll(agg_llists, lid1);
889:           }
890:           /* send requested data tag2 *[lid0, n, n*[gid1] ] */
891:           MPI_Isend(sbuff, count2, MPIU_INT, proc, tag2, comm, request);
892:         }

894:         /* recieve tag2 *[lid0, n, n*[gid] ] */
895:         for (kk=0; kk<nSend1; kk++) {
896:           PetscMPIInt count;
897:           MPI_Request *request;
898:           PetscInt    *pt, *pt2;
899:           request = rreqs2[kk]; /* no need to free -- buffer is in 'sbuffs1' */
900:           MPI_Wait(request, &status);
901:           MPI_Get_count(&status, MPIU_INT, &count);
902:           pt      = pt2 = (PetscInt*)(request+1);
903:           while (pt-pt2 < count) {
904:             PetscInt lid0 = *pt++, n = *pt++;
905:             while (n--) {
906:               PetscInt gid1 = *pt++;
907:               PetscCDAppendID(agg_llists, lid0, gid1);
908:             }
909:           }
910:         }

912:         /* wait for tag1 isends */
913:         while (nSend1--) {
914:           MPI_Request *request;
915:           request = (MPI_Request*)sbuffs1[nSend1];
916:           MPI_Wait(request, &status);
917:           PetscFree(request);
918:         }
919:         PetscFree(sbuffs1);

921:         /* wait for tag2 isends */
922:         while (nSend2--) {
923:           MPI_Request *request = sreqs2[nSend2];
924:           MPI_Wait(request, &status);
925:           PetscFree(request);
926:         }

928:         VecRestoreArray(ghostMaxEdge, &cpcol_max_ew);
929:         VecRestoreArray(ghostMaxPE, &cpcol_max_pe);

931:         /* get 'cpcol_matched' - use locMaxPE, ghostMaxEdge, cpcol_max_ew */
932:         for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
933:           PetscScalar vval = lid_matched[kk] ? 1.0 : 0.0;
934:           VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
935:         }
936:         VecAssemblyBegin(locMaxPE);
937:         VecAssemblyEnd(locMaxPE);
938:         VecScatterBegin(mpimat->Mvctx,locMaxPE,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
939:           VecScatterEnd(mpimat->Mvctx,locMaxPE,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
940:         VecGetArray(ghostMaxEdge, &cpcol_max_ew);
941:         VecGetLocalSize(mpimat->lvec, &n);
942:         for (kk=0; kk<n; kk++) {
943:           cpcol_matched[kk] = (PetscBool)(PetscRealPart(cpcol_max_ew[kk]) != 0.0);
944:         }

946:         VecRestoreArray(ghostMaxEdge, &cpcol_max_ew);
947:       } /* size > 1 */

949:       /* compute 'locMaxEdge' */
950:       VecRestoreArray(locMaxEdge, &lid_max_ew);
951:       for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
952:         PetscReal max_e = 0.,tt;
953:         PetscScalar vval;
954:         PetscInt lid = kk;
955:         if (lid_matched[lid]) vval = 0.;
956:         else {
957:           ii = matA->i; n = ii[lid+1] - ii[lid]; idx = matA->j + ii[lid];
958:           ap = matA->a + ii[lid];
959:           for (jj=0; jj<n; jj++) {
960:             PetscInt lidj = idx[jj];
961:             if (lid_matched[lidj]) continue; /* this is new - can change local max */
962:             if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
963:           }
964:           if (lid_cprowID && (ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
965:             ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
966:             ap  = matB->a + ii[ix];
967:             idx = matB->j + ii[ix];
968:             for (jj=0; jj<n; jj++) {
969:               PetscInt lidj = idx[jj];
970:               if (cpcol_matched[lidj]) continue;
971:               if ((tt=PetscRealPart(ap[jj])) > max_e) max_e = tt;
972:             }
973:           }
974:         }
975:         vval = (PetscScalar)max_e;
976:         VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
977:       }
978:       VecAssemblyBegin(locMaxEdge);
979:       VecAssemblyEnd(locMaxEdge);

981:       if (size>1 && sub_it != n_sub_its-1) {
982:         /* compute 'cpcol_max_ew' */
983:         VecScatterBegin(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
984:           VecScatterEnd(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
985:         VecGetArray(ghostMaxEdge, &cpcol_max_ew);
986:         VecGetArray(locMaxEdge, &lid_max_ew);

988:         /* compute 'cpcol_max_pe' */
989:         for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
990:           PetscInt lid = kk;
991:           PetscReal ew,v1_max_e,v0_max_e=PetscRealPart(lid_max_ew[lid]);
992:           PetscScalar vval;
993:           PetscMPIInt max_pe=rank,pe;
994:           if (lid_matched[lid]) vval = (PetscScalar)rank;
995:           else if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
996:             ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
997:             ap  = matB->a + ii[ix];
998:             idx = matB->j + ii[ix];
999:             for (jj=0; jj<n; jj++) {
1000:               PetscInt lidj = idx[jj];
1001:               if (cpcol_matched[lidj]) continue;
1002:               ew = PetscRealPart(ap[jj]); v1_max_e = PetscRealPart(cpcol_max_ew[lidj]);
1003:               /* get max pe that has a max_e == to this edge w */
1004:               if ((pe=cpcol_pe[idx[jj]]) > max_pe && ew > v1_max_e - 1.e-12 && ew > v0_max_e - 1.e-12) max_pe = pe;
1005:             }
1006:             vval = (PetscScalar)max_pe;
1007:           }
1008:           VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES);
1009:         }
1010:         VecAssemblyBegin(locMaxPE);
1011:         VecAssemblyEnd(locMaxPE);

1013:         VecScatterBegin(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
1014:           VecScatterEnd(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
1015:         VecGetArray(ghostMaxPE, &cpcol_max_pe);
1016:         VecRestoreArray(locMaxEdge, &lid_max_ew);
1017:       } /* deal with deleted ghost */
1018:       PetscInfo3(a_Gmat,"\t %D.%D: %D active edges.\n",iter,sub_it,nactive_edges);
1019:       if (!nactive_edges) break;
1020:     } /* sub_it loop */

1022:     /* clean up iteration */
1023:     PetscFree(Edges);
1024:     if (mpimat) {
1025:       VecRestoreArray(ghostMaxEdge, &cpcol_max_ew);
1026:       VecDestroy(&ghostMaxEdge);
1027:       VecRestoreArray(ghostMaxPE, &cpcol_max_pe);
1028:       VecDestroy(&ghostMaxPE);
1029:       PetscFree(cpcol_pe);
1030:       PetscFree(cpcol_matched);
1031:     }

1033:     VecDestroy(&locMaxEdge);
1034:     VecDestroy(&locMaxPE);

1036:     if (mpimat) {
1037:       VecRestoreArray(mpimat->lvec, &cpcol_gid);
1038:     }

1040:     /* create next G if needed */
1041:     if (iter == n_iter) { /* hard wired test - need to look at full surrounded nodes or something */
1042:       MatDestroy(&P);
1043:       MatDestroy(&cMat);
1044:       break;
1045:     } else {
1046:       Vec diag;
1047:       /* add identity for unmatched vertices so they stay alive */
1048:       for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
1049:         if (!lid_matched[kk]) {
1050:           gid  = kk+my0;
1051:           MatGetRow(cMat,gid,&n,0,0);
1052:           if (n>1) {
1053:             MatSetValues(P,1,&gid,1,&gid,&one,INSERT_VALUES);
1054:           }
1055:           MatRestoreRow(cMat,gid,&n,0,0);
1056:         }
1057:       }
1058:       MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
1059:       MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);

1061:       /* project to make new graph with colapsed edges */
1062:       MatPtAP(cMat,P,MAT_INITIAL_MATRIX,1.0,&tMat);
1063:       MatDestroy(&P);
1064:       MatDestroy(&cMat);
1065:       cMat = tMat;
1066:       MatCreateVecs(cMat, &diag, 0);
1067:       MatGetDiagonal(cMat, diag); /* effectively PCJACOBI */
1068:       VecReciprocal(diag);
1069:       VecSqrtAbs(diag);
1070:       MatDiagonalScale(cMat, diag, diag);
1071:       VecDestroy(&diag);
1072:     }
1073:   } /* coarsen iterator */

1075:   /* make fake matrix */
1076:   if (size>1) {
1077:     Mat mat;
1078:     PetscCDPos pos;
1079:     PetscInt gid, NN, MM, jj = 0, mxsz = 0;

1081:     for (kk=0; kk<nloc; kk++) {
1082:       PetscCDSizeAt(agg_llists, kk, &jj);
1083:       if (jj > mxsz) mxsz = jj;
1084:     }
1085:     MatGetSize(a_Gmat, &MM, &NN);
1086:     if (mxsz > MM-nloc) mxsz = MM-nloc;

1088:     MatCreateAIJ(comm, nloc, nloc,PETSC_DETERMINE, PETSC_DETERMINE,0, 0, mxsz, 0, &mat);

1090:     /* */
1091:     for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
1092:       /* for (pos=PetscCDGetHeadPos(agg_llists,kk) ; pos ; pos=PetscCDGetNextPos(agg_llists,kk,pos)) { */
1093:       PetscCDGetHeadPos(agg_llists,kk,&pos);
1094:       while (pos) {
1095:         PetscInt gid1;
1096:         PetscLLNGetID(pos, &gid1);
1097:         PetscCDGetNextPos(agg_llists,kk,&pos);

1099:         if (gid1 < my0 || gid1 >= my0+nloc) {
1100:           MatSetValues(mat,1,&gid,1,&gid1,&one,ADD_VALUES);
1101:         }
1102:       }
1103:     }
1104:     MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1105:     MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

1107:     PetscCDSetMat(agg_llists, mat);
1108:   }

1110:   PetscFree(lid_cprowID);
1111:   PetscFree(lid_gid);
1112:   PetscFree(lid_matched);
1113:   PetscCDDestroy(deleted_list);
1114:   return(0);
1115: }

1117: typedef struct {
1118:   int dummy;
1119: } MatCoarsen_HEM;
1120: /*
1121:    HEM coarsen, simple greedy.
1122: */
1125: static PetscErrorCode MatCoarsenApply_HEM(MatCoarsen coarse)
1126: {
1127:   /* MatCoarsen_HEM *HEM = (MatCoarsen_HEM*)coarse->subctx; */
1129:   Mat            mat = coarse->graph;

1133:   if (!coarse->perm) {
1134:     IS       perm;
1135:     PetscInt n,m;
1136: 
1137:     MatGetLocalSize(mat, &m, &n);
1138:     ISCreateStride(PetscObjectComm((PetscObject)mat), m, 0, 1, &perm);
1139:     heavyEdgeMatchAgg(perm, mat, &coarse->agg_lists);
1140:     ISDestroy(&perm);
1141:   } else {
1142:     heavyEdgeMatchAgg(coarse->perm, mat, &coarse->agg_lists);
1143:   }
1144:   return(0);
1145: }

1149: static PetscErrorCode MatCoarsenView_HEM(MatCoarsen coarse,PetscViewer viewer)
1150: {
1151:   /* MatCoarsen_HEM *HEM = (MatCoarsen_HEM*)coarse->subctx; */
1153:   PetscMPIInt    rank;
1154:   PetscBool      iascii;

1158:   MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
1159:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1160:   if (iascii) {
1161:     PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] HEM aggregator\n",rank);
1162:     PetscViewerFlush(viewer);
1163:     PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
1164:   }
1165:   return(0);
1166: }

1170: static PetscErrorCode MatCoarsenDestroy_HEM(MatCoarsen coarse)
1171: {
1172:   MatCoarsen_HEM *HEM = (MatCoarsen_HEM*)coarse->subctx;

1177:   PetscFree(HEM);
1178:   return(0);
1179: }

1181: /*MC
1182:    MATCOARSENHEM - A coarsener that uses HEM a simple greedy coarsener

1184:    Level: beginner

1186: .keywords: Coarsen, create, context

1188: .seealso: MatCoarsenSetType(), MatCoarsenType, MatCoarsenCreate()

1190: M*/

1194: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_HEM(MatCoarsen coarse)
1195: {
1197:   MatCoarsen_HEM *HEM;

1200:   PetscNewLog(coarse,&HEM);
1201:   coarse->subctx       = (void*)HEM;
1202:   coarse->ops->apply   = MatCoarsenApply_HEM;
1203:   coarse->ops->view    = MatCoarsenView_HEM;
1204:   coarse->ops->destroy = MatCoarsenDestroy_HEM;
1205:   /* coarse->ops->setfromoptions = MatCoarsenSetFromOptions_HEM; */
1206:   return(0);
1207: }