Actual source code: geo.c
1: /*
2: GAMG geometric-algebraic multigrid PC - Mark Adams 2011
3: */
5: #include <../src/ksp/pc/impls/gamg/gamg.h>
6: #include <petscsf.h>
8: #if PetscDefined(HAVE_TRIANGLE)
9: #if !defined(ANSI_DECLARATORS)
10: #define ANSI_DECLARATORS
11: #endif
12: #include <triangle.h>
13: #endif
15: #include <petscblaslapack.h>
17: /* Private context for the GAMG preconditioner */
18: typedef struct {
19: PetscInt lid; /* local vertex index */
20: PetscInt degree; /* vertex degree */
21: } GAMGNode;
23: static inline int petsc_geo_mg_compare(const void *a, const void *b)
24: {
25: return (int)(((GAMGNode *)a)->degree - ((GAMGNode *)b)->degree);
26: }
28: // PetscClangLinter pragma disable: -fdoc-sowing-chars
29: /*
30: PCSetCoordinates_GEO
32: Input Parameter:
33: . pc - the preconditioner context
34: */
35: static PetscErrorCode PCSetCoordinates_GEO(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
36: {
37: PC_MG *mg = (PC_MG *)pc->data;
38: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
39: PetscInt arrsz, bs, my0, kk, ii, nloc, Iend, aloc;
40: Mat Amat = pc->pmat;
42: PetscFunctionBegin;
44: PetscCall(MatGetBlockSize(Amat, &bs));
45: PetscCall(MatGetOwnershipRange(Amat, &my0, &Iend));
46: aloc = (Iend - my0);
47: nloc = (Iend - my0) / bs;
49: PetscCheck(nloc == a_nloc || aloc == a_nloc, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of local blocks %" PetscInt_FMT " must be %" PetscInt_FMT " or %" PetscInt_FMT ".", a_nloc, nloc, aloc);
51: pc_gamg->data_cell_rows = 1;
52: PetscCheck(coords || (nloc <= 0), 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: PetscCall(PetscFree(pc_gamg->data));
60: PetscCall(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++) pc_gamg->data[ii * nloc + kk] = coords[kk * ndm + ii];
68: }
69: } else { /* assumes the coordinates are blocked */
70: for (kk = 0; kk < nloc; kk++) {
71: for (ii = 0; ii < ndm; ii++) pc_gamg->data[ii * nloc + kk] = coords[bs * kk * ndm + ii];
72: }
73: }
74: PetscCheck(pc_gamg->data[arrsz] == -99., PETSC_COMM_SELF, PETSC_ERR_PLIB, "pc_gamg->data[arrsz %" PetscInt_FMT "] %g != -99.", arrsz, (double)pc_gamg->data[arrsz]);
75: pc_gamg->data_sz = arrsz;
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: // PetscClangLinter pragma disable: -fdoc-sowing-chars
80: /*
81: PCSetData_GEO
83: Input Parameter:
84: . pc -
85: */
86: static PetscErrorCode PCSetData_GEO(PC pc, Mat m)
87: {
88: PetscFunctionBegin;
89: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "GEO MG needs coordinates");
90: }
92: static PetscErrorCode PCSetFromOptions_GEO(PC pc, PetscOptionItems PetscOptionsObject)
93: {
94: PetscFunctionBegin;
95: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-GEO options");
96: {
97: /* -pc_gamg_sa_nsmooths */
98: /* pc_gamg_sa->smooths = 0; */
99: /* ierr = PetscOptionsInt("-pc_gamg_agg_nsmooths", */
100: /* "smoothing steps for smoothed aggregation, usually 1 (0)", */
101: /* "PCGAMGSetNSmooths_AGG", */
102: /* pc_gamg_sa->smooths, */
103: /* &pc_gamg_sa->smooths, */
104: /* &flag); */
105: }
106: PetscOptionsHeadEnd();
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: // PetscClangLinter pragma disable: -fdoc-sowing-chars
111: /*
112: triangulateAndFormProl
114: Input Parameter:
115: . selected_2 - list of selected local ID, includes selected ghosts
116: . data_stride -
117: . coords[2*data_stride] - column vector of local coordinates w/ ghosts
118: . nselected_1 - selected IDs that go with base (1) graph includes selected ghosts
119: . clid_lid_1[nselected_1] - lids of selected (c) nodes ???????????
120: . agg_lists_1 - list of aggregates selected_1 vertices of aggregate unselected vertices
121: . crsGID[selected.size()] - global index for prolongation operator
122: . bs - block size
123: Output Parameter:
124: . a_Prol - prolongation operator
125: . a_worst_best - measure of worst missed fine vertex, 0 is no misses
126: */
127: static PetscErrorCode triangulateAndFormProl(IS selected_2, PetscInt data_stride, PetscReal coords[], PetscInt nselected_1, const PetscInt clid_lid_1[], const PetscCoarsenData *agg_lists_1, const PetscInt crsGID[], PetscInt bs, Mat a_Prol, PetscReal *a_worst_best)
128: {
129: #if PetscDefined(HAVE_TRIANGLE)
130: PetscInt jj, tid, tt, idx, nselected_2;
131: struct triangulateio in, mid;
132: const PetscInt *selected_idx_2;
133: PetscMPIInt rank;
134: PetscInt Istart, Iend, nFineLoc, myFine0;
135: int kk, nPlotPts, sid;
136: MPI_Comm comm;
138: PetscFunctionBegin;
139: PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
140: PetscCallMPI(MPI_Comm_rank(comm, &rank));
141: PetscCall(ISGetSize(selected_2, &nselected_2));
142: if (nselected_2 == 1 || nselected_2 == 2) { /* 0 happens on idle processors */
143: *a_worst_best = 100.0; /* this will cause a stop, but not globalized (should not happen) */
144: } else *a_worst_best = 0.0;
145: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, a_worst_best, 1, MPIU_REAL, MPIU_MAX, comm));
146: if (*a_worst_best > 0.0) {
147: *a_worst_best = 100.0;
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
150: PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
151: nFineLoc = (Iend - Istart) / bs;
152: myFine0 = Istart / bs;
153: PetscCall(PetscCIntCast(nFineLoc, &nPlotPts)); /* locals */
154: /* triangle */
155: /* Define input points - in */
156: PetscCall(PetscCIntCast(nselected_2, &in.numberofpoints));
157: in.numberofpointattributes = 0;
158: /* get nselected points */
159: PetscCall(PetscMalloc1(2 * nselected_2, &in.pointlist));
160: PetscCall(ISGetIndices(selected_2, &selected_idx_2));
162: for (kk = 0, sid = 0; kk < nselected_2; kk++, sid += 2) {
163: PetscInt lid = selected_idx_2[kk];
164: in.pointlist[sid] = coords[lid];
165: in.pointlist[sid + 1] = coords[data_stride + lid];
166: if (lid >= nFineLoc) nPlotPts++;
167: }
168: PetscCheck(sid == 2 * nselected_2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "sid %d != 2*nselected_2 %" PetscInt_FMT, sid, nselected_2);
170: in.numberofsegments = 0;
171: in.numberofedges = 0;
172: in.numberofholes = 0;
173: in.numberofregions = 0;
174: in.trianglelist = NULL;
175: in.segmentmarkerlist = NULL;
176: in.pointattributelist = NULL;
177: in.pointmarkerlist = NULL;
178: in.triangleattributelist = NULL;
179: in.trianglearealist = NULL;
180: in.segmentlist = NULL;
181: in.holelist = NULL;
182: in.regionlist = NULL;
183: in.edgelist = NULL;
184: in.edgemarkerlist = NULL;
185: in.normlist = NULL;
187: /* triangulate */
188: mid.pointlist = NULL; /* Not needed if -N switch used. */
189: /* Not needed if -N switch used or number of point attributes is zero: */
190: mid.pointattributelist = NULL;
191: mid.pointmarkerlist = NULL; /* Not needed if -N or -B switch used. */
192: mid.trianglelist = NULL; /* Not needed if -E switch used. */
193: /* Not needed if -E switch used or number of triangle attributes is zero: */
194: mid.triangleattributelist = NULL;
195: mid.neighborlist = NULL; /* Needed only if -n switch used. */
196: /* Needed only if segments are output (-p or -c) and -P not used: */
197: mid.segmentlist = NULL;
198: /* Needed only if segments are output (-p or -c) and -P and -B not used: */
199: mid.segmentmarkerlist = NULL;
200: mid.edgelist = NULL; /* Needed only if -e switch used. */
201: mid.edgemarkerlist = NULL; /* Needed if -e used and -B not used. */
202: mid.numberoftriangles = 0;
204: /* Triangulate the points. Switches are chosen to read and write a */
205: /* PSLG (p), preserve the convex hull (c), number everything from */
206: /* zero (z), assign a regional attribute to each element (A), and */
207: /* produce an edge list (e), a Voronoi diagram (v), and a triangle */
208: /* neighbor list (n). */
209: if (nselected_2 != 0) { /* inactive processor */
210: char args[] = "npczQ"; /* c is needed ? */
211: triangulate(args, &in, &mid, (struct triangulateio *)NULL);
212: /* output .poly files for 'showme' */
213: if (!PETSC_TRUE) {
214: static int level = 1;
215: FILE *file;
216: char fname[32];
218: PetscCall(PetscSNPrintf(fname, PETSC_STATIC_ARRAY_LENGTH(fname), "C%d_%d.poly", level, rank));
219: file = fopen(fname, "w");
220: /* First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> */
221: fprintf(file, "%d %d %d %d\n", in.numberofpoints, 2, 0, 0);
222: /* Following lines: <vertex #> <x> <y> */
223: for (kk = 0, sid = 0; kk < in.numberofpoints; kk++, sid += 2) fprintf(file, "%d %e %e\n", kk, in.pointlist[sid], in.pointlist[sid + 1]);
224: /* One line: <# of segments> <# of boundary markers (0 or 1)> */
225: fprintf(file, "%d %d\n", 0, 0);
226: /* Following lines: <segment #> <endpoint> <endpoint> [boundary marker] */
227: /* One line: <# of holes> */
228: fprintf(file, "%d\n", 0);
229: /* Following lines: <hole #> <x> <y> */
230: /* Optional line: <# of regional attributes and/or area constraints> */
231: /* Optional following lines: <region #> <x> <y> <attribute> <maximum area> */
232: fclose(file);
234: /* elems */
235: PetscCall(PetscSNPrintf(fname, PETSC_STATIC_ARRAY_LENGTH(fname), "C%d_%d.ele", level, rank));
236: file = fopen(fname, "w");
237: /* First line: <# of triangles> <nodes per triangle> <# of attributes> */
238: fprintf(file, "%d %d %d\n", mid.numberoftriangles, 3, 0);
239: /* Remaining lines: <triangle #> <node> <node> <node> ... [attributes] */
240: for (kk = 0, sid = 0; kk < mid.numberoftriangles; kk++, sid += 3) fprintf(file, "%d %d %d %d\n", kk, mid.trianglelist[sid], mid.trianglelist[sid + 1], mid.trianglelist[sid + 2]);
241: fclose(file);
243: PetscCall(PetscSNPrintf(fname, PETSC_STATIC_ARRAY_LENGTH(fname), "C%d_%d.node", level, rank));
244: file = fopen(fname, "w");
245: /* First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> */
246: /* fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0); */
247: fprintf(file, "%d %d %d %d\n", nPlotPts, 2, 0, 0);
248: /* Following lines: <vertex #> <x> <y> */
249: for (kk = 0, sid = 0; kk < in.numberofpoints; kk++, sid += 2) fprintf(file, "%d %e %e\n", kk, in.pointlist[sid], in.pointlist[sid + 1]);
251: sid /= 2;
252: for (jj = 0; jj < nFineLoc; jj++) {
253: PetscBool sel = PETSC_TRUE;
254: for (kk = 0; kk < nselected_2 && sel; kk++) {
255: PetscInt lid = selected_idx_2[kk];
256: if (lid == jj) sel = PETSC_FALSE;
257: }
258: if (sel) fprintf(file, "%d %e %e\n", sid++, (double)coords[jj], (double)coords[data_stride + jj]);
259: }
260: fclose(file);
261: PetscCheck(sid == nPlotPts, PETSC_COMM_SELF, PETSC_ERR_PLIB, "sid %d != nPlotPts %d", sid, nPlotPts);
262: level++;
263: }
264: }
265: { /* form P - setup some maps */
266: PetscInt clid, mm, *nTri, *node_tri;
268: PetscCall(PetscMalloc2(nselected_2, &node_tri, nselected_2, &nTri));
270: /* need list of triangles on node */
271: for (kk = 0; kk < nselected_2; kk++) nTri[kk] = 0;
272: for (tid = 0, kk = 0; tid < mid.numberoftriangles; tid++) {
273: for (jj = 0; jj < 3; jj++) {
274: PetscInt cid = mid.trianglelist[kk++];
275: if (nTri[cid] == 0) node_tri[cid] = tid;
276: nTri[cid]++;
277: }
278: }
279: #define EPS 1.e-12
280: /* find points and set prolongation */
281: for (mm = clid = 0; mm < nFineLoc; mm++) {
282: PetscBool ise;
283: PetscCall(PetscCDIsEmptyAt(agg_lists_1, mm, &ise));
284: if (!ise) {
285: const PetscInt lid = mm;
286: PetscScalar AA[3][3];
287: PetscBLASInt N = 3, NRHS = 1, LDA = 3, IPIV[3], LDB = 3;
288: PetscCDIntNd *pos;
290: PetscCall(PetscCDGetHeadPos(agg_lists_1, lid, &pos));
291: while (pos) {
292: PetscInt flid;
293: PetscCall(PetscCDIntNdGetID(pos, &flid));
294: PetscCall(PetscCDGetNextPos(agg_lists_1, lid, &pos));
296: if (flid < nFineLoc) { /* could be a ghost */
297: PetscInt bestTID = -1;
298: PetscReal best_alpha = 1.e10;
299: const PetscInt fgid = flid + myFine0;
300: /* compute shape function for gid */
301: const PetscReal fcoord[3] = {coords[flid], coords[data_stride + flid], 1.0};
302: PetscBool haveit = PETSC_FALSE;
303: PetscScalar alpha[3];
304: PetscInt clids[3];
306: /* look for it */
307: for (tid = node_tri[clid], jj = 0; jj < 5 && !haveit && tid != -1; jj++) {
308: for (tt = 0; tt < 3; tt++) {
309: PetscInt cid2 = mid.trianglelist[3 * tid + tt];
310: PetscInt lid2 = selected_idx_2[cid2];
311: AA[tt][0] = coords[lid2];
312: AA[tt][1] = coords[data_stride + lid2];
313: AA[tt][2] = 1.0;
314: clids[tt] = cid2; /* store for interp */
315: }
317: for (tt = 0; tt < 3; tt++) alpha[tt] = (PetscScalar)fcoord[tt];
319: PetscCallLAPACKInfo("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &info));
320: {
321: PetscBool have = PETSC_TRUE;
322: PetscReal lowest = 1.e10;
323: for (tt = 0, idx = 0; tt < 3; tt++) {
324: if (PetscRealPart(alpha[tt]) > (1.0 + EPS) || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE;
325: if (PetscRealPart(alpha[tt]) < lowest) {
326: lowest = PetscRealPart(alpha[tt]);
327: idx = tt;
328: }
329: }
330: haveit = have;
331: }
332: tid = mid.neighborlist[3 * tid + idx];
333: }
335: if (!haveit) {
336: /* brute force */
337: for (tid = 0; tid < mid.numberoftriangles && !haveit; tid++) {
338: for (tt = 0; tt < 3; tt++) {
339: PetscInt cid2 = mid.trianglelist[3 * tid + tt];
340: PetscInt lid2 = selected_idx_2[cid2];
341: AA[tt][0] = coords[lid2];
342: AA[tt][1] = coords[data_stride + lid2];
343: AA[tt][2] = 1.0;
344: clids[tt] = cid2; /* store for interp */
345: }
346: for (tt = 0; tt < 3; tt++) alpha[tt] = fcoord[tt];
347: PetscCallLAPACKInfo("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &info));
348: {
349: PetscBool have = PETSC_TRUE;
350: PetscReal worst = 0.0, v;
351: for (tt = 0; tt < 3 && have; tt++) {
352: if (PetscRealPart(alpha[tt]) > 1.0 + EPS || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE;
353: if ((v = PetscAbs(PetscRealPart(alpha[tt]) - 0.5)) > worst) worst = v;
354: }
355: if (worst < best_alpha) {
356: best_alpha = worst;
357: bestTID = tid;
358: }
359: haveit = have;
360: }
361: }
362: }
363: if (!haveit) {
364: if (best_alpha > *a_worst_best) *a_worst_best = best_alpha;
365: /* use best one */
366: for (tt = 0; tt < 3; tt++) {
367: PetscInt cid2 = mid.trianglelist[3 * bestTID + tt];
368: PetscInt lid2 = selected_idx_2[cid2];
369: AA[tt][0] = coords[lid2];
370: AA[tt][1] = coords[data_stride + lid2];
371: AA[tt][2] = 1.0;
372: clids[tt] = cid2; /* store for interp */
373: }
374: for (tt = 0; tt < 3; tt++) alpha[tt] = fcoord[tt];
375: PetscCallLAPACKInfo("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &info));
376: }
378: /* put in row of P */
379: for (idx = 0; idx < 3; idx++) {
380: PetscScalar shp = alpha[idx];
381: if (PetscAbs(PetscRealPart(shp)) > 1.e-6) {
382: PetscInt cgid = crsGID[clids[idx]];
383: PetscInt jj = cgid * bs, ii = fgid * bs; /* need to gloalize */
384: for (tt = 0; tt < bs; tt++, ii++, jj++) PetscCall(MatSetValues(a_Prol, 1, &ii, 1, &jj, &shp, INSERT_VALUES));
385: }
386: }
387: }
388: } /* aggregates iterations */
389: clid++;
390: } /* a coarse agg */
391: } /* for all fine nodes */
393: PetscCall(ISRestoreIndices(selected_2, &selected_idx_2));
394: PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
395: PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
397: PetscCall(PetscFree2(node_tri, nTri));
398: }
399: free(mid.trianglelist);
400: free(mid.neighborlist);
401: free(mid.segmentlist);
402: free(mid.segmentmarkerlist);
403: free(mid.pointlist);
404: free(mid.pointmarkerlist);
405: PetscCall(PetscFree(in.pointlist));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: #else
408: SETERRQ(PetscObjectComm((PetscObject)a_Prol), PETSC_ERR_PLIB, "configure with TRIANGLE to use geometric MG");
409: #endif
410: }
412: // PetscClangLinter pragma disable: -fdoc-sowing-chars
413: /*
414: getGIDsOnSquareGraph - square graph, get
416: Input Parameter:
417: . nselected_1 - selected local indices (includes ghosts in input Gmat1)
418: . clid_lid_1 - [nselected_1] lids of selected nodes
419: . Gmat1 - graph that goes with 'selected_1'
420: Output Parameter:
421: . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2)
422: . a_Gmat_2 - graph that is squared of 'Gmat_1'
423: . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes
424: */
425: static PetscErrorCode getGIDsOnSquareGraph(PC pc, PetscInt nselected_1, const PetscInt clid_lid_1[], const Mat Gmat1, IS *a_selected_2, Mat *a_Gmat_2, PetscInt **a_crsGID)
426: {
427: PetscMPIInt size;
428: PetscInt *crsGID, kk, my0, Iend, nloc;
429: MPI_Comm comm;
431: PetscFunctionBegin;
432: PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
433: PetscCallMPI(MPI_Comm_size(comm, &size));
434: PetscCall(MatGetOwnershipRange(Gmat1, &my0, &Iend)); /* AIJ */
435: nloc = Iend - my0; /* this does not change */
437: if (size == 1) { /* not much to do in serial */
438: PetscCall(PetscMalloc1(nselected_1, &crsGID));
439: for (kk = 0; kk < nselected_1; kk++) crsGID[kk] = kk;
440: *a_Gmat_2 = NULL;
441: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nselected_1, clid_lid_1, PETSC_COPY_VALUES, a_selected_2));
442: } else {
443: PetscInt idx, num_fine_ghosts, num_crs_ghost, myCrs0;
444: PetscInt *lid_state, *cpcol_state;
445: Mat_MPIAIJ *mpimat2;
446: Mat Gmat2;
447: PetscSF sf;
449: /* scan my coarse zero gid, set 'lid_state' with coarse GID */
450: kk = nselected_1;
451: PetscCallMPI(MPI_Scan(&kk, &myCrs0, 1, MPIU_INT, MPI_SUM, comm));
452: myCrs0 -= nselected_1;
454: if (a_Gmat_2) { /* output */
455: /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */
456: PetscCall(PCGAMGSquareGraph_GAMG(pc, Gmat1, &Gmat2));
457: *a_Gmat_2 = Gmat2; /* output */
458: } else Gmat2 = Gmat1; /* use local to get crsGIDs at least */
459: /* get coarse grid GIDS for selected (locals and ghosts) */
460: mpimat2 = (Mat_MPIAIJ *)Gmat2->data;
461: PetscCall(VecGetLocalSize(mpimat2->lvec, &num_fine_ghosts));
462: PetscCall(PetscMalloc2(nloc, &lid_state, num_fine_ghosts, &cpcol_state));
463: for (kk = 0; kk < nloc; kk++) lid_state[kk] = -1; /* set with UNKNOWN state */
464: for (kk = 0; kk < nselected_1; kk++) lid_state[clid_lid_1[kk]] = kk + myCrs0; /* set with coarse GID */
465: PetscCall(MatGetMultPetscSF(Gmat2, &sf));
466: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE));
467: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE));
468: for (kk = 0, num_crs_ghost = 0; kk < num_fine_ghosts; kk++) {
469: if (cpcol_state[kk] != -1) num_crs_ghost++;
470: }
471: PetscCall(PetscMalloc1(nselected_1 + num_crs_ghost, &crsGID)); /* output */
472: {
473: PetscInt *selected_set;
474: PetscCall(PetscMalloc1(nselected_1 + num_crs_ghost, &selected_set));
475: /* do ghost of 'crsGID' */
476: for (kk = 0, idx = nselected_1; kk < num_fine_ghosts; kk++) {
477: if (cpcol_state[kk] != -1) {
478: PetscInt cgid = cpcol_state[kk];
479: selected_set[idx] = nloc + kk;
480: crsGID[idx++] = cgid;
481: }
482: }
483: PetscCheck(idx == (nselected_1 + num_crs_ghost), PETSC_COMM_SELF, PETSC_ERR_PLIB, "idx %" PetscInt_FMT " != (nselected_1 %" PetscInt_FMT " + num_crs_ghost %" PetscInt_FMT ")", idx, nselected_1, num_crs_ghost);
484: /* do locals in 'crsGID' */
485: for (kk = 0, idx = 0; kk < nloc; kk++) {
486: if (lid_state[kk] != -1) {
487: PetscInt cgid = lid_state[kk];
488: selected_set[idx] = kk;
489: crsGID[idx++] = cgid;
490: }
491: }
492: PetscCheck(idx == nselected_1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "idx %" PetscInt_FMT " != nselected_1 %" PetscInt_FMT, idx, nselected_1);
494: if (a_selected_2 != NULL) { /* output */
495: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nselected_1 + num_crs_ghost, selected_set, PETSC_OWN_POINTER, a_selected_2));
496: } else {
497: PetscCall(PetscFree(selected_set));
498: }
499: }
500: PetscCall(PetscFree2(lid_state, cpcol_state));
501: }
502: *a_crsGID = crsGID; /* output */
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode PCGAMGCreateGraph_GEO(PC pc, Mat Amat, Mat *a_Gmat)
507: {
508: PC_MG *mg = (PC_MG *)pc->data;
509: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
510: const PetscReal vfilter = pc_gamg->threshold[0];
512: PetscFunctionBegin;
513: PetscCall(MatCreateGraph(Amat, PETSC_TRUE, PETSC_TRUE, vfilter, 0, NULL, a_Gmat));
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: static PetscErrorCode PCGAMGCoarsen_GEO(PC a_pc, Mat *a_Gmat, PetscCoarsenData **a_llist_parent)
518: {
519: PetscInt Istart, Iend, nloc, kk, Ii, ncols;
520: IS perm;
521: GAMGNode *gnodes;
522: PetscInt *permute;
523: Mat Gmat = *a_Gmat;
524: MPI_Comm comm;
525: MatCoarsen crs;
527: PetscFunctionBegin;
528: PetscCall(PetscObjectGetComm((PetscObject)a_pc, &comm));
530: PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
531: nloc = (Iend - Istart);
533: /* create random permutation with sort for geo-mg */
534: PetscCall(PetscMalloc1(nloc, &gnodes));
535: PetscCall(PetscMalloc1(nloc, &permute));
537: for (Ii = Istart; Ii < Iend; Ii++) { /* locals only? */
538: PetscCall(MatGetRow(Gmat, Ii, &ncols, NULL, NULL));
539: {
540: PetscInt lid = Ii - Istart;
541: gnodes[lid].lid = lid;
542: gnodes[lid].degree = ncols;
543: }
544: PetscCall(MatRestoreRow(Gmat, Ii, &ncols, NULL, NULL));
545: }
546: if (PETSC_TRUE) {
547: PetscRandom rand;
548: PetscBool *bIndexSet;
549: PetscReal rr;
550: PetscInt iSwapIndex;
552: PetscCall(PetscRandomCreate(comm, &rand));
553: PetscCall(PetscCalloc1(nloc, &bIndexSet));
554: for (Ii = 0; Ii < nloc; Ii++) {
555: PetscCall(PetscRandomGetValueReal(rand, &rr));
556: iSwapIndex = (PetscInt)(rr * nloc);
557: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
558: GAMGNode iTemp = gnodes[iSwapIndex];
559: gnodes[iSwapIndex] = gnodes[Ii];
560: gnodes[Ii] = iTemp;
561: bIndexSet[Ii] = PETSC_TRUE;
562: bIndexSet[iSwapIndex] = PETSC_TRUE;
563: }
564: }
565: PetscCall(PetscRandomDestroy(&rand));
566: PetscCall(PetscFree(bIndexSet));
567: }
568: /* only sort locals */
569: if (gnodes) qsort(gnodes, nloc, sizeof(GAMGNode), petsc_geo_mg_compare);
570: /* create IS of permutation */
571: for (kk = 0; kk < nloc; kk++) permute[kk] = gnodes[kk].lid; /* locals only */
572: PetscCall(PetscFree(gnodes));
573: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_OWN_POINTER, &perm));
575: /* get MIS aggs */
577: PetscCall(MatCoarsenCreate(comm, &crs));
578: PetscCall(MatCoarsenSetType(crs, MATCOARSENMIS));
579: PetscCall(MatCoarsenSetGreedyOrdering(crs, perm));
580: PetscCall(MatCoarsenSetAdjacency(crs, Gmat));
581: PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_FALSE));
582: PetscCall(MatCoarsenApply(crs));
583: PetscCall(MatCoarsenGetData(crs, a_llist_parent));
584: PetscCall(MatCoarsenDestroy(&crs));
586: PetscCall(ISDestroy(&perm));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: static PetscErrorCode PCGAMGProlongator_GEO(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
591: {
592: PC_MG *mg = (PC_MG *)pc->data;
593: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
594: const PetscInt dim = pc_gamg->data_cell_cols, data_cols = pc_gamg->data_cell_cols;
595: PetscInt Istart, Iend, nloc, my0, jj, kk, ncols, nLocalSelected, bs, *clid_flid;
596: Mat Prol, Gmat;
597: PetscMPIInt rank, size;
598: MPI_Comm comm;
599: IS selected_2, selected_1;
600: const PetscInt *selected_idx;
601: MatType mtype;
603: PetscFunctionBegin;
604: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
606: PetscCallMPI(MPI_Comm_rank(comm, &rank));
607: PetscCallMPI(MPI_Comm_size(comm, &size));
608: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
609: PetscCall(MatGetBlockSize(Amat, &bs));
610: nloc = (Iend - Istart) / bs;
611: my0 = Istart / bs;
612: PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") %% bs %" PetscInt_FMT, Iend, Istart, bs);
614: /* get 'nLocalSelected' */
615: PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges
616: PetscCall(PetscCDGetNonemptyIS(agg_lists, &selected_1));
617: PetscCall(ISGetSize(selected_1, &jj));
618: PetscCall(PetscMalloc1(jj, &clid_flid));
619: PetscCall(ISGetIndices(selected_1, &selected_idx));
620: for (kk = 0, nLocalSelected = 0; kk < jj; kk++) {
621: PetscInt lid = selected_idx[kk];
622: if (lid < nloc) {
623: PetscCall(MatGetRow(Gmat, lid + my0, &ncols, NULL, NULL));
624: if (ncols > 1) clid_flid[nLocalSelected++] = lid; /* filter out singletons */
625: PetscCall(MatRestoreRow(Gmat, lid + my0, &ncols, NULL, NULL));
626: }
627: }
628: PetscCall(ISRestoreIndices(selected_1, &selected_idx));
629: PetscCall(ISDestroy(&selected_1)); /* this is selected_1 in serial */
631: /* create prolongator matrix */
632: PetscCall(MatGetType(Amat, &mtype));
633: PetscCall(MatCreate(comm, &Prol));
634: PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * bs, PETSC_DETERMINE, PETSC_DETERMINE));
635: PetscCall(MatSetBlockSizes(Prol, bs, bs));
636: PetscCall(MatSetType(Prol, mtype));
637: PetscCall(MatSeqAIJSetPreallocation(Prol, 3 * data_cols, NULL));
638: PetscCall(MatMPIAIJSetPreallocation(Prol, 3 * data_cols, NULL, 3 * data_cols, NULL));
640: /* can get all points "removed" - but not on geomg */
641: PetscCall(MatGetSize(Prol, &kk, &jj));
642: if (!jj) {
643: PetscCall(PetscInfo(pc, "ERROE: no selected points on coarse grid\n"));
644: PetscCall(PetscFree(clid_flid));
645: PetscCall(MatDestroy(&Prol));
646: *a_P_out = NULL; /* out */
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: {
651: PetscReal *coords;
652: PetscInt data_stride;
653: PetscInt *crsGID = NULL;
654: Mat Gmat2;
656: PetscCheck(dim == data_cols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dim %" PetscInt_FMT " != data_cols %" PetscInt_FMT, dim, data_cols);
657: /* grow ghost data for better coarse grid cover of fine grid */
658: /* messy method, squares graph and gets some data */
659: PetscCall(getGIDsOnSquareGraph(pc, nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID));
660: /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */
661: /* create global vector of coorindates in 'coords' */
662: if (size > 1) {
663: PetscCall(PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords));
664: } else {
665: coords = pc_gamg->data;
666: data_stride = pc_gamg->data_sz / pc_gamg->data_cell_cols;
667: }
668: PetscCall(MatDestroy(&Gmat2));
669: /* triangulate */
670: {
671: PetscReal metric, tm;
673: PetscCheck(dim == 2, comm, PETSC_ERR_PLIB, "3D not implemented for 'geo' AMG");
674: PetscCall(triangulateAndFormProl(selected_2, data_stride, coords, nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric));
675: PetscCall(PetscFree(crsGID));
677: /* clean up and create coordinates for coarse grid (output) */
678: if (size > 1) PetscCall(PetscFree(coords));
680: PetscCallMPI(MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm));
681: if (tm > 1.) { /* needs to be globalized - should not happen */
682: PetscCall(PetscInfo(pc, " failed metric for coarse grid %e\n", (double)tm));
683: PetscCall(MatDestroy(&Prol));
684: } else if (metric > .0) {
685: PetscCall(PetscInfo(pc, "worst metric for coarse grid = %e\n", (double)metric));
686: }
687: }
688: { /* create next coords - output */
689: PetscReal *crs_crds;
690: PetscCall(PetscMalloc1(dim * nLocalSelected, &crs_crds));
691: for (kk = 0; kk < nLocalSelected; kk++) { /* grab local select nodes to promote - output */
692: PetscInt lid = clid_flid[kk];
693: for (jj = 0; jj < dim; jj++) crs_crds[jj * nLocalSelected + kk] = pc_gamg->data[jj * nloc + lid];
694: }
696: PetscCall(PetscFree(pc_gamg->data));
697: pc_gamg->data = crs_crds; /* out */
698: pc_gamg->data_sz = dim * nLocalSelected;
699: }
700: PetscCall(ISDestroy(&selected_2));
701: }
703: *a_P_out = Prol; /* out */
704: PetscCall(PetscFree(clid_flid));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode PCDestroy_GAMG_GEO(PC pc)
709: {
710: PetscFunctionBegin;
711: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: PetscErrorCode PCCreateGAMG_GEO(PC pc)
716: {
717: PC_MG *mg = (PC_MG *)pc->data;
718: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
720: PetscFunctionBegin;
721: pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO;
722: pc_gamg->ops->destroy = PCDestroy_GAMG_GEO;
723: /* reset does not do anything; setup not virtual */
725: /* set internal function pointers */
726: pc_gamg->ops->creategraph = PCGAMGCreateGraph_GEO;
727: pc_gamg->ops->coarsen = PCGAMGCoarsen_GEO;
728: pc_gamg->ops->prolongator = PCGAMGProlongator_GEO;
729: pc_gamg->ops->optprolongator = NULL;
730: pc_gamg->ops->createdefaultdata = PCSetData_GEO;
732: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_GEO));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }