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: }