Actual source code: geo.c

petsc-master 2019-06-26
Report Typos and Errors
  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: #define REAL PetscReal
  9: #include <triangle.h>
 10: #endif

 12:  #include <petscblaslapack.h>

 14: /* Private context for the GAMG preconditioner */
 15: typedef struct {
 16:   PetscInt lid;            /* local vertex index */
 17:   PetscInt degree;         /* vertex degree */
 18: } GAMGNode;

 20: PETSC_STATIC_INLINE int petsc_geo_mg_compare(const void *a, const void *b)
 21: {
 22:   return (((GAMGNode*)a)->degree - ((GAMGNode*)b)->degree);
 23: }

 25: /* -------------------------------------------------------------------------- */
 26: /*
 27:    PCSetCoordinates_GEO

 29:    Input Parameter:
 30:    .  pc - the preconditioner context
 31: */
 32: PetscErrorCode PCSetCoordinates_GEO(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
 33: {
 34:   PC_MG          *mg      = (PC_MG*)pc->data;
 35:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
 37:   PetscInt       arrsz,bs,my0,kk,ii,nloc,Iend,aloc;
 38:   Mat            Amat = pc->pmat;

 42:   MatGetBlockSize(Amat, &bs);
 43:   MatGetOwnershipRange(Amat, &my0, &Iend);
 44:   aloc = (Iend-my0);
 45:   nloc = (Iend-my0)/bs;

 47:   if (nloc!=a_nloc && aloc!=a_nloc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Number of local blocks %D must be %D or %D.",a_nloc,nloc,aloc);

 49:   pc_gamg->data_cell_rows = 1;
 50:   if (!coords && nloc > 0) SETERRQ(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:     PetscFree(pc_gamg->data);
 58:     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++) {
 66:         pc_gamg->data[ii*nloc + kk] =  coords[kk*ndm + ii];
 67:       }
 68:     }
 69:   } else { /* assumes the coordinates are blocked */
 70:     for (kk = 0; kk < nloc; kk++) {
 71:       for (ii = 0; ii < ndm; ii++) {
 72:         pc_gamg->data[ii*nloc + kk] =  coords[bs*kk*ndm + ii];
 73:       }
 74:     }
 75:   }
 76:   if (pc_gamg->data[arrsz] != -99.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->data[arrsz %D] %g != -99.",arrsz,pc_gamg->data[arrsz]);
 77:   pc_gamg->data_sz = arrsz;
 78:   return(0);
 79: }

 81: /* -------------------------------------------------------------------------- */
 82: /*
 83:    PCSetData_GEO

 85:   Input Parameter:
 86:    . pc -
 87: */
 88: PetscErrorCode PCSetData_GEO(PC pc, Mat m)
 89: {
 91:   SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"GEO MG needs coordinates");
 92: }

 94: /* -------------------------------------------------------------------------- */
 95: /*
 96:    PCSetFromOptions_GEO

 98:   Input Parameter:
 99:    . pc -
100: */
101: PetscErrorCode PCSetFromOptions_GEO(PetscOptionItems *PetscOptionsObject,PC pc)
102: {

106:   PetscOptionsHead(PetscOptionsObject,"GAMG-GEO options");
107:   {
108:     /* -pc_gamg_sa_nsmooths */
109:     /* pc_gamg_sa->smooths = 0; */
110:     /* PetscOptionsInt("-pc_gamg_agg_nsmooths", */
111:     /*                        "smoothing steps for smoothed aggregation, usually 1 (0)", */
112:     /*                        "PCGAMGSetNSmooths_AGG", */
113:     /*                        pc_gamg_sa->smooths, */
114:     /*                        &pc_gamg_sa->smooths, */
115:     /*                        &flag);  */
116:     /*  */
117:   }
118:   PetscOptionsTail();
119:   return(0);
120: }

