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