Actual source code: geo.c
petsc-3.13.6 2020-09-29
1: /*
2: GAMG geometric-algebric multiogrid PC - Mark Adams 2011
3: */
5: #include <../src/ksp/pc/impls/gamg/gamg.h>
7: #if defined(PETSC_HAVE_TRIANGLE)
8: #if !defined(ANSI_DECLARATORS)
9: #define ANSI_DECLARATORS
10: #endif
11: #include <triangle.h>
12: #endif
14: #include <petscblaslapack.h>
16: /* Private context for the GAMG preconditioner */
17: typedef struct {
18: PetscInt lid; /* local vertex index */
19: PetscInt degree; /* vertex degree */
20: } GAMGNode;
22: PETSC_STATIC_INLINE int petsc_geo_mg_compare(const void *a, const void *b)
23: {
24: return (((GAMGNode*)a)->degree - ((GAMGNode*)b)->degree);
25: }
27: /* -------------------------------------------------------------------------- */
28: /*
29: PCSetCoordinates_GEO
31: Input Parameter:
32: . pc - the preconditioner context
33: */
34: PetscErrorCode PCSetCoordinates_GEO(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
35: {
36: PC_MG *mg = (PC_MG*)pc->data;
37: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
39: PetscInt arrsz,bs,my0,kk,ii,nloc,Iend,aloc;
40: Mat Amat = pc->pmat;
44: MatGetBlockSize(Amat, &bs);
45: MatGetOwnershipRange(Amat, &my0, &Iend);
46: aloc = (Iend-my0);
47: nloc = (Iend-my0)/bs;
49: if (nloc!=a_nloc && aloc!=a_nloc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Number of local blocks %D must be %D or %D.",a_nloc,nloc,aloc);
51: pc_gamg->data_cell_rows = 1;
52: if (!coords && nloc > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Need coordinates for pc_gamg_type 'geo'.");
53: pc_gamg->data_cell_cols = ndm; /* coordinates */
55: arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;
57: /* create data - syntactic sugar that should be refactored at some point */
58: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
59: PetscFree(pc_gamg->data);
60: PetscMalloc1(arrsz+1, &pc_gamg->data);
61: }
62: for (kk=0; kk<arrsz; kk++) pc_gamg->data[kk] = -999.;
63: pc_gamg->data[arrsz] = -99.;
64: /* copy data in - column oriented */
65: if (nloc == a_nloc) {
66: for (kk = 0; kk < nloc; kk++) {
67: for (ii = 0; ii < ndm; ii++) {
68: pc_gamg->data[ii*nloc + kk] = coords[kk*ndm + ii];
69: }
70: }
71: } else { /* assumes the coordinates are blocked */
72: for (kk = 0; kk < nloc; kk++) {
73: for (ii = 0; ii < ndm; ii++) {
74: pc_gamg->data[ii*nloc + kk] = coords[bs*kk*ndm + ii];
75: }
76: }
77: }
78: if (pc_gamg->data[arrsz] != -99.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->data[arrsz %D] %g != -99.",arrsz,pc_gamg->data[arrsz]);
79: pc_gamg->data_sz = arrsz;
80: return(0);
81: }
83: /* -------------------------------------------------------------------------- */
84: /*
85: PCSetData_GEO
87: Input Parameter:
88: . pc -
89: */
90: PetscErrorCode PCSetData_GEO(PC pc, Mat m)
91: {
93: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"GEO MG needs coordinates");
94: }
96: /* -------------------------------------------------------------------------- */
97: /*
98: PCSetFromOptions_GEO
100: Input Parameter:
101: . pc -
102: */
103: PetscErrorCode PCSetFromOptions_GEO(PetscOptionItems *PetscOptionsObject,PC pc)
104: {
108: PetscOptionsHead(PetscOptionsObject,"GAMG-GEO options");
109: {
110: /* -pc_gamg_sa_nsmooths */
111: /* pc_gamg_sa->smooths = 0; */
112: /* PetscOptionsInt("-pc_gamg_agg_nsmooths", */
113: /* "smoothing steps for smoothed aggregation, usually 1 (0)", */
114: /* "PCGAMGSetNSmooths_AGG", */
115: /* pc_gamg_sa->smooths, */
116: /* &pc_gamg_sa->smooths, */
117: /* &flag); */
118: /* */
119: }
120: PetscOptionsTail();
121: return(0);
122: }
124: /* -------------------------------------------------------------------------- */
125: /*
126: triangulateAndFormProl
128: Input Parameter:
129: . selected_2 - list of selected local ID, includes selected ghosts
130: . data_stride -
131: . coords[2*data_stride] - column vector of local coordinates w/ ghosts
132: . nselected_1 - selected IDs that go with base (1) graph includes selected ghosts
133: . clid_lid_1[nselected_1] - lids of selected (c) nodes ???????????
134: . agg_lists_1 - list of aggregates selected_1 vertices of aggregate unselected vertices
135: . crsGID[selected.size()] - global index for prolongation operator
136: . bs - block size
137: Output Parameter:
138: . a_Prol - prolongation operator
139: . a_worst_best - measure of worst missed fine vertex, 0 is no misses
140: */
141: static PetscErrorCode triangulateAndFormProl(IS selected_2,PetscInt data_stride,PetscReal coords[],PetscInt nselected_1,const PetscInt clid_lid_1[],const PetscCoarsenData *agg_lists_1,
142: const PetscInt crsGID[],PetscInt bs,Mat a_Prol,PetscReal *a_worst_best)
143: {
144: #if defined(PETSC_HAVE_TRIANGLE)
145: PetscErrorCode ierr;
146: PetscInt jj,tid,tt,idx,nselected_2;
147: struct triangulateio in,mid;
148: const PetscInt *selected_idx_2;
149: PetscMPIInt rank;
150: PetscInt Istart,Iend,nFineLoc,myFine0;
151: int kk,nPlotPts,sid;
152: MPI_Comm comm;
153: PetscReal tm;
156: PetscObjectGetComm((PetscObject)a_Prol,&comm);
157: MPI_Comm_rank(comm,&rank);
158: ISGetSize(selected_2, &nselected_2);
159: if (nselected_2 == 1 || nselected_2 == 2) { /* 0 happens on idle processors */
160: *a_worst_best = 100.0; /* this will cause a stop, but not globalized (should not happen) */
161: } else *a_worst_best = 0.0;
162: MPIU_Allreduce(a_worst_best, &tm, 1, MPIU_REAL, MPIU_MAX, comm);
163: if (tm > 0.0) {
164: *a_worst_best = 100.0;
165: return(0);
166: }
167: MatGetOwnershipRange(a_Prol, &Istart, &Iend);
168: nFineLoc = (Iend-Istart)/bs; myFine0 = Istart/bs;
169: nPlotPts = nFineLoc; /* locals */
170: /* triangle */
171: /* Define input points - in*/
172: in.numberofpoints = nselected_2;
173: in.numberofpointattributes = 0;
174: /* get nselected points */
175: PetscMalloc1(2*nselected_2, &in.pointlist);
176: ISGetIndices(selected_2, &selected_idx_2);
178: for (kk=0,sid=0; kk<nselected_2; kk++,sid += 2) {
179: PetscInt lid = selected_idx_2[kk];
180: in.pointlist[sid] = coords[lid];
181: in.pointlist[sid+1] = coords[data_stride + lid];
182: if (lid>=nFineLoc) nPlotPts++;
183: }
184: if (sid != 2*nselected_2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %D != 2*nselected_2 %D",sid,nselected_2);
186: in.numberofsegments = 0;
187: in.numberofedges = 0;
188: in.numberofholes = 0;
189: in.numberofregions = 0;
190: in.trianglelist = 0;
191: in.segmentmarkerlist = 0;
192: in.pointattributelist = 0;
193: in.pointmarkerlist = 0;
194: in.triangleattributelist = 0;
195: in.trianglearealist = 0;
196: in.segmentlist = 0;
197: in.holelist = 0;
198: in.regionlist = 0;
199: in.edgelist = 0;
200: in.edgemarkerlist = 0;
201: in.normlist = 0;
203: /* triangulate */
204: mid.pointlist = 0; /* Not needed if -N switch used. */
205: /* Not needed if -N switch used or number of point attributes is zero: */
206: mid.pointattributelist = 0;
207: mid.pointmarkerlist = 0; /* Not needed if -N or -B switch used. */
208: mid.trianglelist = 0; /* Not needed if -E switch used. */
209: /* Not needed if -E switch used or number of triangle attributes is zero: */
210: mid.triangleattributelist = 0;
211: mid.neighborlist = 0; /* Needed only if -n switch used. */
212: /* Needed only if segments are output (-p or -c) and -P not used: */
213: mid.segmentlist = 0;
214: /* Needed only if segments are output (-p or -c) and -P and -B not used: */
215: mid.segmentmarkerlist = 0;
216: mid.edgelist = 0; /* Needed only if -e switch used. */
217: mid.edgemarkerlist = 0; /* Needed if -e used and -B not used. */
218: mid.numberoftriangles = 0;
220: /* Triangulate the points. Switches are chosen to read and write a */
221: /* PSLG (p), preserve the convex hull (c), number everything from */
222: /* zero (z), assign a regional attribute to each element (A), and */
223: /* produce an edge list (e), a Voronoi diagram (v), and a triangle */
224: /* neighbor list (n). */
225: if (nselected_2 != 0) { /* inactive processor */
226: char args[] = "npczQ"; /* c is needed ? */
227: triangulate(args, &in, &mid, (struct triangulateio*) NULL);
228: /* output .poly files for 'showme' */
229: if (!PETSC_TRUE) {
230: static int level = 1;
231: FILE *file; char fname[32];
233: sprintf(fname,"C%d_%d.poly",level,rank); file = fopen(fname, "w");
234: /*First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>*/
235: fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0);
236: /*Following lines: <vertex #> <x> <y> */
237: for (kk=0,sid=0; kk<in.numberofpoints; kk++,sid += 2) {
238: fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]);
239: }
240: /*One line: <# of segments> <# of boundary markers (0 or 1)> */
241: fprintf(file, "%d %d\n",0,0);
242: /*Following lines: <segment #> <endpoint> <endpoint> [boundary marker] */
243: /* One line: <# of holes> */
244: fprintf(file, "%d\n",0);
245: /* Following lines: <hole #> <x> <y> */
246: /* Optional line: <# of regional attributes and/or area constraints> */
247: /* Optional following lines: <region #> <x> <y> <attribute> <maximum area> */
248: fclose(file);
250: /* elems */
251: sprintf(fname,"C%d_%d.ele",level,rank); file = fopen(fname, "w");
252: /* First line: <# of triangles> <nodes per triangle> <# of attributes> */
253: fprintf(file, "%d %d %d\n",mid.numberoftriangles,3,0);
254: /* Remaining lines: <triangle #> <node> <node> <node> ... [attributes] */
255: for (kk=0,sid=0; kk<mid.numberoftriangles; kk++,sid += 3) {
256: fprintf(file, "%d %d %d %d\n",kk,mid.trianglelist[sid],mid.trianglelist[sid+1],mid.trianglelist[sid+2]);
257: }
258: fclose(file);
260: sprintf(fname,"C%d_%d.node",level,rank); file = fopen(fname, "w");
261: /* First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> */
262: /* fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0); */
263: fprintf(file, "%d %d %d %d\n",nPlotPts,2,0,0);
264: /*Following lines: <vertex #> <x> <y> */
265: for (kk=0,sid=0; kk<in.numberofpoints; kk++,sid+=2) {
266: fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]);
267: }
269: sid /= 2;
270: for (jj=0; jj<nFineLoc; jj++) {
271: PetscBool sel = PETSC_TRUE;
272: for (kk=0; kk<nselected_2 && sel; kk++) {
273: PetscInt lid = selected_idx_2[kk];
274: if (lid == jj) sel = PETSC_FALSE;
275: }
276: if (sel) fprintf(file, "%d %e %e\n",sid++,coords[jj],coords[data_stride + jj]);
277: }
278: fclose(file);
279: if (sid != nPlotPts) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %D != nPlotPts %D",sid,nPlotPts);
280: level++;
281: }
282: }
283: #if defined PETSC_GAMG_USE_LOG
284: PetscLogEventBegin(petsc_gamg_setup_events[FIND_V],0,0,0,0);
285: #endif
286: { /* form P - setup some maps */
287: PetscInt clid,mm,*nTri,*node_tri;
289: PetscMalloc2(nselected_2, &node_tri,nselected_2, &nTri);
291: /* need list of triangles on node */
292: for (kk=0; kk<nselected_2; kk++) nTri[kk] = 0;
293: for (tid=0,kk=0; tid<mid.numberoftriangles; tid++) {
294: for (jj=0; jj<3; jj++) {
295: PetscInt cid = mid.trianglelist[kk++];
296: if (nTri[cid] == 0) node_tri[cid] = tid;
297: nTri[cid]++;
298: }
299: }
300: #define EPS 1.e-12
301: /* find points and set prolongation */
302: for (mm = clid = 0; mm < nFineLoc; mm++) {
303: PetscBool ise;
304: PetscCDEmptyAt(agg_lists_1,mm,&ise);
305: if (!ise) {
306: const PetscInt lid = mm;
307: PetscScalar AA[3][3];
308: PetscBLASInt N=3,NRHS=1,LDA=3,IPIV[3],LDB=3,INFO;
309: PetscCDIntNd *pos;
311: PetscCDGetHeadPos(agg_lists_1,lid,&pos);
312: while (pos) {
313: PetscInt flid;
314: PetscCDIntNdGetID(pos, &flid);
315: PetscCDGetNextPos(agg_lists_1,lid,&pos);
317: if (flid < nFineLoc) { /* could be a ghost */
318: PetscInt bestTID = -1; PetscReal best_alpha = 1.e10;
319: const PetscInt fgid = flid + myFine0;
320: /* compute shape function for gid */
321: const PetscReal fcoord[3] = {coords[flid],coords[data_stride+flid],1.0};
322: PetscBool haveit =PETSC_FALSE; PetscScalar alpha[3]; PetscInt clids[3];
324: /* look for it */
325: for (tid = node_tri[clid], jj=0;
326: jj < 5 && !haveit && tid != -1;
327: jj++) {
328: for (tt=0; tt<3; tt++) {
329: PetscInt cid2 = mid.trianglelist[3*tid + tt];
330: PetscInt lid2 = selected_idx_2[cid2];
331: AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
332: clids[tt] = cid2; /* store for interp */
333: }
335: for (tt=0; tt<3; tt++) alpha[tt] = (PetscScalar)fcoord[tt];
337: /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
338: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
339: {
340: PetscBool have=PETSC_TRUE; PetscReal lowest=1.e10;
341: for (tt = 0, idx = 0; tt < 3; tt++) {
342: if (PetscRealPart(alpha[tt]) > (1.0+EPS) || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE;
343: if (PetscRealPart(alpha[tt]) < lowest) {
344: lowest = PetscRealPart(alpha[tt]);
345: idx = tt;
346: }
347: }
348: haveit = have;
349: }
350: tid = mid.neighborlist[3*tid + idx];
351: }
353: if (!haveit) {
354: /* brute force */
355: for (tid=0; tid<mid.numberoftriangles && !haveit; tid++) {
356: for (tt=0; tt<3; tt++) {
357: PetscInt cid2 = mid.trianglelist[3*tid + tt];
358: PetscInt lid2 = selected_idx_2[cid2];
359: AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
360: clids[tt] = cid2; /* store for interp */
361: }
362: for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt];
363: /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
364: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
365: {
366: PetscBool have=PETSC_TRUE; PetscReal worst=0.0, v;
367: for (tt=0; tt<3 && have; tt++) {
368: if (PetscRealPart(alpha[tt]) > 1.0+EPS || PetscRealPart(alpha[tt]) < -EPS) have=PETSC_FALSE;
369: if ((v=PetscAbs(PetscRealPart(alpha[tt])-0.5)) > worst) worst = v;
370: }
371: if (worst < best_alpha) {
372: best_alpha = worst; bestTID = tid;
373: }
374: haveit = have;
375: }
376: }
377: }
378: if (!haveit) {
379: if (best_alpha > *a_worst_best) *a_worst_best = best_alpha;
380: /* use best one */
381: for (tt=0; tt<3; tt++) {
382: PetscInt cid2 = mid.trianglelist[3*bestTID + tt];
383: PetscInt lid2 = selected_idx_2[cid2];
384: AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
385: clids[tt] = cid2; /* store for interp */
386: }
387: for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt];
388: /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
389: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
390: }
392: /* put in row of P */
393: for (idx=0; idx<3; idx++) {
394: PetscScalar shp = alpha[idx];
395: if (PetscAbs(PetscRealPart(shp)) > 1.e-6) {
396: PetscInt cgid = crsGID[clids[idx]];
397: PetscInt jj = cgid*bs, ii = fgid*bs; /* need to gloalize */
398: for (tt=0; tt < bs; tt++, ii++, jj++) {
399: MatSetValues(a_Prol,1,&ii,1,&jj,&shp,INSERT_VALUES);
400: }
401: }
402: }
403: }
404: } /* aggregates iterations */
405: clid++;
406: } /* a coarse agg */
407: } /* for all fine nodes */
409: ISRestoreIndices(selected_2, &selected_idx_2);
410: MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
411: MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);
413: PetscFree2(node_tri,nTri);
414: }
415: #if defined PETSC_GAMG_USE_LOG
416: PetscLogEventEnd(petsc_gamg_setup_events[FIND_V],0,0,0,0);
417: #endif
418: free(mid.trianglelist);
419: free(mid.neighborlist);
420: free(mid.segmentlist);
421: free(mid.segmentmarkerlist);
422: free(mid.pointlist);
423: free(mid.pointmarkerlist);
424: PetscFree(in.pointlist);
425: return(0);
426: #else
427: SETERRQ(PetscObjectComm((PetscObject)a_Prol),PETSC_ERR_PLIB,"configure with TRIANGLE to use geometric MG");
428: #endif
429: }
430: /* -------------------------------------------------------------------------- */
431: /*
432: getGIDsOnSquareGraph - square graph, get
434: Input Parameter:
435: . nselected_1 - selected local indices (includes ghosts in input Gmat1)
436: . clid_lid_1 - [nselected_1] lids of selected nodes
437: . Gmat1 - graph that goes with 'selected_1'
438: Output Parameter:
439: . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2)
440: . a_Gmat_2 - graph that is squared of 'Gmat_1'
441: . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes
442: */
443: static PetscErrorCode getGIDsOnSquareGraph(PetscInt nselected_1,const PetscInt clid_lid_1[],const Mat Gmat1,IS *a_selected_2,Mat *a_Gmat_2,PetscInt **a_crsGID)
444: {
446: PetscMPIInt size;
447: PetscInt *crsGID, kk,my0,Iend,nloc;
448: MPI_Comm comm;
451: PetscObjectGetComm((PetscObject)Gmat1,&comm);
452: MPI_Comm_size(comm,&size);
453: MatGetOwnershipRange(Gmat1,&my0,&Iend); /* AIJ */
454: nloc = Iend - my0; /* this does not change */
456: if (size == 1) { /* not much to do in serial */
457: PetscMalloc1(nselected_1, &crsGID);
458: for (kk=0; kk<nselected_1; kk++) crsGID[kk] = kk;
459: *a_Gmat_2 = 0;
460: ISCreateGeneral(PETSC_COMM_SELF,nselected_1,clid_lid_1,PETSC_COPY_VALUES,a_selected_2);
461: } else {
462: PetscInt idx,num_fine_ghosts,num_crs_ghost,myCrs0;
463: Mat_MPIAIJ *mpimat2;
464: Mat Gmat2;
465: Vec locState;
466: PetscScalar *cpcol_state;
468: /* scan my coarse zero gid, set 'lid_state' with coarse GID */
469: kk = nselected_1;
470: MPI_Scan(&kk, &myCrs0, 1, MPIU_INT, MPI_SUM, comm);
471: myCrs0 -= nselected_1;
473: if (a_Gmat_2) { /* output */
474: /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */
475: MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);
476: *a_Gmat_2 = Gmat2; /* output */
477: } else Gmat2 = Gmat1; /* use local to get crsGIDs at least */
478: /* get coarse grid GIDS for selected (locals and ghosts) */
479: mpimat2 = (Mat_MPIAIJ*)Gmat2->data;
480: MatCreateVecs(Gmat2, &locState, 0);
481: VecSet(locState, (PetscScalar)(PetscReal)(-1)); /* set with UNKNOWN state */
482: for (kk=0; kk<nselected_1; kk++) {
483: PetscInt fgid = clid_lid_1[kk] + my0;
484: PetscScalar v = (PetscScalar)(kk+myCrs0);
485: VecSetValues(locState, 1, &fgid, &v, INSERT_VALUES); /* set with PID */
486: }
487: VecAssemblyBegin(locState);
488: VecAssemblyEnd(locState);
489: VecScatterBegin(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);
490: VecScatterEnd(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);
491: VecGetLocalSize(mpimat2->lvec, &num_fine_ghosts);
492: VecGetArray(mpimat2->lvec, &cpcol_state);
493: for (kk=0,num_crs_ghost=0; kk<num_fine_ghosts; kk++) {
494: if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) num_crs_ghost++;
495: }
496: PetscMalloc1(nselected_1+num_crs_ghost, &crsGID); /* output */
497: {
498: PetscInt *selected_set;
499: PetscMalloc1(nselected_1+num_crs_ghost, &selected_set);
500: /* do ghost of 'crsGID' */
501: for (kk=0,idx=nselected_1; kk<num_fine_ghosts; kk++) {
502: if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
503: PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
504: selected_set[idx] = nloc + kk;
505: crsGID[idx++] = cgid;
506: }
507: }
508: if (idx != (nselected_1+num_crs_ghost)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"idx %D != (nselected_1 %D + num_crs_ghost %D)",idx,nselected_1,num_crs_ghost);
509: VecRestoreArray(mpimat2->lvec, &cpcol_state);
510: /* do locals in 'crsGID' */
511: VecGetArray(locState, &cpcol_state);
512: for (kk=0,idx=0; kk<nloc; kk++) {
513: if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
514: PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
515: selected_set[idx] = kk;
516: crsGID[idx++] = cgid;
517: }
518: }
519: if (idx != nselected_1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"idx %D != nselected_1 %D",idx,nselected_1);
520: VecRestoreArray(locState, &cpcol_state);
522: if (a_selected_2 != 0) { /* output */
523: ISCreateGeneral(PETSC_COMM_SELF,(nselected_1+num_crs_ghost),selected_set,PETSC_OWN_POINTER,a_selected_2);
524: } else {
525: PetscFree(selected_set);
526: }
527: }
528: VecDestroy(&locState);
529: }
530: *a_crsGID = crsGID; /* output */
531: return(0);
532: }
534: /* -------------------------------------------------------------------------- */
535: /*
536: PCGAMGGraph_GEO
538: Input Parameter:
539: . pc - this
540: . Amat - matrix on this fine level
541: Output Parameter:
542: . a_Gmat
543: */
544: PetscErrorCode PCGAMGGraph_GEO(PC pc,Mat Amat,Mat *a_Gmat)
545: {
546: PetscErrorCode ierr;
547: PC_MG *mg = (PC_MG*)pc->data;
548: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
549: const PetscReal vfilter = pc_gamg->threshold[0];
550: MPI_Comm comm;
551: Mat Gmat;
552: PetscBool set,flg,symm;
555: PetscObjectGetComm((PetscObject)Amat,&comm);
556: PetscLogEventBegin(PC_GAMGGraph_GEO,0,0,0,0);
558: MatIsSymmetricKnown(Amat, &set, &flg);
559: symm = (PetscBool)!(set && flg);
561: PCGAMGCreateGraph(Amat, &Gmat);
562: PCGAMGFilterGraph(&Gmat, vfilter, symm);
564: *a_Gmat = Gmat;
565: PetscLogEventEnd(PC_GAMGGraph_GEO,0,0,0,0);
566: return(0);
567: }
569: /* -------------------------------------------------------------------------- */
570: /*
571: PCGAMGCoarsen_GEO
573: Input Parameter:
574: . a_pc - this
575: . a_Gmat - graph
576: Output Parameter:
577: . a_llist_parent - linked list from selected indices for data locality only
578: */
579: PetscErrorCode PCGAMGCoarsen_GEO(PC a_pc,Mat *a_Gmat,PetscCoarsenData **a_llist_parent)
580: {
582: PetscInt Istart,Iend,nloc,kk,Ii,ncols;
583: IS perm;
584: GAMGNode *gnodes;
585: PetscInt *permute;
586: Mat Gmat = *a_Gmat;
587: MPI_Comm comm;
588: MatCoarsen crs;
591: PetscObjectGetComm((PetscObject)a_pc,&comm);
592: PetscLogEventBegin(PC_GAMGCoarsen_GEO,0,0,0,0);
593: MatGetOwnershipRange(Gmat, &Istart, &Iend);
594: nloc = (Iend-Istart);
596: /* create random permutation with sort for geo-mg */
597: PetscMalloc1(nloc, &gnodes);
598: PetscMalloc1(nloc, &permute);
600: for (Ii=Istart; Ii<Iend; Ii++) { /* locals only? */
601: MatGetRow(Gmat,Ii,&ncols,0,0);
602: {
603: PetscInt lid = Ii - Istart;
604: gnodes[lid].lid = lid;
605: gnodes[lid].degree = ncols;
606: }
607: MatRestoreRow(Gmat,Ii,&ncols,0,0);
608: }
609: if (PETSC_TRUE) {
610: PetscRandom rand;
611: PetscBool *bIndexSet;
612: PetscReal rr;
613: PetscInt iSwapIndex;
615: PetscRandomCreate(comm,&rand);
616: PetscCalloc1(nloc, &bIndexSet);
617: for (Ii = 0; Ii < nloc; Ii++) {
618: PetscRandomGetValueReal(rand,&rr);
619: iSwapIndex = (PetscInt) (rr*nloc);
620: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
621: GAMGNode iTemp = gnodes[iSwapIndex];
622: gnodes[iSwapIndex] = gnodes[Ii];
623: gnodes[Ii] = iTemp;
624: bIndexSet[Ii] = PETSC_TRUE;
625: bIndexSet[iSwapIndex] = PETSC_TRUE;
626: }
627: }
628: PetscRandomDestroy(&rand);
629: PetscFree(bIndexSet);
630: }
631: /* only sort locals */
632: qsort(gnodes, nloc, sizeof(GAMGNode), petsc_geo_mg_compare);
633: /* create IS of permutation */
634: for (kk=0; kk<nloc; kk++) permute[kk] = gnodes[kk].lid; /* locals only */
635: PetscFree(gnodes);
636: ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_OWN_POINTER, &perm);
638: /* get MIS aggs */
640: MatCoarsenCreate(comm, &crs);
641: MatCoarsenSetType(crs, MATCOARSENMIS);
642: MatCoarsenSetGreedyOrdering(crs, perm);
643: MatCoarsenSetAdjacency(crs, Gmat);
644: MatCoarsenSetStrictAggs(crs, PETSC_FALSE);
645: MatCoarsenApply(crs);
646: MatCoarsenGetData(crs, a_llist_parent);
647: MatCoarsenDestroy(&crs);
649: ISDestroy(&perm);
650: PetscLogEventEnd(PC_GAMGCoarsen_GEO,0,0,0,0);
651: return(0);
652: }
654: /* -------------------------------------------------------------------------- */
655: /*
656: PCGAMGProlongator_GEO
658: Input Parameter:
659: . pc - this
660: . Amat - matrix on this fine level
661: . Graph - used to get ghost data for nodes in
662: . selected_1 - [nselected]
663: . agg_lists - [nselected]
664: Output Parameter:
665: . a_P_out - prolongation operator to the next level
666: */
667: PetscErrorCode PCGAMGProlongator_GEO(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
668: {
669: PC_MG *mg = (PC_MG*)pc->data;
670: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
671: const PetscInt dim = pc_gamg->data_cell_cols, data_cols = pc_gamg->data_cell_cols;
673: PetscInt Istart,Iend,nloc,my0,jj,kk,ncols,nLocalSelected,bs,*clid_flid;
674: Mat Prol;
675: PetscMPIInt rank, size;
676: MPI_Comm comm;
677: IS selected_2,selected_1;
678: const PetscInt *selected_idx;
679: MatType mtype;
682: PetscObjectGetComm((PetscObject)Amat,&comm);
683: PetscLogEventBegin(PC_GAMGProlongator_GEO,0,0,0,0);
684: MPI_Comm_rank(comm,&rank);
685: MPI_Comm_size(comm,&size);
686: MatGetOwnershipRange(Amat, &Istart, &Iend);
687: MatGetBlockSize(Amat, &bs);
688: nloc = (Iend-Istart)/bs; my0 = Istart/bs;
689: if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) % bs %D",Iend,Istart,bs);
691: /* get 'nLocalSelected' */
692: PetscCDGetMIS(agg_lists, &selected_1);
693: ISGetSize(selected_1, &jj);
694: PetscMalloc1(jj, &clid_flid);
695: ISGetIndices(selected_1, &selected_idx);
696: for (kk=0,nLocalSelected=0; kk<jj; kk++) {
697: PetscInt lid = selected_idx[kk];
698: if (lid<nloc) {
699: MatGetRow(Gmat,lid+my0,&ncols,0,0);
700: if (ncols>1) clid_flid[nLocalSelected++] = lid; /* fiter out singletons */
701: MatRestoreRow(Gmat,lid+my0,&ncols,0,0);
702: }
703: }
704: ISRestoreIndices(selected_1, &selected_idx);
705: ISDestroy(&selected_1); /* this is selected_1 in serial */
707: /* create prolongator matrix */
708: MatGetType(Amat,&mtype);
709: MatCreate(comm, &Prol);
710: MatSetSizes(Prol,nloc*bs,nLocalSelected*bs,PETSC_DETERMINE,PETSC_DETERMINE);
711: MatSetBlockSizes(Prol, bs, bs);
712: MatSetType(Prol, mtype);
713: MatSeqAIJSetPreallocation(Prol,3*data_cols,NULL);
714: MatMPIAIJSetPreallocation(Prol,3*data_cols,NULL,3*data_cols,NULL);
716: /* can get all points "removed" - but not on geomg */
717: MatGetSize(Prol, &kk, &jj);
718: if (!jj) {
719: PetscInfo(pc,"ERROE: no selected points on coarse grid\n");
720: PetscFree(clid_flid);
721: MatDestroy(&Prol);
722: *a_P_out = NULL; /* out */
723: return(0);
724: }
726: {
727: PetscReal *coords;
728: PetscInt data_stride;
729: PetscInt *crsGID = NULL;
730: Mat Gmat2;
732: if (dim != data_cols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"dim %D != data_cols %D",dim,data_cols);
733: /* grow ghost data for better coarse grid cover of fine grid */
734: #if defined PETSC_GAMG_USE_LOG
735: PetscLogEventBegin(petsc_gamg_setup_events[SET5],0,0,0,0);
736: #endif
737: /* messy method, squares graph and gets some data */
738: getGIDsOnSquareGraph(nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID);
739: /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */
740: #if defined PETSC_GAMG_USE_LOG
741: PetscLogEventEnd(petsc_gamg_setup_events[SET5],0,0,0,0);
742: #endif
743: /* create global vector of coorindates in 'coords' */
744: if (size > 1) {
745: PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords);
746: } else {
747: coords = (PetscReal*)pc_gamg->data;
748: data_stride = pc_gamg->data_sz/pc_gamg->data_cell_cols;
749: }
750: MatDestroy(&Gmat2);
752: /* triangulate */
753: if (dim == 2) {
754: PetscReal metric,tm;
755: #if defined PETSC_GAMG_USE_LOG
756: PetscLogEventBegin(petsc_gamg_setup_events[SET6],0,0,0,0);
757: #endif
758: triangulateAndFormProl(selected_2, data_stride, coords,nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric);
759: #if defined PETSC_GAMG_USE_LOG
760: PetscLogEventEnd(petsc_gamg_setup_events[SET6],0,0,0,0);
761: #endif
762: PetscFree(crsGID);
764: /* clean up and create coordinates for coarse grid (output) */
765: if (size > 1) PetscFree(coords);
767: MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm);
768: if (tm > 1.) { /* needs to be globalized - should not happen */
769: PetscInfo1(pc," failed metric for coarse grid %e\n",(double)tm);
770: MatDestroy(&Prol);
771: } else if (metric > .0) {
772: PetscInfo1(pc,"worst metric for coarse grid = %e\n",(double)metric);
773: }
774: } else SETERRQ(comm,PETSC_ERR_PLIB,"3D not implemented for 'geo' AMG");
775: { /* create next coords - output */
776: PetscReal *crs_crds;
777: PetscMalloc1(dim*nLocalSelected, &crs_crds);
778: for (kk=0; kk<nLocalSelected; kk++) { /* grab local select nodes to promote - output */
779: PetscInt lid = clid_flid[kk];
780: for (jj=0; jj<dim; jj++) crs_crds[jj*nLocalSelected + kk] = pc_gamg->data[jj*nloc + lid];
781: }
783: PetscFree(pc_gamg->data);
784: pc_gamg->data = crs_crds; /* out */
785: pc_gamg->data_sz = dim*nLocalSelected;
786: }
787: ISDestroy(&selected_2);
788: }
790: *a_P_out = Prol; /* out */
791: PetscFree(clid_flid);
792: PetscLogEventEnd(PC_GAMGProlongator_GEO,0,0,0,0);
793: return(0);
794: }
796: static PetscErrorCode PCDestroy_GAMG_GEO(PC pc)
797: {
801: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
802: return(0);
803: }
805: /* -------------------------------------------------------------------------- */
806: /*
807: PCCreateGAMG_GEO
809: Input Parameter:
810: . pc -
811: */
812: PetscErrorCode PCCreateGAMG_GEO(PC pc)
813: {
815: PC_MG *mg = (PC_MG*)pc->data;
816: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
819: pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO;
820: pc_gamg->ops->destroy = PCDestroy_GAMG_GEO;
821: /* reset does not do anything; setup not virtual */
823: /* set internal function pointers */
824: pc_gamg->ops->graph = PCGAMGGraph_GEO;
825: pc_gamg->ops->coarsen = PCGAMGCoarsen_GEO;
826: pc_gamg->ops->prolongator = PCGAMGProlongator_GEO;
827: pc_gamg->ops->optprolongator = NULL;
828: pc_gamg->ops->createdefaultdata = PCSetData_GEO;
830: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_GEO);
831: return(0);
832: }