122: /* -------------------------------------------------------------------------- */
123: /*
124:  triangulateAndFormProl

126:    Input Parameter:
127:    . selected_2 - list of selected local ID, includes selected ghosts
128:    . data_stride -
129:    . coords[2*data_stride] - column vector of local coordinates w/ ghosts
130:    . nselected_1 - selected IDs that go with base (1) graph includes selected ghosts
131:    . clid_lid_1[nselected_1] - lids of selected (c) nodes   ???????????
132:    . agg_lists_1 - list of aggregates selected_1 vertices of aggregate unselected vertices
133:    . crsGID[selected.size()] - global index for prolongation operator
134:    . bs - block size
135:   Output Parameter:
136:    . a_Prol - prolongation operator
137:    . a_worst_best - measure of worst missed fine vertex, 0 is no misses
138: */
139: static PetscErrorCode triangulateAndFormProl(IS selected_2,PetscInt data_stride,PetscReal coords[],PetscInt nselected_1,const PetscInt clid_lid_1[],const PetscCoarsenData *agg_lists_1,
140:                                              const PetscInt crsGID[],PetscInt bs,Mat a_Prol,PetscReal *a_worst_best)
141: {
142: #if defined(PETSC_HAVE_TRIANGLE)
143:   PetscErrorCode       ierr;
144:   PetscInt             jj,tid,tt,idx,nselected_2;
145:   struct triangulateio in,mid;
146:   const PetscInt       *selected_idx_2;
147:   PetscMPIInt          rank;
148:   PetscInt             Istart,Iend,nFineLoc,myFine0;
149:   int                  kk,nPlotPts,sid;
150:   MPI_Comm             comm;
151:   PetscReal            tm;

154:   PetscObjectGetComm((PetscObject)a_Prol,&comm);
155:   MPI_Comm_rank(comm,&rank);
156:   ISGetSize(selected_2, &nselected_2);
157:   if (nselected_2 == 1 || nselected_2 == 2) { /* 0 happens on idle processors */
158:     *a_worst_best = 100.0; /* this will cause a stop, but not globalized (should not happen) */
159:   } else *a_worst_best = 0.0;
160:   MPIU_Allreduce(a_worst_best, &tm, 1, MPIU_REAL, MPIU_MAX, comm);
161:   if (tm > 0.0) {
162:     *a_worst_best = 100.0;
163:     return(0);
164:   }
165:   MatGetOwnershipRange(a_Prol, &Istart, &Iend);
166:   nFineLoc = (Iend-Istart)/bs; myFine0 = Istart/bs;
167:   nPlotPts = nFineLoc; /* locals */
168:   /* traingle */
169:   /* Define input points - in*/
170:   in.numberofpoints          = nselected_2;
171:   in.numberofpointattributes = 0;
172:   /* get nselected points */
173:   PetscMalloc1(2*nselected_2, &in.pointlist);
174:   ISGetIndices(selected_2, &selected_idx_2);

176:   for (kk=0,sid=0; kk<nselected_2; kk++,sid += 2) {
177:     PetscInt lid = selected_idx_2[kk];
178:     in.pointlist[sid]   = coords[lid];
179:     in.pointlist[sid+1] = coords[data_stride + lid];
180:     if (lid>=nFineLoc) nPlotPts++;
181:   }
182:   if (sid != 2*nselected_2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %D != 2*nselected_2 %D",sid,nselected_2);

184:   in.numberofsegments      = 0;
185:   in.numberofedges         = 0;
186:   in.numberofholes         = 0;
187:   in.numberofregions       = 0;
188:   in.trianglelist          = 0;
189:   in.segmentmarkerlist     = 0;
190:   in.pointattributelist    = 0;
191:   in.pointmarkerlist       = 0;
192:   in.triangleattributelist = 0;
193:   in.trianglearealist      = 0;
194:   in.segmentlist           = 0;
195:   in.holelist              = 0;
196:   in.regionlist            = 0;
197:   in.edgelist              = 0;
198:   in.edgemarkerlist        = 0;
199:   in.normlist              = 0;

201:   /* triangulate */
202:   mid.pointlist = 0;            /* Not needed if -N switch used. */
203:   /* Not needed if -N switch used or number of point attributes is zero: */
204:   mid.pointattributelist = 0;
205:   mid.pointmarkerlist    = 0; /* Not needed if -N or -B switch used. */
206:   mid.trianglelist       = 0;    /* Not needed if -E switch used. */
207:   /* Not needed if -E switch used or number of triangle attributes is zero: */
208:   mid.triangleattributelist = 0;
209:   mid.neighborlist          = 0; /* Needed only if -n switch used. */
210:   /* Needed only if segments are output (-p or -c) and -P not used: */
211:   mid.segmentlist = 0;
212:   /* Needed only if segments are output (-p or -c) and -P and -B not used: */
213:   mid.segmentmarkerlist = 0;
214:   mid.edgelist          = 0;    /* Needed only if -e switch used. */
215:   mid.edgemarkerlist    = 0; /* Needed if -e used and -B not used. */
216:   mid.numberoftriangles = 0;

218:   /* Triangulate the points.  Switches are chosen to read and write a  */
219:   /*   PSLG (p), preserve the convex hull (c), number everything from  */
220:   /*   zero (z), assign a regional attribute to each element (A), and  */
221:   /*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
222:   /*   neighbor list (n).                                            */
223:   if (nselected_2 != 0) { /* inactive processor */
224:     char args[] = "npczQ"; /* c is needed ? */
225:     triangulate(args, &in, &mid, (struct triangulateio*) NULL);
226:     /* output .poly files for 'showme' */
227:     if (!PETSC_TRUE) {
228:       static int level = 1;
229:       FILE       *file; char fname[32];

231:       sprintf(fname,"C%d_%d.poly",level,rank); file = fopen(fname, "w");
232:       /*First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>*/
233:       fprintf(file, "%d  %d  %d  %d\n",in.numberofpoints,2,0,0);
234:       /*Following lines: <vertex #> <x> <y> */
235:       for (kk=0,sid=0; kk<in.numberofpoints; kk++,sid += 2) {
236:         fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]);
237:       }
238:       /*One line: <# of segments> <# of boundary markers (0 or 1)> */
239:       fprintf(file, "%d  %d\n",0,0);
240:       /*Following lines: <segment #> <endpoint> <endpoint> [boundary marker] */
241:       /* One line: <# of holes> */
242:       fprintf(file, "%d\n",0);
243:       /* Following lines: <hole #> <x> <y> */
244:       /* Optional line: <# of regional attributes and/or area constraints> */
245:       /* Optional following lines: <region #> <x> <y> <attribute> <maximum area> */
246:       fclose(file);

248:       /* elems */
249:       sprintf(fname,"C%d_%d.ele",level,rank); file = fopen(fname, "w");
250:       /* First line: <# of triangles> <nodes per triangle> <# of attributes> */
251:       fprintf(file, "%d %d %d\n",mid.numberoftriangles,3,0);
252:       /* Remaining lines: <triangle #> <node> <node> <node> ... [attributes] */
253:       for (kk=0,sid=0; kk<mid.numberoftriangles; kk++,sid += 3) {
254:         fprintf(file, "%d %d %d %d\n",kk,mid.trianglelist[sid],mid.trianglelist[sid+1],mid.trianglelist[sid+2]);
255:       }
256:       fclose(file);

258:       sprintf(fname,"C%d_%d.node",level,rank); file = fopen(fname, "w");
259:       /* First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> */
260:       /* fprintf(file, "%d  %d  %d  %d\n",in.numberofpoints,2,0,0); */
261:       fprintf(file, "%d  %d  %d  %d\n",nPlotPts,2,0,0);
262:       /*Following lines: <vertex #> <x> <y> */
263:       for (kk=0,sid=0; kk<in.numberofpoints; kk++,sid+=2) {
264:         fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]);
265:       }

