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