Actual source code: geo.c

      1: /*
      2:  GAMG geometric-algebraic multigrid 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 (int)(((GAMGNode *)a)->degree - ((GAMGNode *)b)->degree);
     25: }

     27: // PetscClangLinter pragma disable: -fdoc-sowing-chars
     28: /*
     29:    PCSetCoordinates_GEO

     31:    Input Parameter:
     32:    .  pc - the preconditioner context
     33: */
     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;

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

     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);

     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 */

     54:   arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;

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

     78: // PetscClangLinter pragma disable: -fdoc-sowing-chars
     79: /*
     80:    PCSetData_GEO

     82:   Input Parameter:
     83:    . pc -
     84: */
     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: }

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

    109: // PetscClangLinter pragma disable: -fdoc-sowing-chars
    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:   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));

    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);

    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;

    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;

    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];

    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);

    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);

    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]);

    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;

    268:     PetscCall(PetscMalloc2(nselected_2, &node_tri, nselected_2, &nTri));

    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;

    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));

    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];

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

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

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

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

    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 */

    396:     PetscCall(ISRestoreIndices(selected_2, &selected_idx_2));
    397:     PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
    398:     PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));

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

    415: // PetscClangLinter pragma disable: -fdoc-sowing-chars
    416: /*
    417:    getGIDsOnSquareGraph - square graph, get

    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: */
    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;

    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 */

    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;

    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;

    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));

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

    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];

    524:   PetscFunctionBegin;
    525:   PetscCall(MatCreateGraph(Amat, PETSC_TRUE, PETSC_TRUE, vfilter, 0, NULL, a_Gmat));
    526:   PetscFunctionReturn(PETSC_SUCCESS);
    527: }

    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;

    539:   PetscFunctionBegin;
    540:   PetscCall(PetscObjectGetComm((PetscObject)a_pc, &comm));

    542:   PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
    543:   nloc = (Iend - Istart);

    545:   /* create random permutation with sort for geo-mg */
    546:   PetscCall(PetscMalloc1(nloc, &gnodes));
    547:   PetscCall(PetscMalloc1(nloc, &permute));

    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;

    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));

    587:   /* get MIS aggs */

    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));

    598:   PetscCall(ISDestroy(&perm));
    599:   PetscFunctionReturn(PETSC_SUCCESS);
    600: }

    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;

    615:   PetscFunctionBegin;
    616:   PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));

    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);

    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 */

    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));

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

    662:   {
    663:     PetscReal *coords;
    664:     PetscInt   data_stride;
    665:     PetscInt  *crsGID = NULL;
    666:     Mat        Gmat2;

    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;

    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));

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

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

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

    715:   *a_P_out = Prol; /* out */
    716:   PetscCall(PetscFree(clid_flid));
    717:   PetscFunctionReturn(PETSC_SUCCESS);
    718: }

    720: static PetscErrorCode PCDestroy_GAMG_GEO(PC pc)
    721: {
    722:   PetscFunctionBegin;
    723:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
    724:   PetscFunctionReturn(PETSC_SUCCESS);
    725: }

    727: PetscErrorCode PCCreateGAMG_GEO(PC pc)
    728: {
    729:   PC_MG   *mg      = (PC_MG *)pc->data;
    730:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    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 */

    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;

    744:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_GEO));
    745:   PetscFunctionReturn(PETSC_SUCCESS);
    746: }