267:       sid /= 2;
268:       for (jj=0; jj<nFineLoc; jj++) {
269:         PetscBool sel = PETSC_TRUE;
270:         for (kk=0; kk<nselected_2 && sel; kk++) {
271:           PetscInt lid = selected_idx_2[kk];
272:           if (lid == jj) sel = PETSC_FALSE;
273:         }
274:         if (sel) fprintf(file, "%d %e %e\n",sid++,coords[jj],coords[data_stride + jj]);
275:       }
276:       fclose(file);
277:       if (sid != nPlotPts) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %D != nPlotPts %D",sid,nPlotPts);
278:       level++;
279:     }
280:   }
281: #if defined PETSC_GAMG_USE_LOG
282:   PetscLogEventBegin(petsc_gamg_setup_events[FIND_V],0,0,0,0);
283: #endif
284:   { /* form P - setup some maps */
285:     PetscInt clid,mm,*nTri,*node_tri;

287:     PetscMalloc2(nselected_2, &node_tri,nselected_2, &nTri);

289:     /* need list of triangles on node */
290:     for (kk=0; kk<nselected_2; kk++) nTri[kk] = 0;
291:     for (tid=0,kk=0; tid<mid.numberoftriangles; tid++) {
292:       for (jj=0; jj<3; jj++) {
293:         PetscInt cid = mid.trianglelist[kk++];
294:         if (nTri[cid] == 0) node_tri[cid] = tid;
295:         nTri[cid]++;
296:       }
297:     }
298: #define EPS 1.e-12
299:     /* find points and set prolongation */
300:     for (mm = clid = 0; mm < nFineLoc; mm++) {
301:       PetscBool ise;
302:       PetscCDEmptyAt(agg_lists_1,mm,&ise);
303:       if (!ise) {
304:         const PetscInt lid = mm;
305:         PetscScalar    AA[3][3];
306:         PetscBLASInt   N=3,NRHS=1,LDA=3,IPIV[3],LDB=3,INFO;
307:         PetscCDIntNd   *pos;

309:         PetscCDGetHeadPos(agg_lists_1,lid,&pos);
310:         while (pos) {
311:           PetscInt flid;
312:           PetscCDIntNdGetID(pos, &flid);
313:           PetscCDGetNextPos(agg_lists_1,lid,&pos);

315:           if (flid < nFineLoc) {  /* could be a ghost */
316:             PetscInt       bestTID = -1; PetscReal best_alpha = 1.e10;
317:             const PetscInt fgid    = flid + myFine0;
318:             /* compute shape function for gid */
319:             const PetscReal fcoord[3] = {coords[flid],coords[data_stride+flid],1.0};
320:             PetscBool       haveit    =PETSC_FALSE; PetscScalar alpha[3]; PetscInt clids[3];

322:             /* look for it */
323:             for (tid = node_tri[clid], jj=0;
324:                  jj < 5 && !haveit && tid != -1;
325:                  jj++) {
326:               for (tt=0; tt<3; tt++) {
327:                 PetscInt cid2 = mid.trianglelist[3*tid + tt];
328:                 PetscInt lid2 = selected_idx_2[cid2];
329:                 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
330:                 clids[tt] = cid2; /* store for interp */
331:               }

333:               for (tt=0; tt<3; tt++) alpha[tt] = (PetscScalar)fcoord[tt];

335:               /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
336:               PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
337:               {
338:                 PetscBool have=PETSC_TRUE;  PetscReal lowest=1.e10;
339:                 for (tt = 0, idx = 0; tt < 3; tt++) {
340:                   if (PetscRealPart(alpha[tt]) > (1.0+EPS) || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE;
341:                   if (PetscRealPart(alpha[tt]) < lowest) {
342:                     lowest = PetscRealPart(alpha[tt]);
343:                     idx    = tt;
344:                   }
345:                 }
346:                 haveit = have;
347:               }
348:               tid = mid.neighborlist[3*tid + idx];
349:             }

351:             if (!haveit) {
352:               /* brute force */
353:               for (tid=0; tid<mid.numberoftriangles && !haveit; tid++) {
354:                 for (tt=0; tt<3; tt++) {
355:                   PetscInt cid2 = mid.trianglelist[3*tid + tt];
356:                   PetscInt lid2 = selected_idx_2[cid2];
357:                   AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
358:                   clids[tt] = cid2; /* store for interp */
359:                 }
360:                 for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt];
361:                 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
362:                 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
363:                 {
364:                   PetscBool have=PETSC_TRUE;  PetscReal worst=0.0, v;
365:                   for (tt=0; tt<3 && have; tt++) {
366:                     if (PetscRealPart(alpha[tt]) > 1.0+EPS || PetscRealPart(alpha[tt]) < -EPS) have=PETSC_FALSE;
367:                     if ((v=PetscAbs(PetscRealPart(alpha[tt])-0.5)) > worst) worst = v;
368:                   }
369:                   if (worst < best_alpha) {
370:                     best_alpha = worst; bestTID = tid;
371:                   }
372:                   haveit = have;
373:                 }
374:               }
375:             }
376:             if (!haveit) {
377:               if (best_alpha > *a_worst_best) *a_worst_best = best_alpha;
378:               /* use best one */
379:               for (tt=0; tt<3; tt++) {
380:                 PetscInt cid2 = mid.trianglelist[3*bestTID + tt];
381:                 PetscInt lid2 = selected_idx_2[cid2];
382:                 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
383:                 clids[tt] = cid2; /* store for interp */
384:               }
385:               for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt];
386:               /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
387:               PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
388:             }

390:             /* put in row of P */
391:             for (idx=0; idx<3; idx++) {
392:               PetscScalar shp = alpha[idx];
393:               if (PetscAbs(PetscRealPart(shp)) > 1.e-6) {
394:                 PetscInt cgid = crsGID[clids[idx]];
395:                 PetscInt jj   = cgid*bs, ii = fgid*bs; /* need to gloalize */
396:                 for (tt=0; tt < bs; tt++, ii++, jj++) {
397:                   MatSetValues(a_Prol,1,&ii,1,&jj,&shp,INSERT_VALUES);
398:                 }
399:               }
400:             }
401:           }
402:         } /* aggregates iterations */
403:         clid++;
404:       } /* a coarse agg */
405:     } /* for all fine nodes */

