Actual source code: hem.c
petsc-3.4.5 2014-06-29
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: PetscMalloc(sizeof(PetscCoarsenData), &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: PetscMalloc(a_size*sizeof(PetscCDIntNd*), &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: PetscMalloc(ail->chk_sz*sizeof(PetscCDIntNd), &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: PetscMalloc(kk*sizeof(PetscInt), &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: PetscMalloc(kk*sizeof(IS*), &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: PetscMalloc(a_bs*lsz*sizeof(PetscInt), &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: 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)
452: . verbose -
453: Output Parameter:
454: . a_locals_llist - array of list of local nodes rooted at local node
455: */
458: PetscErrorCode heavyEdgeMatchAgg(IS perm,Mat a_Gmat,PetscInt verbose,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: PetscMalloc(nloc*sizeof(PetscInt), &lid_gid); /* explicit array needed */
485: PetscMalloc(nloc*sizeof(PetscInt), &lid_cprowID);
486: PetscMalloc(nloc*sizeof(PetscBool), &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: /* force compressed storage of B */
519: matB->compressedrow.check = PETSC_TRUE;
521: MatCheckCompressedRow(mpimat->B,&matB->compressedrow,matB->i,cMat->rmap->n,-1.0);
522: if (!matB->compressedrow.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"matB must have compressed row usage");
523: } else {
524: matA = (Mat_SeqAIJ*)cMat->data;
525: }
527: /* set max edge on nodes */
528: MatGetVecs(cMat, &locMaxEdge, 0);
529: MatGetVecs(cMat, &locMaxPE, 0);
531: /* get 'cpcol_pe' & 'cpcol_gid' & init. 'cpcol_matched' using 'mpimat->lvec' */
532: if (mpimat) {
533: Vec vec;
534: PetscScalar vval;
536: MatGetVecs(cMat, &vec, 0);
537: /* cpcol_pe */
538: vval = (PetscScalar)(rank);
539: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
540: VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
541: }
542: VecAssemblyBegin(vec);
543: VecAssemblyEnd(vec);
544: VecScatterBegin(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
545: VecScatterEnd(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
546: VecGetArray(mpimat->lvec, &cpcol_gid); /* get proc ID in 'cpcol_gid' */
547: VecGetLocalSize(mpimat->lvec, &n);
548: PetscMalloc(n*sizeof(PetscInt), &cpcol_pe);
549: for (kk=0; kk<n; kk++) cpcol_pe[kk] = (PetscMPIInt)PetscRealPart(cpcol_gid[kk]);
550: VecRestoreArray(mpimat->lvec, &cpcol_gid);
552: /* cpcol_gid */
553: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
554: vval = (PetscScalar)(gid);
555: VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
556: }
557: VecAssemblyBegin(vec);
558: VecAssemblyEnd(vec);
559: VecScatterBegin(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
560: VecScatterEnd(mpimat->Mvctx,vec,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
561: VecDestroy(&vec);
562: VecGetArray(mpimat->lvec, &cpcol_gid); /* get proc ID in 'cpcol_gid' */
564: /* cpcol_matched */
565: VecGetLocalSize(mpimat->lvec, &n);
566: PetscMalloc(n*sizeof(PetscBool), &cpcol_matched);
567: for (kk=0; kk<n; kk++) cpcol_matched[kk] = PETSC_FALSE;
568: }
570: /* need an inverse map - locals */
571: for (kk=0; kk<nloc; kk++) lid_cprowID[kk] = -1;
572: /* set index into compressed row 'lid_cprowID' */
573: if (matB) {
574: ii = matB->compressedrow.i;
575: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
576: lid_cprowID[matB->compressedrow.rindex[ix]] = ix;
577: }
578: }
580: /* get removed IS, use '' */
581: /* if (iter==1) { */
582: /* PetscInt *lid_rem,idx; */
583: /* PetscMalloc(nloc*sizeof(PetscInt), &lid_rem); */
584: /* for (kk=idx=0;kk<nloc;kk++) { */
585: /* PetscInt nn,lid=kk; */
586: /* ii = matA->i; nn = ii[lid+1] - ii[lid]; */
587: /* if ((ix=lid_cprowID[lid]) != -1) { /\* if I have any ghost neighbors *\/ */
588: /* ii = matB->compressedrow.i; */
589: /* nn += ii[ix+1] - ii[ix]; */
590: /* } */
591: /* if (nn < 2) { */
592: /* lid_rem[idx++] = kk + my0; */
593: /* } */
594: /* } */
595: /* PetscCDSetRemovedIS(agg_llists, comm, idx, lid_rem); */
596: /* PetscFree(lid_rem); */
597: /* } */
599: /* compute 'locMaxEdge' & 'locMaxPE', and create list of edges, count edges' */
600: for (nEdges=0,kk=0,gid=my0; kk<nloc; kk++,gid++) {
601: PetscReal max_e = 0., tt;
602: PetscScalar vval;
603: PetscInt lid = kk;
604: PetscMPIInt max_pe=rank,pe;
605: ii = matA->i; n = ii[lid+1] - ii[lid]; idx = matA->j + ii[lid];
606: ap = matA->a + ii[lid];
607: for (jj=0; jj<n; jj++) {
608: PetscInt lidj = idx[jj];
609: if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
610: if (lidj > lid) nEdges++;
611: }
612: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
613: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
614: ap = matB->a + ii[ix];
615: idx = matB->j + ii[ix];
616: for (jj=0; jj<n; jj++) {
617: if ((tt=PetscRealPart(ap[jj])) > max_e) max_e = tt;
618: nEdges++;
619: if ((pe=cpcol_pe[idx[jj]]) > max_pe) max_pe = pe;
620: }
621: }
622: vval = max_e;
623: VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES);
625: vval = (PetscScalar)max_pe;
626: VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES);
627: }
628: VecAssemblyBegin(locMaxEdge);
629: VecAssemblyEnd(locMaxEdge);
630: VecAssemblyBegin(locMaxPE);
631: VecAssemblyEnd(locMaxPE);
633: /* get 'cpcol_max_ew' & 'cpcol_max_pe' */
634: if (mpimat) {
635: VecDuplicate(mpimat->lvec, &ghostMaxEdge);
636: VecScatterBegin(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
637: VecScatterEnd(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
638: VecGetArray(ghostMaxEdge, &cpcol_max_ew);
640: VecDuplicate(mpimat->lvec, &ghostMaxPE);
641: VecScatterBegin(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
642: VecScatterEnd(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
643: VecGetArray(ghostMaxPE, &cpcol_max_pe);
644: }
646: /* setup sorted list of edges */
647: PetscMalloc(nEdges*sizeof(Edge), &Edges);
648: ISGetIndices(perm, &perm_ix);
649: for (nEdges=n_nz_row=kk=0; kk<nloc; kk++) {
650: PetscInt nn, lid = perm_ix[kk];
651: ii = matA->i; nn = n = ii[lid+1] - ii[lid]; idx = matA->j + ii[lid];
652: ap = matA->a + ii[lid];
653: for (jj=0; jj<n; jj++) {
654: PetscInt lidj = idx[jj];
655: if (lidj > lid) {
656: Edges[nEdges].lid0 = lid;
657: Edges[nEdges].gid1 = lidj + my0;
658: Edges[nEdges].cpid1 = -1;
659: Edges[nEdges].weight = PetscRealPart(ap[jj]);
660: nEdges++;
661: }
662: }
663: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
664: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
665: ap = matB->a + ii[ix];
666: idx = matB->j + ii[ix];
667: nn += n;
668: for (jj=0; jj<n; jj++) {
669: Edges[nEdges].lid0 = lid;
670: Edges[nEdges].gid1 = (PetscInt)PetscRealPart(cpcol_gid[idx[jj]]);
671: Edges[nEdges].cpid1 = idx[jj];
672: Edges[nEdges].weight = PetscRealPart(ap[jj]);
673: nEdges++;
674: }
675: }
676: if (nn > 1) n_nz_row++;
677: else if (iter == 1) {
678: /* should select this because it is technically in the MIS but lets not */
679: PetscCDRemoveAll(agg_llists, lid);
680: }
681: }
682: ISRestoreIndices(perm,&perm_ix);
684: qsort(Edges, nEdges, sizeof(Edge), gamg_hem_compare);
686: /* projection matrix */
687: MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 1, 0, 1, 0, &P);
689: /* clear matched flags */
690: for (kk=0; kk<nloc; kk++) lid_matched[kk] = PETSC_FALSE;
691: /* process - communicate - process */
692: for (sub_it=0; sub_it<n_sub_its; sub_it++) {
693: PetscInt nactive_edges;
695: VecGetArray(locMaxEdge, &lid_max_ew);
696: for (kk=nactive_edges=0; kk<nEdges; kk++) {
697: /* HEM */
698: const Edge *e = &Edges[kk];
699: const PetscInt lid0 =e->lid0,gid1=e->gid1,cpid1=e->cpid1,gid0=lid0+my0,lid1=gid1-my0;
700: PetscBool isOK = PETSC_TRUE;
702: /* skip if either (local) vertex is done already */
703: if (lid_matched[lid0] || (gid1>=my0 && gid1<Iend && lid_matched[gid1-my0])) continue;
705: /* skip if ghost vertex is done */
706: if (cpid1 != -1 && cpcol_matched[cpid1]) continue;
708: nactive_edges++;
709: /* skip if I have a bigger edge someplace (lid_max_ew gets updated) */
710: if (PetscRealPart(lid_max_ew[lid0]) > e->weight + 1.e-12) continue;
712: if (cpid1 == -1) {
713: if (PetscRealPart(lid_max_ew[lid1]) > e->weight + 1.e-12) continue;
714: } else {
715: /* see if edge might get matched on other proc */
716: PetscReal g_max_e = PetscRealPart(cpcol_max_ew[cpid1]);
717: if (g_max_e > e->weight + 1.e-12) continue;
718: else if (e->weight > g_max_e - 1.e-12 && (PetscMPIInt)PetscRealPart(cpcol_max_pe[cpid1]) > rank) {
719: /* check for max_e == to this edge and larger processor that will deal with this */
720: continue;
721: }
722: }
724: /* check ghost for v0 */
725: if (isOK) {
726: PetscReal max_e,ew;
727: if ((ix=lid_cprowID[lid0]) != -1) { /* if I have any ghost neighbors */
728: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
729: ap = matB->a + ii[ix];
730: idx = matB->j + ii[ix];
731: for (jj=0; jj<n && isOK; jj++) {
732: PetscInt lidj = idx[jj];
733: if (cpcol_matched[lidj]) continue;
734: ew = PetscRealPart(ap[jj]); max_e = PetscRealPart(cpcol_max_ew[lidj]);
735: /* check for max_e == to this edge and larger processor that will deal with this */
736: if (ew > max_e - 1.e-12 && ew > PetscRealPart(lid_max_ew[lid0]) - 1.e-12 && (PetscMPIInt)PetscRealPart(cpcol_max_pe[lidj]) > rank) {
737: isOK = PETSC_FALSE;
738: }
739: }
740: }
742: /* for v1 */
743: if (cpid1 == -1 && isOK) {
744: if ((ix=lid_cprowID[lid1]) != -1) { /* if I have any ghost neighbors */
745: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
746: ap = matB->a + ii[ix];
747: idx = matB->j + ii[ix];
748: for (jj=0; jj<n && isOK; jj++) {
749: PetscInt lidj = idx[jj];
750: if (cpcol_matched[lidj]) continue;
751: ew = PetscRealPart(ap[jj]); max_e = PetscRealPart(cpcol_max_ew[lidj]);
752: /* check for max_e == to this edge and larger processor that will deal with this */
753: if (ew > max_e - 1.e-12 && ew > PetscRealPart(lid_max_ew[lid1]) - 1.e-12 && (PetscMPIInt)PetscRealPart(cpcol_max_pe[lidj]) > rank) {
754: isOK = PETSC_FALSE;
755: }
756: }
757: }
758: }
759: }
761: /* do it */
762: if (isOK) {
763: if (cpid1 == -1) {
764: lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
765: PetscCDAppendRemove(agg_llists, lid0, lid1);
766: } else if (sub_it != n_sub_its-1) {
767: /* add gid1 to list of ghost deleted by me -- I need their children */
768: proc = cpcol_pe[cpid1];
770: cpcol_matched[cpid1] = PETSC_TRUE; /* messing with VecGetArray array -- needed??? */
772: PetscCDAppendID(deleted_list, proc, cpid1); /* cache to send messages */
773: PetscCDAppendID(deleted_list, proc, lid0);
774: } else continue;
776: lid_matched[lid0] = PETSC_TRUE; /* keep track of what we've done this round */
777: /* set projection */
778: MatSetValues(P,1,&gid0,1,&gid0,&one,INSERT_VALUES);
779: MatSetValues(P,1,&gid1,1,&gid0,&one,INSERT_VALUES);
780: } /* matched */
781: } /* edge loop */
783: /* deal with deleted ghost on first pass */
784: if (size>1 && sub_it != n_sub_its-1) {
785: PetscCDPos pos; PetscBool ise = PETSC_FALSE;
786: PetscInt nSend1, **sbuffs1,nSend2;
787: #define REQ_BF_SIZE 100
788: MPI_Request *sreqs2[REQ_BF_SIZE],*rreqs2[REQ_BF_SIZE];
789: MPI_Status status;
791: /* send request */
792: for (proc=0,nSend1=0; proc<size; proc++) {
793: PetscCDEmptyAt(deleted_list,proc,&ise);
794: if (!ise) nSend1++;
795: }
796: PetscMalloc(nSend1*sizeof(PetscInt*), &sbuffs1);
797: /* PetscMalloc4(nSend1, PetscInt*, sbuffs1, nSend1, PetscInt*, rbuffs1, nSend1, MPI_Request*, sreqs1, nSend1, MPI_Request*, rreqs1); */
798: /* PetscFree4(sbuffs1,rbuffs1,sreqs1,rreqs1); */
799: for (proc=0,nSend1=0; proc<size; proc++) {
800: /* count ghosts */
801: PetscCDSizeAt(deleted_list,proc,&n);
802: if (n>0) {
803: #define CHUNCK_SIZE 100
804: PetscInt *sbuff,*pt;
805: MPI_Request *request;
806: n /= 2;
807: PetscMalloc((2 + 2*n + n*CHUNCK_SIZE)*sizeof(PetscInt) + 2*sizeof(MPI_Request), &sbuff);
808: /* PetscMalloc4(2+2*n,PetscInt,sbuffs1[nSend1],n*CHUNCK_SIZE,PetscInt,rbuffs1[nSend1],1,MPI_Request,rreqs2[nSend1],1,MPI_Request,sreqs2[nSend1]); */
809: /* save requests */
810: sbuffs1[nSend1] = sbuff;
811: request = (MPI_Request*)sbuff;
812: sbuff = pt = (PetscInt*)(request+1);
813: *pt++ = n; *pt++ = rank;
815: PetscCDGetHeadPos(deleted_list,proc,&pos);
816: while (pos) {
817: PetscInt lid0, cpid, gid;
818: PetscLLNGetID(pos, &cpid);
819: gid = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
820: PetscCDGetNextPos(deleted_list,proc,&pos);
821: PetscLLNGetID(pos, &lid0);
822: PetscCDGetNextPos(deleted_list,proc,&pos);
823: *pt++ = gid; *pt++ = lid0;
824: }
825: /* send request tag1 [n, proc, n*[gid1,lid0] ] */
826: MPI_Isend(sbuff, 2*n+2, MPIU_INT, proc, tag1, comm, request);
827: /* post recieve */
828: request = (MPI_Request*)pt;
829: rreqs2[nSend1] = request; /* cache recv request */
830: pt = (PetscInt*)(request+1);
831: MPI_Irecv(pt, n*CHUNCK_SIZE, MPIU_INT, proc, tag2, comm, request);
832: /* clear list */
833: PetscCDRemoveAll(deleted_list, proc);
834: nSend1++;
835: }
836: }
837: /* recieve requests, send response, clear lists */
838: kk = nactive_edges;
839: MPI_Allreduce(&kk,&nactive_edges,1,MPIU_INT,MPI_SUM,comm); /* not correct syncronization and global */
840: nSend2 = 0;
841: while (1) {
842: #define BF_SZ 10000
843: PetscMPIInt flag,count;
844: PetscInt rbuff[BF_SZ],*pt,*pt2,*pt3,count2,*sbuff,count3;
845: MPI_Request *request;
846: MPI_Iprobe(MPI_ANY_SOURCE, tag1, comm, &flag, &status);
847: if (!flag) break;
848: MPI_Get_count(&status, MPIU_INT, &count);
849: if (count > BF_SZ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"buffer too small for recieve: %d",count);
850: proc = status.MPI_SOURCE;
851: /* recieve request tag1 [n, proc, n*[gid1,lid0] ] */
852: MPI_Recv(rbuff, count, MPIU_INT, proc, tag1, comm, &status);
853: /* count sends */
854: pt = rbuff; count3 = count2 = 0;
855: n = *pt++; kk = *pt++;
856: while (n--) {
857: PetscInt gid1=*pt++, lid1=gid1-my0; kk=*pt++;
858: if (lid_matched[lid1]) {
859: 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);
860: PetscSleep(1);
861: }
862: lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
863: PetscCDSizeAt(agg_llists, lid1, &kk);
864: count2 += kk + 2;
865: count3++; /* number of verts requested (n) */
866: }
867: if (count2 > count3*CHUNCK_SIZE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Irecv will be too small: %d",count2);
868: /* send tag2 *[lid0, n, n*[gid] ] */
869: PetscMalloc(count2*sizeof(PetscInt) + sizeof(MPI_Request), &sbuff);
870: request = (MPI_Request*)sbuff;
871: sreqs2[nSend2++] = request; /* cache request */
872: if (nSend2==REQ_BF_SIZE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"buffer too small for requests: %d",nSend2);
873: pt2 = sbuff = (PetscInt*)(request+1);
874: pt = rbuff;
875: n = *pt++; kk = *pt++;
876: while (n--) {
877: /* read [n, proc, n*[gid1,lid0] */
878: PetscInt gid1=*pt++, lid1=gid1-my0, lid0=*pt++;
879: /* write [lid0, n, n*[gid] ] */
880: *pt2++ = lid0;
881: pt3 = pt2++; /* save pointer for later */
882: /* for (pos=PetscCDGetHeadPos(agg_llists,lid1) ; pos ; pos=PetscCDGetNextPos(agg_llists,lid1,pos)) { */
883: PetscCDGetHeadPos(agg_llists,lid1,&pos);
884: while (pos) {
885: PetscInt gid;
886: PetscLLNGetID(pos, &gid);
887: PetscCDGetNextPos(agg_llists,lid1,&pos);
888: *pt2++ = gid;
889: }
890: *pt3 = (pt2-pt3)-1;
891: /* clear list */
892: PetscCDRemoveAll(agg_llists, lid1);
893: }
894: /* send requested data tag2 *[lid0, n, n*[gid1] ] */
895: MPI_Isend(sbuff, count2, MPIU_INT, proc, tag2, comm, request);
896: }
898: /* recieve tag2 *[lid0, n, n*[gid] ] */
899: for (kk=0; kk<nSend1; kk++) {
900: PetscMPIInt count;
901: MPI_Request *request;
902: PetscInt *pt, *pt2;
903: request = rreqs2[kk]; /* no need to free -- buffer is in 'sbuffs1' */
904: MPI_Wait(request, &status);
905: MPI_Get_count(&status, MPIU_INT, &count);
906: pt = pt2 = (PetscInt*)(request+1);
907: while (pt-pt2 < count) {
908: PetscInt lid0 = *pt++, n = *pt++;
909: while (n--) {
910: PetscInt gid1 = *pt++;
911: PetscCDAppendID(agg_llists, lid0, gid1);
912: }
913: }
914: }
916: /* wait for tag1 isends */
917: while (nSend1--) {
918: MPI_Request *request;
919: request = (MPI_Request*)sbuffs1[nSend1];
920: MPI_Wait(request, &status);
921: PetscFree(request);
922: }
923: PetscFree(sbuffs1);
925: /* wait for tag2 isends */
926: while (nSend2--) {
927: MPI_Request *request = sreqs2[nSend2];
928: MPI_Wait(request, &status);
929: PetscFree(request);
930: }
932: VecRestoreArray(ghostMaxEdge, &cpcol_max_ew);
933: VecRestoreArray(ghostMaxPE, &cpcol_max_pe);
935: /* get 'cpcol_matched' - use locMaxPE, ghostMaxEdge, cpcol_max_ew */
936: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
937: PetscScalar vval = lid_matched[kk] ? 1.0 : 0.0;
938: VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
939: }
940: VecAssemblyBegin(locMaxPE);
941: VecAssemblyEnd(locMaxPE);
942: VecScatterBegin(mpimat->Mvctx,locMaxPE,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
943: VecScatterEnd(mpimat->Mvctx,locMaxPE,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
944: VecGetArray(ghostMaxEdge, &cpcol_max_ew);
945: VecGetLocalSize(mpimat->lvec, &n);
946: for (kk=0; kk<n; kk++) {
947: cpcol_matched[kk] = (PetscBool)(PetscRealPart(cpcol_max_ew[kk]) != 0.0);
948: }
950: VecRestoreArray(ghostMaxEdge, &cpcol_max_ew);
951: } /* size > 1 */
953: /* compute 'locMaxEdge' */
954: VecRestoreArray(locMaxEdge, &lid_max_ew);
955: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
956: PetscReal max_e = 0.,tt;
957: PetscScalar vval;
958: PetscInt lid = kk;
959: if (lid_matched[lid]) vval = 0.;
960: else {
961: ii = matA->i; n = ii[lid+1] - ii[lid]; idx = matA->j + ii[lid];
962: ap = matA->a + ii[lid];
963: for (jj=0; jj<n; jj++) {
964: PetscInt lidj = idx[jj];
965: if (lid_matched[lidj]) continue; /* this is new - can change local max */
966: if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
967: }
968: if (lid_cprowID && (ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
969: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
970: ap = matB->a + ii[ix];
971: idx = matB->j + ii[ix];
972: for (jj=0; jj<n; jj++) {
973: PetscInt lidj = idx[jj];
974: if (cpcol_matched[lidj]) continue;
975: if ((tt=PetscRealPart(ap[jj])) > max_e) max_e = tt;
976: }
977: }
978: }
979: vval = (PetscScalar)max_e;
980: VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES); /* set with GID */
981: }
982: VecAssemblyBegin(locMaxEdge);
983: VecAssemblyEnd(locMaxEdge);
985: if (size>1 && sub_it != n_sub_its-1) {
986: /* compute 'cpcol_max_ew' */
987: VecScatterBegin(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
988: VecScatterEnd(mpimat->Mvctx,locMaxEdge,ghostMaxEdge,INSERT_VALUES,SCATTER_FORWARD);
989: VecGetArray(ghostMaxEdge, &cpcol_max_ew);
990: VecGetArray(locMaxEdge, &lid_max_ew);
992: /* compute 'cpcol_max_pe' */
993: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
994: PetscInt lid = kk;
995: PetscReal ew,v1_max_e,v0_max_e=PetscRealPart(lid_max_ew[lid]);
996: PetscScalar vval;
997: PetscMPIInt max_pe=rank,pe;
998: if (lid_matched[lid]) vval = (PetscScalar)rank;
999: else if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
1000: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
1001: ap = matB->a + ii[ix];
1002: idx = matB->j + ii[ix];
1003: for (jj=0; jj<n; jj++) {
1004: PetscInt lidj = idx[jj];
1005: if (cpcol_matched[lidj]) continue;
1006: ew = PetscRealPart(ap[jj]); v1_max_e = PetscRealPart(cpcol_max_ew[lidj]);
1007: /* get max pe that has a max_e == to this edge w */
1008: 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;
1009: }
1010: vval = (PetscScalar)max_pe;
1011: }
1012: VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES);
1013: }
1014: VecAssemblyBegin(locMaxPE);
1015: VecAssemblyEnd(locMaxPE);
1017: VecScatterBegin(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
1018: VecScatterEnd(mpimat->Mvctx,locMaxPE,ghostMaxPE,INSERT_VALUES,SCATTER_FORWARD);
1019: VecGetArray(ghostMaxPE, &cpcol_max_pe);
1020: VecRestoreArray(locMaxEdge, &lid_max_ew);
1021: } /* deal with deleted ghost */
1022: if (verbose>2) PetscPrintf(comm,"\t[%d]%s %d.%d: %d active edges.\n",rank,__FUNCT__,iter,sub_it,nactive_edges);
1023: if (!nactive_edges) break;
1024: } /* sub_it loop */
1026: /* clean up iteration */
1027: PetscFree(Edges);
1028: if (mpimat) {
1029: VecRestoreArray(ghostMaxEdge, &cpcol_max_ew);
1030: VecDestroy(&ghostMaxEdge);
1031: VecRestoreArray(ghostMaxPE, &cpcol_max_pe);
1032: VecDestroy(&ghostMaxPE);
1033: PetscFree(cpcol_pe);
1034: PetscFree(cpcol_matched);
1035: }
1037: VecDestroy(&locMaxEdge);
1038: VecDestroy(&locMaxPE);
1040: if (mpimat) {
1041: VecRestoreArray(mpimat->lvec, &cpcol_gid);
1042: }
1044: /* create next G if needed */
1045: if (iter == n_iter) { /* hard wired test - need to look at full surrounded nodes or something */
1046: MatDestroy(&P);
1047: MatDestroy(&cMat);
1048: break;
1049: } else {
1050: Vec diag;
1051: /* add identity for unmatched vertices so they stay alive */
1052: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
1053: if (!lid_matched[kk]) {
1054: gid = kk+my0;
1055: MatGetRow(cMat,gid,&n,0,0);
1056: if (n>1) {
1057: MatSetValues(P,1,&gid,1,&gid,&one,INSERT_VALUES);
1058: }
1059: MatRestoreRow(cMat,gid,&n,0,0);
1060: }
1061: }
1062: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
1063: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
1065: /* project to make new graph with colapsed edges */
1066: MatPtAP(cMat,P,MAT_INITIAL_MATRIX,1.0,&tMat);
1067: MatDestroy(&P);
1068: MatDestroy(&cMat);
1069: cMat = tMat;
1070: MatGetVecs(cMat, &diag, 0);
1071: MatGetDiagonal(cMat, diag); /* effectively PCJACOBI */
1072: VecReciprocal(diag);
1073: VecSqrtAbs(diag);
1074: MatDiagonalScale(cMat, diag, diag);
1075: VecDestroy(&diag);
1076: }
1077: } /* coarsen iterator */
1079: /* make fake matrix */
1080: if (size>1) {
1081: Mat mat;
1082: PetscCDPos pos;
1083: PetscInt gid, NN, MM, jj = 0, mxsz = 0;
1085: for (kk=0; kk<nloc; kk++) {
1086: PetscCDSizeAt(agg_llists, kk, &jj);
1087: if (jj > mxsz) mxsz = jj;
1088: }
1089: MatGetSize(a_Gmat, &MM, &NN);
1090: if (mxsz > MM-nloc) mxsz = MM-nloc;
1092: MatCreateAIJ(comm, nloc, nloc,PETSC_DETERMINE, PETSC_DETERMINE,0, 0, mxsz, 0, &mat);
1094: /* */
1095: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
1096: /* for (pos=PetscCDGetHeadPos(agg_llists,kk) ; pos ; pos=PetscCDGetNextPos(agg_llists,kk,pos)) { */
1097: PetscCDGetHeadPos(agg_llists,kk,&pos);
1098: while (pos) {
1099: PetscInt gid1;
1100: PetscLLNGetID(pos, &gid1);
1101: PetscCDGetNextPos(agg_llists,kk,&pos);
1103: if (gid1 < my0 || gid1 >= my0+nloc) {
1104: MatSetValues(mat,1,&gid,1,&gid1,&one,ADD_VALUES);
1105: }
1106: }
1107: }
1108: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1109: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1111: PetscCDSetMat(agg_llists, mat);
1112: }
1114: PetscFree(lid_cprowID);
1115: PetscFree(lid_gid);
1116: PetscFree(lid_matched);
1117: PetscCDDestroy(deleted_list);
1118: return(0);
1119: }
1121: typedef struct {
1122: int dummy;
1123: } MatCoarsen_HEM;
1124: /*
1125: HEM coarsen, simple greedy.
1126: */
1129: static PetscErrorCode MatCoarsenApply_HEM(MatCoarsen coarse)
1130: {
1131: /* MatCoarsen_HEM *HEM = (MatCoarsen_HEM*)coarse->subctx; */
1133: Mat mat = coarse->graph;
1137: if (!coarse->perm) {
1138: IS perm;
1139: PetscInt n,m;
1140:
1141: MatGetLocalSize(mat, &m, &n);
1142: ISCreateStride(PetscObjectComm((PetscObject)mat), m, 0, 1, &perm);
1143: heavyEdgeMatchAgg(perm, mat, coarse->verbose, &coarse->agg_lists);
1144: ISDestroy(&perm);
1145: } else {
1146: heavyEdgeMatchAgg(coarse->perm, mat, coarse->verbose, &coarse->agg_lists);
1147: }
1148: return(0);
1149: }
1153: PetscErrorCode MatCoarsenView_HEM(MatCoarsen coarse,PetscViewer viewer)
1154: {
1155: /* MatCoarsen_HEM *HEM = (MatCoarsen_HEM*)coarse->subctx; */
1157: PetscMPIInt rank;
1158: PetscBool iascii;
1162: MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
1163: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1164: if (iascii) {
1165: PetscViewerASCIISynchronizedPrintf(viewer," [%d] HEM aggregator\n",rank);
1166: PetscViewerFlush(viewer);
1167: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
1168: }
1169: return(0);
1170: }
1174: PetscErrorCode MatCoarsenDestroy_HEM(MatCoarsen coarse)
1175: {
1176: MatCoarsen_HEM *HEM = (MatCoarsen_HEM*)coarse->subctx;
1181: PetscFree(HEM);
1182: return(0);
1183: }
1185: /*MC
1186: MATCOARSENHEM - Creates a coarsen context via the external package HEM.
1188: Collective on MPI_Comm
1190: Input Parameter:
1191: . coarse - the coarsen context
1193: Options Database Keys:
1194: + -mat_coarsen_HEM_xxx -
1196: Level: beginner
1198: .keywords: Coarsen, create, context
1200: .seealso: MatCoarsenSetType(), MatCoarsenType
1202: M*/
1206: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_HEM(MatCoarsen coarse)
1207: {
1209: MatCoarsen_HEM *HEM;
1212: PetscNewLog(coarse, MatCoarsen_HEM, &HEM);
1213: coarse->subctx = (void*)HEM;
1214: coarse->ops->apply = MatCoarsenApply_HEM;
1215: coarse->ops->view = MatCoarsenView_HEM;
1216: coarse->ops->destroy = MatCoarsenDestroy_HEM;
1217: /* coarse->ops->setfromoptions = MatCoarsenSetFromOptions_HEM; */
1218: return(0);
1219: }