1 /*
2 GAMG geometric-algebraic multigrid PC - Mark Adams 2011
3 */
4
5 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/
6
7 #if defined(PETSC_HAVE_TRIANGLE)
8 #if !defined(ANSI_DECLARATORS)
9 #define ANSI_DECLARATORS
10 #endif
11 #include <triangle.h>
12 #endif
13
14 #include <petscblaslapack.h>
15
16 /* Private context for the GAMG preconditioner */
17 typedef struct {
18 PetscInt lid; /* local vertex index */
19 PetscInt degree; /* vertex degree */
20 } GAMGNode;
21
petsc_geo_mg_compare(const void * a,const void * b)22 static inline int petsc_geo_mg_compare(const void *a, const void *b)
23 {
24 return (int)(((GAMGNode *)a)->degree - ((GAMGNode *)b)->degree);
25 }
26
27 // PetscClangLinter pragma disable: -fdoc-sowing-chars
28 /*
29 PCSetCoordinates_GEO
30
31 Input Parameter:
32 . pc - the preconditioner context
33 */
PCSetCoordinates_GEO(PC pc,PetscInt ndm,PetscInt a_nloc,PetscReal * coords)34 static 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;
38 PetscInt arrsz, bs, my0, kk, ii, nloc, Iend, aloc;
39 Mat Amat = pc->pmat;
40
41 PetscFunctionBegin;
42 PetscValidHeaderSpecific(Amat, MAT_CLASSID, 1);
43 PetscCall(MatGetBlockSize(Amat, &bs));
44 PetscCall(MatGetOwnershipRange(Amat, &my0, &Iend));
45 aloc = (Iend - my0);
46 nloc = (Iend - my0) / bs;
47
48 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
50 pc_gamg->data_cell_rows = 1;
51 PetscCheck(coords || (nloc <= 0), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need coordinates for pc_gamg_type 'geo'.");
52 pc_gamg->data_cell_cols = ndm; /* coordinates */
53
54 arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;
55
56 /* create data - syntactic sugar that should be refactored at some point */
57 if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
58 PetscCall(PetscFree(pc_gamg->data));
59 PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
60 }
61 for (kk = 0; kk < arrsz; kk++) pc_gamg->data[kk] = -999.;
62 pc_gamg->data[arrsz] = -99.;
63 /* copy data in - column-oriented */
64 if (nloc == a_nloc) {
65 for (kk = 0; kk < nloc; kk++) {
66 for (ii = 0; ii < ndm; ii++) pc_gamg->data[ii * nloc + kk] = coords[kk * ndm + ii];
67 }
68 } else { /* assumes the coordinates are blocked */
69 for (kk = 0; kk < nloc; kk++) {
70 for (ii = 0; ii < ndm; ii++) pc_gamg->data[ii * nloc + kk] = coords[bs * kk * ndm + ii];
71 }
72 }
73 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]);
74 pc_gamg->data_sz = arrsz;
75 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
78 // PetscClangLinter pragma disable: -fdoc-sowing-chars
79 /*
80 PCSetData_GEO
81
82 Input Parameter:
83 . pc -
84 */
PCSetData_GEO(PC pc,Mat m)85 static PetscErrorCode PCSetData_GEO(PC pc, Mat m)
86 {
87 PetscFunctionBegin;
88 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "GEO MG needs coordinates");
89 }
90
PCSetFromOptions_GEO(PC pc,PetscOptionItems PetscOptionsObject)91 static PetscErrorCode PCSetFromOptions_GEO(PC pc, PetscOptionItems PetscOptionsObject)
92 {
93 PetscFunctionBegin;
94 PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-GEO options");
95 {
96 /* -pc_gamg_sa_nsmooths */
97 /* pc_gamg_sa->smooths = 0; */
98 /* ierr = PetscOptionsInt("-pc_gamg_agg_nsmooths", */
99 /* "smoothing steps for smoothed aggregation, usually 1 (0)", */
100 /* "PCGAMGSetNSmooths_AGG", */
101 /* pc_gamg_sa->smooths, */
102 /* &pc_gamg_sa->smooths, */
103 /* &flag); */
104 }
105 PetscOptionsHeadEnd();
106 PetscFunctionReturn(PETSC_SUCCESS);
107 }
108
109 // PetscClangLinter pragma disable: -fdoc-sowing-chars
110 /*
111 triangulateAndFormProl
112
113 Input Parameter:
114 . selected_2 - list of selected local ID, includes selected ghosts
115 . data_stride -
116 . coords[2*data_stride] - column vector of local coordinates w/ ghosts
117 . nselected_1 - selected IDs that go with base (1) graph includes selected ghosts
118 . clid_lid_1[nselected_1] - lids of selected (c) nodes ???????????
119 . agg_lists_1 - list of aggregates selected_1 vertices of aggregate unselected vertices
120 . crsGID[selected.size()] - global index for prolongation operator
121 . bs - block size
122 Output Parameter:
123 . a_Prol - prolongation operator
124 . a_worst_best - measure of worst missed fine vertex, 0 is no misses
125 */
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)126 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)
127 {
128 #if defined(PETSC_HAVE_TRIANGLE)
129 PetscInt jj, tid, tt, idx, nselected_2;
130 struct triangulateio in, mid;
131 const PetscInt *selected_idx_2;
132 PetscMPIInt rank;
133 PetscInt Istart, Iend, nFineLoc, myFine0;
134 int kk, nPlotPts, sid;
135 MPI_Comm comm;
136 PetscReal tm;
137
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(a_worst_best, &tm, 1, MPIU_REAL, MPIU_MAX, comm));
146 if (tm > 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));
161
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);
169
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;
186
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;
203
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];
217
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);
233
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);
242
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]);
250
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;
267
268 PetscCall(PetscMalloc2(nselected_2, &node_tri, nselected_2, &nTri));
269
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, INFO;
288 PetscCDIntNd *pos;
289
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));
295
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];
305
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 }
316
317 for (tt = 0; tt < 3; tt++) alpha[tt] = (PetscScalar)fcoord[tt];
318
319 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
320 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &INFO));
321 {
322 PetscBool have = PETSC_TRUE;
323 PetscReal lowest = 1.e10;
324 for (tt = 0, idx = 0; tt < 3; tt++) {
325 if (PetscRealPart(alpha[tt]) > (1.0 + EPS) || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE;
326 if (PetscRealPart(alpha[tt]) < lowest) {
327 lowest = PetscRealPart(alpha[tt]);
328 idx = tt;
329 }
330 }
331 haveit = have;
332 }
333 tid = mid.neighborlist[3 * tid + idx];
334 }
335
336 if (!haveit) {
337 /* brute force */
338 for (tid = 0; tid < mid.numberoftriangles && !haveit; tid++) {
339 for (tt = 0; tt < 3; tt++) {
340 PetscInt cid2 = mid.trianglelist[3 * tid + tt];
341 PetscInt lid2 = selected_idx_2[cid2];
342 AA[tt][0] = coords[lid2];
343 AA[tt][1] = coords[data_stride + lid2];
344 AA[tt][2] = 1.0;
345 clids[tt] = cid2; /* store for interp */
346 }
347 for (tt = 0; tt < 3; tt++) alpha[tt] = fcoord[tt];
348 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
349 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &INFO));
350 {
351 PetscBool have = PETSC_TRUE;
352 PetscReal worst = 0.0, v;
353 for (tt = 0; tt < 3 && have; tt++) {
354 if (PetscRealPart(alpha[tt]) > 1.0 + EPS || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE;
355 if ((v = PetscAbs(PetscRealPart(alpha[tt]) - 0.5)) > worst) worst = v;
356 }
357 if (worst < best_alpha) {
358 best_alpha = worst;
359 bestTID = tid;
360 }
361 haveit = have;
362 }
363 }
364 }
365 if (!haveit) {
366 if (best_alpha > *a_worst_best) *a_worst_best = best_alpha;
367 /* use best one */
368 for (tt = 0; tt < 3; tt++) {
369 PetscInt cid2 = mid.trianglelist[3 * bestTID + tt];
370 PetscInt lid2 = selected_idx_2[cid2];
371 AA[tt][0] = coords[lid2];
372 AA[tt][1] = coords[data_stride + lid2];
373 AA[tt][2] = 1.0;
374 clids[tt] = cid2; /* store for interp */
375 }
376 for (tt = 0; tt < 3; tt++) alpha[tt] = fcoord[tt];
377 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
378 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &INFO));
379 }
380
381 /* put in row of P */
382 for (idx = 0; idx < 3; idx++) {
383 PetscScalar shp = alpha[idx];
384 if (PetscAbs(PetscRealPart(shp)) > 1.e-6) {
385 PetscInt cgid = crsGID[clids[idx]];
386 PetscInt jj = cgid * bs, ii = fgid * bs; /* need to gloalize */
387 for (tt = 0; tt < bs; tt++, ii++, jj++) PetscCall(MatSetValues(a_Prol, 1, &ii, 1, &jj, &shp, INSERT_VALUES));
388 }
389 }
390 }
391 } /* aggregates iterations */
392 clid++;
393 } /* a coarse agg */
394 } /* for all fine nodes */
395
396 PetscCall(ISRestoreIndices(selected_2, &selected_idx_2));
397 PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
398 PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
399
400 PetscCall(PetscFree2(node_tri, nTri));
401 }
402 free(mid.trianglelist);
403 free(mid.neighborlist);
404 free(mid.segmentlist);
405 free(mid.segmentmarkerlist);
406 free(mid.pointlist);
407 free(mid.pointmarkerlist);
408 PetscCall(PetscFree(in.pointlist));
409 PetscFunctionReturn(PETSC_SUCCESS);
410 #else
411 SETERRQ(PetscObjectComm((PetscObject)a_Prol), PETSC_ERR_PLIB, "configure with TRIANGLE to use geometric MG");
412 #endif
413 }
414
415 // PetscClangLinter pragma disable: -fdoc-sowing-chars
416 /*
417 getGIDsOnSquareGraph - square graph, get
418
419 Input Parameter:
420 . nselected_1 - selected local indices (includes ghosts in input Gmat1)
421 . clid_lid_1 - [nselected_1] lids of selected nodes
422 . Gmat1 - graph that goes with 'selected_1'
423 Output Parameter:
424 . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2)
425 . a_Gmat_2 - graph that is squared of 'Gmat_1'
426 . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes
427 */
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)428 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)
429 {
430 PetscMPIInt size;
431 PetscInt *crsGID, kk, my0, Iend, nloc;
432 MPI_Comm comm;
433
434 PetscFunctionBegin;
435 PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
436 PetscCallMPI(MPI_Comm_size(comm, &size));
437 PetscCall(MatGetOwnershipRange(Gmat1, &my0, &Iend)); /* AIJ */
438 nloc = Iend - my0; /* this does not change */
439
440 if (size == 1) { /* not much to do in serial */
441 PetscCall(PetscMalloc1(nselected_1, &crsGID));
442 for (kk = 0; kk < nselected_1; kk++) crsGID[kk] = kk;
443 *a_Gmat_2 = NULL;
444 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nselected_1, clid_lid_1, PETSC_COPY_VALUES, a_selected_2));
445 } else {
446 PetscInt idx, num_fine_ghosts, num_crs_ghost, myCrs0;
447 Mat_MPIAIJ *mpimat2;
448 Mat Gmat2;
449 Vec locState;
450 PetscScalar *cpcol_state;
451
452 /* scan my coarse zero gid, set 'lid_state' with coarse GID */
453 kk = nselected_1;
454 PetscCallMPI(MPI_Scan(&kk, &myCrs0, 1, MPIU_INT, MPI_SUM, comm));
455 myCrs0 -= nselected_1;
456
457 if (a_Gmat_2) { /* output */
458 /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */
459 PetscCall(PCGAMGSquareGraph_GAMG(pc, Gmat1, &Gmat2));
460 *a_Gmat_2 = Gmat2; /* output */
461 } else Gmat2 = Gmat1; /* use local to get crsGIDs at least */
462 /* get coarse grid GIDS for selected (locals and ghosts) */
463 mpimat2 = (Mat_MPIAIJ *)Gmat2->data;
464 PetscCall(MatCreateVecs(Gmat2, &locState, NULL));
465 PetscCall(VecSet(locState, -1)); /* set with UNKNOWN state */
466 for (kk = 0; kk < nselected_1; kk++) {
467 PetscInt fgid = clid_lid_1[kk] + my0;
468 PetscScalar v = (PetscScalar)(kk + myCrs0);
469 PetscCall(VecSetValues(locState, 1, &fgid, &v, INSERT_VALUES)); /* set with PID */
470 }
471 PetscCall(VecAssemblyBegin(locState));
472 PetscCall(VecAssemblyEnd(locState));
473 PetscCall(VecScatterBegin(mpimat2->Mvctx, locState, mpimat2->lvec, INSERT_VALUES, SCATTER_FORWARD));
474 PetscCall(VecScatterEnd(mpimat2->Mvctx, locState, mpimat2->lvec, INSERT_VALUES, SCATTER_FORWARD));
475 PetscCall(VecGetLocalSize(mpimat2->lvec, &num_fine_ghosts));
476 PetscCall(VecGetArray(mpimat2->lvec, &cpcol_state));
477 for (kk = 0, num_crs_ghost = 0; kk < num_fine_ghosts; kk++) {
478 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) num_crs_ghost++;
479 }
480 PetscCall(PetscMalloc1(nselected_1 + num_crs_ghost, &crsGID)); /* output */
481 {
482 PetscInt *selected_set;
483 PetscCall(PetscMalloc1(nselected_1 + num_crs_ghost, &selected_set));
484 /* do ghost of 'crsGID' */
485 for (kk = 0, idx = nselected_1; kk < num_fine_ghosts; kk++) {
486 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
487 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
488 selected_set[idx] = nloc + kk;
489 crsGID[idx++] = cgid;
490 }
491 }
492 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);
493 PetscCall(VecRestoreArray(mpimat2->lvec, &cpcol_state));
494 /* do locals in 'crsGID' */
495 PetscCall(VecGetArray(locState, &cpcol_state));
496 for (kk = 0, idx = 0; kk < nloc; kk++) {
497 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
498 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
499 selected_set[idx] = kk;
500 crsGID[idx++] = cgid;
501 }
502 }
503 PetscCheck(idx == nselected_1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "idx %" PetscInt_FMT " != nselected_1 %" PetscInt_FMT, idx, nselected_1);
504 PetscCall(VecRestoreArray(locState, &cpcol_state));
505
506 if (a_selected_2 != NULL) { /* output */
507 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nselected_1 + num_crs_ghost, selected_set, PETSC_OWN_POINTER, a_selected_2));
508 } else {
509 PetscCall(PetscFree(selected_set));
510 }
511 }
512 PetscCall(VecDestroy(&locState));
513 }
514 *a_crsGID = crsGID; /* output */
515 PetscFunctionReturn(PETSC_SUCCESS);
516 }
517
PCGAMGCreateGraph_GEO(PC pc,Mat Amat,Mat * a_Gmat)518 static PetscErrorCode PCGAMGCreateGraph_GEO(PC pc, Mat Amat, Mat *a_Gmat)
519 {
520 PC_MG *mg = (PC_MG *)pc->data;
521 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
522 const PetscReal vfilter = pc_gamg->threshold[0];
523
524 PetscFunctionBegin;
525 PetscCall(MatCreateGraph(Amat, PETSC_TRUE, PETSC_TRUE, vfilter, 0, NULL, a_Gmat));
526 PetscFunctionReturn(PETSC_SUCCESS);
527 }
528
PCGAMGCoarsen_GEO(PC a_pc,Mat * a_Gmat,PetscCoarsenData ** a_llist_parent)529 static PetscErrorCode PCGAMGCoarsen_GEO(PC a_pc, Mat *a_Gmat, PetscCoarsenData **a_llist_parent)
530 {
531 PetscInt Istart, Iend, nloc, kk, Ii, ncols;
532 IS perm;
533 GAMGNode *gnodes;
534 PetscInt *permute;
535 Mat Gmat = *a_Gmat;
536 MPI_Comm comm;
537 MatCoarsen crs;
538
539 PetscFunctionBegin;
540 PetscCall(PetscObjectGetComm((PetscObject)a_pc, &comm));
541
542 PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
543 nloc = (Iend - Istart);
544
545 /* create random permutation with sort for geo-mg */
546 PetscCall(PetscMalloc1(nloc, &gnodes));
547 PetscCall(PetscMalloc1(nloc, &permute));
548
549 for (Ii = Istart; Ii < Iend; Ii++) { /* locals only? */
550 PetscCall(MatGetRow(Gmat, Ii, &ncols, NULL, NULL));
551 {
552 PetscInt lid = Ii - Istart;
553 gnodes[lid].lid = lid;
554 gnodes[lid].degree = ncols;
555 }
556 PetscCall(MatRestoreRow(Gmat, Ii, &ncols, NULL, NULL));
557 }
558 if (PETSC_TRUE) {
559 PetscRandom rand;
560 PetscBool *bIndexSet;
561 PetscReal rr;
562 PetscInt iSwapIndex;
563
564 PetscCall(PetscRandomCreate(comm, &rand));
565 PetscCall(PetscCalloc1(nloc, &bIndexSet));
566 for (Ii = 0; Ii < nloc; Ii++) {
567 PetscCall(PetscRandomGetValueReal(rand, &rr));
568 iSwapIndex = (PetscInt)(rr * nloc);
569 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
570 GAMGNode iTemp = gnodes[iSwapIndex];
571 gnodes[iSwapIndex] = gnodes[Ii];
572 gnodes[Ii] = iTemp;
573 bIndexSet[Ii] = PETSC_TRUE;
574 bIndexSet[iSwapIndex] = PETSC_TRUE;
575 }
576 }
577 PetscCall(PetscRandomDestroy(&rand));
578 PetscCall(PetscFree(bIndexSet));
579 }
580 /* only sort locals */
581 if (gnodes) qsort(gnodes, nloc, sizeof(GAMGNode), petsc_geo_mg_compare);
582 /* create IS of permutation */
583 for (kk = 0; kk < nloc; kk++) permute[kk] = gnodes[kk].lid; /* locals only */
584 PetscCall(PetscFree(gnodes));
585 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_OWN_POINTER, &perm));
586
587 /* get MIS aggs */
588
589 PetscCall(MatCoarsenCreate(comm, &crs));
590 PetscCall(MatCoarsenSetType(crs, MATCOARSENMIS));
591 PetscCall(MatCoarsenSetGreedyOrdering(crs, perm));
592 PetscCall(MatCoarsenSetAdjacency(crs, Gmat));
593 PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_FALSE));
594 PetscCall(MatCoarsenApply(crs));
595 PetscCall(MatCoarsenGetData(crs, a_llist_parent));
596 PetscCall(MatCoarsenDestroy(&crs));
597
598 PetscCall(ISDestroy(&perm));
599 PetscFunctionReturn(PETSC_SUCCESS);
600 }
601
PCGAMGProlongator_GEO(PC pc,Mat Amat,PetscCoarsenData * agg_lists,Mat * a_P_out)602 static PetscErrorCode PCGAMGProlongator_GEO(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
603 {
604 PC_MG *mg = (PC_MG *)pc->data;
605 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
606 const PetscInt dim = pc_gamg->data_cell_cols, data_cols = pc_gamg->data_cell_cols;
607 PetscInt Istart, Iend, nloc, my0, jj, kk, ncols, nLocalSelected, bs, *clid_flid;
608 Mat Prol, Gmat;
609 PetscMPIInt rank, size;
610 MPI_Comm comm;
611 IS selected_2, selected_1;
612 const PetscInt *selected_idx;
613 MatType mtype;
614
615 PetscFunctionBegin;
616 PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
617
618 PetscCallMPI(MPI_Comm_rank(comm, &rank));
619 PetscCallMPI(MPI_Comm_size(comm, &size));
620 PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
621 PetscCall(MatGetBlockSize(Amat, &bs));
622 nloc = (Iend - Istart) / bs;
623 my0 = Istart / bs;
624 PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") %% bs %" PetscInt_FMT, Iend, Istart, bs);
625
626 /* get 'nLocalSelected' */
627 PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges
628 PetscCall(PetscCDGetNonemptyIS(agg_lists, &selected_1));
629 PetscCall(ISGetSize(selected_1, &jj));
630 PetscCall(PetscMalloc1(jj, &clid_flid));
631 PetscCall(ISGetIndices(selected_1, &selected_idx));
632 for (kk = 0, nLocalSelected = 0; kk < jj; kk++) {
633 PetscInt lid = selected_idx[kk];
634 if (lid < nloc) {
635 PetscCall(MatGetRow(Gmat, lid + my0, &ncols, NULL, NULL));
636 if (ncols > 1) clid_flid[nLocalSelected++] = lid; /* filter out singletons */
637 PetscCall(MatRestoreRow(Gmat, lid + my0, &ncols, NULL, NULL));
638 }
639 }
640 PetscCall(ISRestoreIndices(selected_1, &selected_idx));
641 PetscCall(ISDestroy(&selected_1)); /* this is selected_1 in serial */
642
643 /* create prolongator matrix */
644 PetscCall(MatGetType(Amat, &mtype));
645 PetscCall(MatCreate(comm, &Prol));
646 PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * bs, PETSC_DETERMINE, PETSC_DETERMINE));
647 PetscCall(MatSetBlockSizes(Prol, bs, bs));
648 PetscCall(MatSetType(Prol, mtype));
649 PetscCall(MatSeqAIJSetPreallocation(Prol, 3 * data_cols, NULL));
650 PetscCall(MatMPIAIJSetPreallocation(Prol, 3 * data_cols, NULL, 3 * data_cols, NULL));
651
652 /* can get all points "removed" - but not on geomg */
653 PetscCall(MatGetSize(Prol, &kk, &jj));
654 if (!jj) {
655 PetscCall(PetscInfo(pc, "ERROE: no selected points on coarse grid\n"));
656 PetscCall(PetscFree(clid_flid));
657 PetscCall(MatDestroy(&Prol));
658 *a_P_out = NULL; /* out */
659 PetscFunctionReturn(PETSC_SUCCESS);
660 }
661
662 {
663 PetscReal *coords;
664 PetscInt data_stride;
665 PetscInt *crsGID = NULL;
666 Mat Gmat2;
667
668 PetscCheck(dim == data_cols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dim %" PetscInt_FMT " != data_cols %" PetscInt_FMT, dim, data_cols);
669 /* grow ghost data for better coarse grid cover of fine grid */
670 /* messy method, squares graph and gets some data */
671 PetscCall(getGIDsOnSquareGraph(pc, nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID));
672 /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */
673 /* create global vector of coorindates in 'coords' */
674 if (size > 1) {
675 PetscCall(PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords));
676 } else {
677 coords = pc_gamg->data;
678 data_stride = pc_gamg->data_sz / pc_gamg->data_cell_cols;
679 }
680 PetscCall(MatDestroy(&Gmat2));
681 /* triangulate */
682 {
683 PetscReal metric, tm;
684
685 PetscCheck(dim == 2, comm, PETSC_ERR_PLIB, "3D not implemented for 'geo' AMG");
686 PetscCall(triangulateAndFormProl(selected_2, data_stride, coords, nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric));
687 PetscCall(PetscFree(crsGID));
688
689 /* clean up and create coordinates for coarse grid (output) */
690 if (size > 1) PetscCall(PetscFree(coords));
691
692 PetscCallMPI(MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm));
693 if (tm > 1.) { /* needs to be globalized - should not happen */
694 PetscCall(PetscInfo(pc, " failed metric for coarse grid %e\n", (double)tm));
695 PetscCall(MatDestroy(&Prol));
696 } else if (metric > .0) {
697 PetscCall(PetscInfo(pc, "worst metric for coarse grid = %e\n", (double)metric));
698 }
699 }
700 { /* create next coords - output */
701 PetscReal *crs_crds;
702 PetscCall(PetscMalloc1(dim * nLocalSelected, &crs_crds));
703 for (kk = 0; kk < nLocalSelected; kk++) { /* grab local select nodes to promote - output */
704 PetscInt lid = clid_flid[kk];
705 for (jj = 0; jj < dim; jj++) crs_crds[jj * nLocalSelected + kk] = pc_gamg->data[jj * nloc + lid];
706 }
707
708 PetscCall(PetscFree(pc_gamg->data));
709 pc_gamg->data = crs_crds; /* out */
710 pc_gamg->data_sz = dim * nLocalSelected;
711 }
712 PetscCall(ISDestroy(&selected_2));
713 }
714
715 *a_P_out = Prol; /* out */
716 PetscCall(PetscFree(clid_flid));
717 PetscFunctionReturn(PETSC_SUCCESS);
718 }
719
PCDestroy_GAMG_GEO(PC pc)720 static PetscErrorCode PCDestroy_GAMG_GEO(PC pc)
721 {
722 PetscFunctionBegin;
723 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
724 PetscFunctionReturn(PETSC_SUCCESS);
725 }
726
PCCreateGAMG_GEO(PC pc)727 PetscErrorCode PCCreateGAMG_GEO(PC pc)
728 {
729 PC_MG *mg = (PC_MG *)pc->data;
730 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
731
732 PetscFunctionBegin;
733 pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO;
734 pc_gamg->ops->destroy = PCDestroy_GAMG_GEO;
735 /* reset does not do anything; setup not virtual */
736
737 /* set internal function pointers */
738 pc_gamg->ops->creategraph = PCGAMGCreateGraph_GEO;
739 pc_gamg->ops->coarsen = PCGAMGCoarsen_GEO;
740 pc_gamg->ops->prolongator = PCGAMGProlongator_GEO;
741 pc_gamg->ops->optprolongator = NULL;
742 pc_gamg->ops->createdefaultdata = PCSetData_GEO;
743
744 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_GEO));
745 PetscFunctionReturn(PETSC_SUCCESS);
746 }
747