407:     ISRestoreIndices(selected_2, &selected_idx_2);
408:     MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
409:     MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);

411:     PetscFree2(node_tri,nTri);
412:   }
413: #if defined PETSC_GAMG_USE_LOG
414:   PetscLogEventEnd(petsc_gamg_setup_events[FIND_V],0,0,0,0);
415: #endif
416:   free(mid.trianglelist);
417:   free(mid.neighborlist);
418:   PetscFree(in.pointlist);
419:   return(0);
420: #else
421:   SETERRQ(PetscObjectComm((PetscObject)a_Prol),PETSC_ERR_PLIB,"configure with TRIANGLE to use geometric MG");
422: #endif
423: }
424: /* -------------------------------------------------------------------------- */
425: /*
426:    getGIDsOnSquareGraph - square graph, get

428:    Input Parameter:
429:    . nselected_1 - selected local indices (includes ghosts in input Gmat1)
430:    . clid_lid_1 - [nselected_1] lids of selected nodes
431:    . Gmat1 - graph that goes with 'selected_1'
432:    Output Parameter:
433:    . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2)
434:    . a_Gmat_2 - graph that is squared of 'Gmat_1'
435:    . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes
436: */
437: static PetscErrorCode getGIDsOnSquareGraph(PetscInt nselected_1,const PetscInt clid_lid_1[],const Mat Gmat1,IS *a_selected_2,Mat *a_Gmat_2,PetscInt **a_crsGID)
438: {
440:   PetscMPIInt    size;
441:   PetscInt       *crsGID, kk,my0,Iend,nloc;
442:   MPI_Comm       comm;

445:   PetscObjectGetComm((PetscObject)Gmat1,&comm);
446:   MPI_Comm_size(comm,&size);
447:   MatGetOwnershipRange(Gmat1,&my0,&Iend); /* AIJ */
448:   nloc = Iend - my0; /* this does not change */

450:   if (size == 1) { /* not much to do in serial */
451:     PetscMalloc1(nselected_1, &crsGID);
452:     for (kk=0; kk<nselected_1; kk++) crsGID[kk] = kk;
453:     *a_Gmat_2 = 0;
454:     ISCreateGeneral(PETSC_COMM_SELF,nselected_1,clid_lid_1,PETSC_COPY_VALUES,a_selected_2);
455:   } else {
456:     PetscInt    idx,num_fine_ghosts,num_crs_ghost,myCrs0;
457:     Mat_MPIAIJ  *mpimat2;
458:     Mat         Gmat2;
459:     Vec         locState;
460:     PetscScalar *cpcol_state;

462:     /* scan my coarse zero gid, set 'lid_state' with coarse GID */
463:     kk = nselected_1;
464:     MPI_Scan(&kk, &myCrs0, 1, MPIU_INT, MPI_SUM, comm);
465:     myCrs0 -= nselected_1;

467:     if (a_Gmat_2) { /* output */
468:       /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */
469:       MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);
470:       *a_Gmat_2 = Gmat2; /* output */
471:     } else Gmat2 = Gmat1;  /* use local to get crsGIDs at least */
472:     /* get coarse grid GIDS for selected (locals and ghosts) */
473:     mpimat2 = (Mat_MPIAIJ*)Gmat2->data;
474:     MatCreateVecs(Gmat2, &locState, 0);
475:     VecSet(locState, (PetscScalar)(PetscReal)(-1)); /* set with UNKNOWN state */
476:     for (kk=0; kk<nselected_1; kk++) {
477:       PetscInt    fgid = clid_lid_1[kk] + my0;
478:       PetscScalar v    = (PetscScalar)(kk+myCrs0);
479:       VecSetValues(locState, 1, &fgid, &v, INSERT_VALUES); /* set with PID */
480:     }
481:     VecAssemblyBegin(locState);
482:     VecAssemblyEnd(locState);
483:     VecScatterBegin(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);
484:       VecScatterEnd(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);
485:     VecGetLocalSize(mpimat2->lvec, &num_fine_ghosts);
486:     VecGetArray(mpimat2->lvec, &cpcol_state);
487:     for (kk=0,num_crs_ghost=0; kk<num_fine_ghosts; kk++) {
488:       if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) num_crs_ghost++;
489:     }
490:     PetscMalloc1(nselected_1+num_crs_ghost, &crsGID); /* output */
491:     {
492:       PetscInt *selected_set;
493:       PetscMalloc1(nselected_1+num_crs_ghost, &selected_set);
494:       /* do ghost of 'crsGID' */
495:       for (kk=0,idx=nselected_1; kk<num_fine_ghosts; kk++) {
496:         if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
497:           PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
498:           selected_set[idx] = nloc + kk;
499:           crsGID[idx++]     = cgid;
500:         }
501:       }
502:       if (idx != (nselected_1+num_crs_ghost)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"idx %D != (nselected_1 %D + num_crs_ghost %D)",idx,nselected_1,num_crs_ghost);
503:       VecRestoreArray(mpimat2->lvec, &cpcol_state);
504:       /* do locals in 'crsGID' */
505:       VecGetArray(locState, &cpcol_state);
506:       for (kk=0,idx=0; kk<nloc; kk++) {
507:         if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
508:           PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
509:           selected_set[idx] = kk;
510:           crsGID[idx++]     = cgid;
511:         }
512:       }
513:       if (idx != nselected_1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"idx %D != nselected_1 %D",idx,nselected_1);
514:       VecRestoreArray(locState, &cpcol_state);

516:       if (a_selected_2 != 0) { /* output */
517:         ISCreateGeneral(PETSC_COMM_SELF,(nselected_1+num_crs_ghost),selected_set,PETSC_OWN_POINTER,a_selected_2);
518:       } else {
519:         PetscFree(selected_set);
520:       }
521:     }
522:     VecDestroy(&locState);
523:   }
524:   *a_crsGID = crsGID; /* output */
525:   return(0);
526: }

528: /* -------------------------------------------------------------------------- */
529: /*
530:    PCGAMGGraph_GEO

532:   Input Parameter:
533:    . pc - this
534:    . Amat - matrix on this fine level
535:   Output Parameter:
536:    . a_Gmat
537: */
538: PetscErrorCode PCGAMGGraph_GEO(PC pc,Mat Amat,Mat *a_Gmat)
539: {
540:   PetscErrorCode  ierr;
541:   PC_MG           *mg      = (PC_MG*)pc->data;
542:   PC_GAMG         *pc_gamg = (PC_GAMG*)mg->innerctx;
543:   const PetscReal vfilter  = pc_gamg->threshold[0];
544:   MPI_Comm        comm;
545:   Mat             Gmat;
546:   PetscBool       set,flg,symm;

549:   PetscObjectGetComm((PetscObject)Amat,&comm);
550:   PetscLogEventBegin(PC_GAMGGraph_GEO,0,0,0,0);

552:   MatIsSymmetricKnown(Amat, &set, &flg);
553:   symm = (PetscBool)!(set && flg);

555:   PCGAMGCreateGraph(Amat, &Gmat);
556:   PCGAMGFilterGraph(&Gmat, vfilter, symm);

558:   *a_Gmat = Gmat;
559:   PetscLogEventEnd(PC_GAMGGraph_GEO,0,0,0,0);
560:   return(0);
561: }

563: /* -------------------------------------------------------------------------- */
564: /*
565:    PCGAMGCoarsen_GEO

567:   Input Parameter:
568:    . a_pc - this
569:    . a_Gmat - graph
570:   Output Parameter:
571:    . a_llist_parent - linked list from selected indices for data locality only
572: */
573: PetscErrorCode PCGAMGCoarsen_GEO(PC a_pc,Mat *a_Gmat,PetscCoarsenData **a_llist_parent)
574: {
576:   PetscInt       Istart,Iend,nloc,kk,Ii,ncols;
577:   IS             perm;
578:   GAMGNode       *gnodes;
579:   PetscInt       *permute;
580:   Mat            Gmat  = *a_Gmat;
581:   MPI_Comm       comm;
582:   MatCoarsen     crs;

585:   PetscObjectGetComm((PetscObject)a_pc,&comm);
586:   PetscLogEventBegin(PC_GAMGCoarsen_GEO,0,0,0,0);
587:   MatGetOwnershipRange(Gmat, &Istart, &Iend);
588:   nloc = (Iend-Istart);

590:   /* create random permutation with sort for geo-mg */
591:   PetscMalloc1(nloc, &gnodes);
592:   PetscMalloc1(nloc, &permute);

594:   for (Ii=Istart; Ii<Iend; Ii++) { /* locals only? */
595:     MatGetRow(Gmat,Ii,&ncols,0,0);
596:     {
597:       PetscInt lid = Ii - Istart;
598:       gnodes[lid].lid    = lid;
599:       gnodes[lid].degree = ncols;
600:     }
601:     MatRestoreRow(Gmat,Ii,&ncols,0,0);
602:   }
603:   if (PETSC_TRUE) {
604:     PetscRandom  rand;
605:     PetscBool    *bIndexSet;
606:     PetscReal    rr;
607:     PetscInt     iSwapIndex;

609:     PetscRandomCreate(comm,&rand);
610:     PetscCalloc1(nloc, &bIndexSet);
611:     for (Ii = 0; Ii < nloc; Ii++) {
612:       PetscRandomGetValueReal(rand,&rr);
613:       iSwapIndex = (PetscInt) (rr*nloc);
614:       if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
615:         GAMGNode iTemp = gnodes[iSwapIndex];
616:         gnodes[iSwapIndex]    = gnodes[Ii];
617:         gnodes[Ii]            = iTemp;
618:         bIndexSet[Ii]         = PETSC_TRUE;
619:         bIndexSet[iSwapIndex] = PETSC_TRUE;
620:       }
621:     }
622:     PetscRandomDestroy(&rand);
623:     PetscFree(bIndexSet);
624:   }
625:   /* only sort locals */
626:   qsort(gnodes, nloc, sizeof(GAMGNode), petsc_geo_mg_compare);
627:   /* create IS of permutation */
628:   for (kk=0; kk<nloc; kk++) permute[kk] = gnodes[kk].lid; /* locals only */
629:   PetscFree(gnodes);
630:   ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_OWN_POINTER, &perm);

632:   /* get MIS aggs */

634:   MatCoarsenCreate(comm, &crs);
635:   MatCoarsenSetType(crs, MATCOARSENMIS);
636:   MatCoarsenSetGreedyOrdering(crs, perm);
637:   MatCoarsenSetAdjacency(crs, Gmat);
638:   MatCoarsenSetStrictAggs(crs, PETSC_FALSE);
639:   MatCoarsenApply(crs);
640:   MatCoarsenGetData(crs, a_llist_parent);
641:   MatCoarsenDestroy(&crs);

643:   ISDestroy(&perm);
644:   PetscLogEventEnd(PC_GAMGCoarsen_GEO,0,0,0,0);
645:   return(0);
646: }

648: /* -------------------------------------------------------------------------- */
649: /*
650:  PCGAMGProlongator_GEO

652:  Input Parameter:
653:  . pc - this
654:  . Amat - matrix on this fine level
655:  . Graph - used to get ghost data for nodes in
656:  . selected_1 - [nselected]
657:  . agg_lists - [nselected]
658:  Output Parameter:
659:  . a_P_out - prolongation operator to the next level
660:  */
661: PetscErrorCode PCGAMGProlongator_GEO(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
662: {
663:   PC_MG          *mg      = (PC_MG*)pc->data;
664:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
665:   const PetscInt dim      = pc_gamg->data_cell_cols, data_cols = pc_gamg->data_cell_cols;
667:   PetscInt       Istart,Iend,nloc,my0,jj,kk,ncols,nLocalSelected,bs,*clid_flid;
668:   Mat            Prol;
669:   PetscMPIInt    rank, size;
670:   MPI_Comm       comm;
671:   IS             selected_2,selected_1;
672:   const PetscInt *selected_idx;
673:   MatType        mtype;

676:   PetscObjectGetComm((PetscObject)Amat,&comm);
677:   PetscLogEventBegin(PC_GAMGProlongator_GEO,0,0,0,0);
678:   MPI_Comm_rank(comm,&rank);
679:   MPI_Comm_size(comm,&size);
680:   MatGetOwnershipRange(Amat, &Istart, &Iend);
681:   MatGetBlockSize(Amat, &bs);
682:   nloc = (Iend-Istart)/bs; my0 = Istart/bs;
683:   if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) % bs %D",Iend,Istart,bs);

685:   /* get 'nLocalSelected' */
686:   PetscCDGetMIS(agg_lists, &selected_1);
687:   ISGetSize(selected_1, &jj);
688:   PetscMalloc1(jj, &clid_flid);
689:   ISGetIndices(selected_1, &selected_idx);
690:   for (kk=0,nLocalSelected=0; kk<jj; kk++) {
691:     PetscInt lid = selected_idx[kk];
692:     if (lid<nloc) {
693:       MatGetRow(Gmat,lid+my0,&ncols,0,0);
694:       if (ncols>1) clid_flid[nLocalSelected++] = lid; /* fiter out singletons */
695:       MatRestoreRow(Gmat,lid+my0,&ncols,0,0);
696:     }
697:   }
698:   ISRestoreIndices(selected_1, &selected_idx);
699:   ISDestroy(&selected_1); /* this is selected_1 in serial */

701:   /* create prolongator  matrix */
702:   MatGetType(Amat,&mtype);
703:   MatCreate(comm, &Prol);
704:   MatSetSizes(Prol,nloc*bs,nLocalSelected*bs,PETSC_DETERMINE,PETSC_DETERMINE);
705:   MatSetBlockSizes(Prol, bs, bs);
706:   MatSetType(Prol, mtype);
707:   MatSeqAIJSetPreallocation(Prol,3*data_cols,NULL);
708:   MatMPIAIJSetPreallocation(Prol,3*data_cols,NULL,3*data_cols,NULL);

710:   /* can get all points "removed" - but not on geomg */
711:    MatGetSize(Prol, &kk, &jj);
712:   if (!jj) {
713:     PetscInfo(pc,"ERROE: no selected points on coarse grid\n");
714:     PetscFree(clid_flid);
715:     MatDestroy(&Prol);
716:     *a_P_out = NULL;  /* out */
717:     return(0);
718:   }

720:   {
721:     PetscReal *coords;
722:     PetscInt  data_stride;
723:     PetscInt  *crsGID = NULL;
724:     Mat       Gmat2;

726:     if (dim != data_cols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"dim %D != data_cols %D",dim,data_cols);
727:     /* grow ghost data for better coarse grid cover of fine grid */
728: #if defined PETSC_GAMG_USE_LOG
729:     PetscLogEventBegin(petsc_gamg_setup_events[SET5],0,0,0,0);
730: #endif
731:     /* messy method, squares graph and gets some data */
732:     getGIDsOnSquareGraph(nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID);
733:     /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */
734: #if defined PETSC_GAMG_USE_LOG
735:     PetscLogEventEnd(petsc_gamg_setup_events[SET5],0,0,0,0);
736: #endif
737:     /* create global vector of coorindates in 'coords' */
738:     if (size > 1) {
739:       PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords);
740:     } else {
741:       coords      = (PetscReal*)pc_gamg->data;
742:       data_stride = pc_gamg->data_sz/pc_gamg->data_cell_cols;
743:     }
744:     MatDestroy(&Gmat2);

746:     /* triangulate */
747:     if (dim == 2) {
748:       PetscReal metric,tm;
749: #if defined PETSC_GAMG_USE_LOG
750:       PetscLogEventBegin(petsc_gamg_setup_events[SET6],0,0,0,0);
751: #endif
752:       triangulateAndFormProl(selected_2, data_stride, coords,nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric);
753: #if defined PETSC_GAMG_USE_LOG
754:       PetscLogEventEnd(petsc_gamg_setup_events[SET6],0,0,0,0);
755: #endif
756:       PetscFree(crsGID);

758:       /* clean up and create coordinates for coarse grid (output) */
759:       if (size > 1) PetscFree(coords);

761:       MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm);
762:       if (tm > 1.) { /* needs to be globalized - should not happen */
763:         PetscInfo1(pc," failed metric for coarse grid %e\n",(double)tm);
764:         MatDestroy(&Prol);
765:       } else if (metric > .0) {
766:         PetscInfo1(pc,"worst metric for coarse grid = %e\n",(double)metric);
767:       }
768:     } else SETERRQ(comm,PETSC_ERR_PLIB,"3D not implemented for 'geo' AMG");
769:     { /* create next coords - output */
770:       PetscReal *crs_crds;
771:       PetscMalloc1(dim*nLocalSelected, &crs_crds);
772:       for (kk=0; kk<nLocalSelected; kk++) { /* grab local select nodes to promote - output */
773:         PetscInt lid = clid_flid[kk];
774:         for (jj=0; jj<dim; jj++) crs_crds[jj*nLocalSelected + kk] = pc_gamg->data[jj*nloc + lid];
775:       }

777:       PetscFree(pc_gamg->data);
778:       pc_gamg->data    = crs_crds; /* out */
779:       pc_gamg->data_sz = dim*nLocalSelected;
780:     }
781:     ISDestroy(&selected_2);
782:   }

784:   *a_P_out = Prol;  /* out */
785:   PetscFree(clid_flid);
786:   PetscLogEventEnd(PC_GAMGProlongator_GEO,0,0,0,0);
787:   return(0);
788: }

790: static PetscErrorCode PCDestroy_GAMG_GEO(PC pc)
791: {

795:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
796:   return(0);
797: }

799: /* -------------------------------------------------------------------------- */
800: /*
801:  PCCreateGAMG_GEO

803:   Input Parameter:
804:    . pc -
805: */
806: PetscErrorCode  PCCreateGAMG_GEO(PC pc)
807: {
809:   PC_MG          *mg      = (PC_MG*)pc->data;
810:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;

813:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO;
814:   pc_gamg->ops->destroy        = PCDestroy_GAMG_GEO;
815:   /* reset does not do anything; setup not virtual */

817:   /* set internal function pointers */
818:   pc_gamg->ops->graph             = PCGAMGGraph_GEO;
819:   pc_gamg->ops->coarsen           = PCGAMGCoarsen_GEO;
820:   pc_gamg->ops->prolongator       = PCGAMGProlongator_GEO;
821:   pc_gamg->ops->optprolongator    = NULL;
822:   pc_gamg->ops->createdefaultdata = PCSetData_GEO;

824:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_GEO);
825:   return(0);
826: }