1 /* 2 GAMG geometric-algebric multiogrid PC - Mark Adams 2011 3 */ 4 5 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 6 7 #if defined(PETSC_HAVE_TRIANGLE) 8 #if !defined(ANSI_DECLARATORS) 9 #define ANSI_DECLARATORS 10 #endif 11 #include <triangle.h> 12 #endif 13 14 #include <petscblaslapack.h> 15 16 /* Private context for the GAMG preconditioner */ 17 typedef struct { 18 PetscInt lid; /* local vertex index */ 19 PetscInt degree; /* vertex degree */ 20 } GAMGNode; 21 22 static inline int petsc_geo_mg_compare(const void *a, const void *b) { 23 return (((GAMGNode *)a)->degree - ((GAMGNode *)b)->degree); 24 } 25 26 /* -------------------------------------------------------------------------- */ 27 /* 28 PCSetCoordinates_GEO 29 30 Input Parameter: 31 . pc - the preconditioner context 32 */ 33 PetscErrorCode PCSetCoordinates_GEO(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords) { 34 PC_MG *mg = (PC_MG *)pc->data; 35 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 36 PetscInt arrsz, bs, my0, kk, ii, nloc, Iend, aloc; 37 Mat Amat = pc->pmat; 38 39 PetscFunctionBegin; 40 PetscValidHeaderSpecific(Amat, MAT_CLASSID, 1); 41 PetscCall(MatGetBlockSize(Amat, &bs)); 42 PetscCall(MatGetOwnershipRange(Amat, &my0, &Iend)); 43 aloc = (Iend - my0); 44 nloc = (Iend - my0) / bs; 45 46 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); 47 48 pc_gamg->data_cell_rows = 1; 49 PetscCheck(coords || (nloc <= 0), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need coordinates for pc_gamg_type 'geo'."); 50 pc_gamg->data_cell_cols = ndm; /* coordinates */ 51 52 arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols; 53 54 /* create data - syntactic sugar that should be refactored at some point */ 55 if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) { 56 PetscCall(PetscFree(pc_gamg->data)); 57 PetscCall(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 } 71 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]); 72 pc_gamg->data_sz = arrsz; 73 PetscFunctionReturn(0); 74 } 75 76 /* -------------------------------------------------------------------------- */ 77 /* 78 PCSetData_GEO 79 80 Input Parameter: 81 . pc - 82 */ 83 PetscErrorCode PCSetData_GEO(PC pc, Mat m) { 84 PetscFunctionBegin; 85 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "GEO MG needs coordinates"); 86 } 87 88 /* -------------------------------------------------------------------------- */ 89 /* 90 PCSetFromOptions_GEO 91 92 Input Parameter: 93 . pc - 94 */ 95 PetscErrorCode PCSetFromOptions_GEO(PC pc, PetscOptionItems *PetscOptionsObject) { 96 PetscFunctionBegin; 97 PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-GEO options"); 98 { 99 /* -pc_gamg_sa_nsmooths */ 100 /* pc_gamg_sa->smooths = 0; */ 101 /* ierr = PetscOptionsInt("-pc_gamg_agg_nsmooths", */ 102 /* "smoothing steps for smoothed aggregation, usually 1 (0)", */ 103 /* "PCGAMGSetNSmooths_AGG", */ 104 /* pc_gamg_sa->smooths, */ 105 /* &pc_gamg_sa->smooths, */ 106 /* &flag); */ 107 } 108 PetscOptionsHeadEnd(); 109 PetscFunctionReturn(0); 110 } 111 112 /* -------------------------------------------------------------------------- */ 113 /* 114 triangulateAndFormProl 115 116 Input Parameter: 117 . selected_2 - list of selected local ID, includes selected ghosts 118 . data_stride - 119 . coords[2*data_stride] - column vector of local coordinates w/ ghosts 120 . nselected_1 - selected IDs that go with base (1) graph includes selected ghosts 121 . clid_lid_1[nselected_1] - lids of selected (c) nodes ??????????? 122 . agg_lists_1 - list of aggregates selected_1 vertices of aggregate unselected vertices 123 . crsGID[selected.size()] - global index for prolongation operator 124 . bs - block size 125 Output Parameter: 126 . a_Prol - prolongation operator 127 . a_worst_best - measure of worst missed fine vertex, 0 is no misses 128 */ 129 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) { 130 #if defined(PETSC_HAVE_TRIANGLE) 131 PetscInt jj, tid, tt, idx, nselected_2; 132 struct triangulateio in, mid; 133 const PetscInt *selected_idx_2; 134 PetscMPIInt rank; 135 PetscInt Istart, Iend, nFineLoc, myFine0; 136 int kk, nPlotPts, sid; 137 MPI_Comm comm; 138 PetscReal tm; 139 140 PetscFunctionBegin; 141 PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm)); 142 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 143 PetscCall(ISGetSize(selected_2, &nselected_2)); 144 if (nselected_2 == 1 || nselected_2 == 2) { /* 0 happens on idle processors */ 145 *a_worst_best = 100.0; /* this will cause a stop, but not globalized (should not happen) */ 146 } else *a_worst_best = 0.0; 147 PetscCall(MPIU_Allreduce(a_worst_best, &tm, 1, MPIU_REAL, MPIU_MAX, comm)); 148 if (tm > 0.0) { 149 *a_worst_best = 100.0; 150 PetscFunctionReturn(0); 151 } 152 PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend)); 153 nFineLoc = (Iend - Istart) / bs; 154 myFine0 = Istart / bs; 155 nPlotPts = nFineLoc; /* locals */ 156 /* triangle */ 157 /* Define input points - in*/ 158 in.numberofpoints = nselected_2; 159 in.numberofpointattributes = 0; 160 /* get nselected points */ 161 PetscCall(PetscMalloc1(2 * nselected_2, &in.pointlist)); 162 PetscCall(ISGetIndices(selected_2, &selected_idx_2)); 163 164 for (kk = 0, sid = 0; kk < nselected_2; kk++, sid += 2) { 165 PetscInt lid = selected_idx_2[kk]; 166 in.pointlist[sid] = coords[lid]; 167 in.pointlist[sid + 1] = coords[data_stride + lid]; 168 if (lid >= nFineLoc) nPlotPts++; 169 } 170 PetscCheck(sid == 2 * nselected_2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "sid %d != 2*nselected_2 %" PetscInt_FMT, sid, nselected_2); 171 172 in.numberofsegments = 0; 173 in.numberofedges = 0; 174 in.numberofholes = 0; 175 in.numberofregions = 0; 176 in.trianglelist = NULL; 177 in.segmentmarkerlist = NULL; 178 in.pointattributelist = NULL; 179 in.pointmarkerlist = NULL; 180 in.triangleattributelist = NULL; 181 in.trianglearealist = NULL; 182 in.segmentlist = NULL; 183 in.holelist = NULL; 184 in.regionlist = NULL; 185 in.edgelist = NULL; 186 in.edgemarkerlist = NULL; 187 in.normlist = NULL; 188 189 /* triangulate */ 190 mid.pointlist = NULL; /* Not needed if -N switch used. */ 191 /* Not needed if -N switch used or number of point attributes is zero: */ 192 mid.pointattributelist = NULL; 193 mid.pointmarkerlist = NULL; /* Not needed if -N or -B switch used. */ 194 mid.trianglelist = NULL; /* Not needed if -E switch used. */ 195 /* Not needed if -E switch used or number of triangle attributes is zero: */ 196 mid.triangleattributelist = NULL; 197 mid.neighborlist = NULL; /* Needed only if -n switch used. */ 198 /* Needed only if segments are output (-p or -c) and -P not used: */ 199 mid.segmentlist = NULL; 200 /* Needed only if segments are output (-p or -c) and -P and -B not used: */ 201 mid.segmentmarkerlist = NULL; 202 mid.edgelist = NULL; /* Needed only if -e switch used. */ 203 mid.edgemarkerlist = NULL; /* Needed if -e used and -B not used. */ 204 mid.numberoftriangles = 0; 205 206 /* Triangulate the points. Switches are chosen to read and write a */ 207 /* PSLG (p), preserve the convex hull (c), number everything from */ 208 /* zero (z), assign a regional attribute to each element (A), and */ 209 /* produce an edge list (e), a Voronoi diagram (v), and a triangle */ 210 /* neighbor list (n). */ 211 if (nselected_2 != 0) { /* inactive processor */ 212 char args[] = "npczQ"; /* c is needed ? */ 213 triangulate(args, &in, &mid, (struct triangulateio *)NULL); 214 /* output .poly files for 'showme' */ 215 if (!PETSC_TRUE) { 216 static int level = 1; 217 FILE *file; 218 char fname[32]; 219 220 sprintf(fname, "C%d_%d.poly", level, rank); 221 file = fopen(fname, "w"); 222 /*First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>*/ 223 fprintf(file, "%d %d %d %d\n", in.numberofpoints, 2, 0, 0); 224 /*Following lines: <vertex #> <x> <y> */ 225 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]); } 226 /*One line: <# of segments> <# of boundary markers (0 or 1)> */ 227 fprintf(file, "%d %d\n", 0, 0); 228 /*Following lines: <segment #> <endpoint> <endpoint> [boundary marker] */ 229 /* One line: <# of holes> */ 230 fprintf(file, "%d\n", 0); 231 /* Following lines: <hole #> <x> <y> */ 232 /* Optional line: <# of regional attributes and/or area constraints> */ 233 /* Optional following lines: <region #> <x> <y> <attribute> <maximum area> */ 234 fclose(file); 235 236 /* elems */ 237 sprintf(fname, "C%d_%d.ele", level, rank); 238 file = fopen(fname, "w"); 239 /* First line: <# of triangles> <nodes per triangle> <# of attributes> */ 240 fprintf(file, "%d %d %d\n", mid.numberoftriangles, 3, 0); 241 /* Remaining lines: <triangle #> <node> <node> <node> ... [attributes] */ 242 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]); } 243 fclose(file); 244 245 sprintf(fname, "C%d_%d.node", level, rank); 246 file = fopen(fname, "w"); 247 /* First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> */ 248 /* fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0); */ 249 fprintf(file, "%d %d %d %d\n", nPlotPts, 2, 0, 0); 250 /*Following lines: <vertex #> <x> <y> */ 251 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]); } 252 253 sid /= 2; 254 for (jj = 0; jj < nFineLoc; jj++) { 255 PetscBool sel = PETSC_TRUE; 256 for (kk = 0; kk < nselected_2 && sel; kk++) { 257 PetscInt lid = selected_idx_2[kk]; 258 if (lid == jj) sel = PETSC_FALSE; 259 } 260 if (sel) fprintf(file, "%d %e %e\n", sid++, coords[jj], coords[data_stride + jj]); 261 } 262 fclose(file); 263 PetscCheck(sid == nPlotPts, PETSC_COMM_SELF, PETSC_ERR_PLIB, "sid %d != nPlotPts %d", sid, nPlotPts); 264 level++; 265 } 266 } 267 { /* form P - setup some maps */ 268 PetscInt clid, mm, *nTri, *node_tri; 269 270 PetscCall(PetscMalloc2(nselected_2, &node_tri, nselected_2, &nTri)); 271 272 /* need list of triangles on node */ 273 for (kk = 0; kk < nselected_2; kk++) nTri[kk] = 0; 274 for (tid = 0, kk = 0; tid < mid.numberoftriangles; tid++) { 275 for (jj = 0; jj < 3; jj++) { 276 PetscInt cid = mid.trianglelist[kk++]; 277 if (nTri[cid] == 0) node_tri[cid] = tid; 278 nTri[cid]++; 279 } 280 } 281 #define EPS 1.e-12 282 /* find points and set prolongation */ 283 for (mm = clid = 0; mm < nFineLoc; mm++) { 284 PetscBool ise; 285 PetscCall(PetscCDEmptyAt(agg_lists_1, mm, &ise)); 286 if (!ise) { 287 const PetscInt lid = mm; 288 PetscScalar AA[3][3]; 289 PetscBLASInt N = 3, NRHS = 1, LDA = 3, IPIV[3], LDB = 3, INFO; 290 PetscCDIntNd *pos; 291 292 PetscCall(PetscCDGetHeadPos(agg_lists_1, lid, &pos)); 293 while (pos) { 294 PetscInt flid; 295 PetscCall(PetscCDIntNdGetID(pos, &flid)); 296 PetscCall(PetscCDGetNextPos(agg_lists_1, lid, &pos)); 297 298 if (flid < nFineLoc) { /* could be a ghost */ 299 PetscInt bestTID = -1; 300 PetscReal best_alpha = 1.e10; 301 const PetscInt fgid = flid + myFine0; 302 /* compute shape function for gid */ 303 const PetscReal fcoord[3] = {coords[flid], coords[data_stride + flid], 1.0}; 304 PetscBool haveit = PETSC_FALSE; 305 PetscScalar alpha[3]; 306 PetscInt clids[3]; 307 308 /* look for it */ 309 for (tid = node_tri[clid], jj = 0; jj < 5 && !haveit && tid != -1; jj++) { 310 for (tt = 0; tt < 3; tt++) { 311 PetscInt cid2 = mid.trianglelist[3 * tid + tt]; 312 PetscInt lid2 = selected_idx_2[cid2]; 313 AA[tt][0] = coords[lid2]; 314 AA[tt][1] = coords[data_stride + lid2]; 315 AA[tt][2] = 1.0; 316 clids[tt] = cid2; /* store for interp */ 317 } 318 319 for (tt = 0; tt < 3; tt++) alpha[tt] = (PetscScalar)fcoord[tt]; 320 321 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */ 322 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &INFO)); 323 { 324 PetscBool have = PETSC_TRUE; 325 PetscReal lowest = 1.e10; 326 for (tt = 0, idx = 0; tt < 3; tt++) { 327 if (PetscRealPart(alpha[tt]) > (1.0 + EPS) || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE; 328 if (PetscRealPart(alpha[tt]) < lowest) { 329 lowest = PetscRealPart(alpha[tt]); 330 idx = tt; 331 } 332 } 333 haveit = have; 334 } 335 tid = mid.neighborlist[3 * tid + idx]; 336 } 337 338 if (!haveit) { 339 /* brute force */ 340 for (tid = 0; tid < mid.numberoftriangles && !haveit; tid++) { 341 for (tt = 0; tt < 3; tt++) { 342 PetscInt cid2 = mid.trianglelist[3 * tid + tt]; 343 PetscInt lid2 = selected_idx_2[cid2]; 344 AA[tt][0] = coords[lid2]; 345 AA[tt][1] = coords[data_stride + lid2]; 346 AA[tt][2] = 1.0; 347 clids[tt] = cid2; /* store for interp */ 348 } 349 for (tt = 0; tt < 3; tt++) alpha[tt] = fcoord[tt]; 350 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */ 351 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &INFO)); 352 { 353 PetscBool have = PETSC_TRUE; 354 PetscReal worst = 0.0, v; 355 for (tt = 0; tt < 3 && have; tt++) { 356 if (PetscRealPart(alpha[tt]) > 1.0 + EPS || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE; 357 if ((v = PetscAbs(PetscRealPart(alpha[tt]) - 0.5)) > worst) worst = v; 358 } 359 if (worst < best_alpha) { 360 best_alpha = worst; 361 bestTID = tid; 362 } 363 haveit = have; 364 } 365 } 366 } 367 if (!haveit) { 368 if (best_alpha > *a_worst_best) *a_worst_best = best_alpha; 369 /* use best one */ 370 for (tt = 0; tt < 3; tt++) { 371 PetscInt cid2 = mid.trianglelist[3 * bestTID + tt]; 372 PetscInt lid2 = selected_idx_2[cid2]; 373 AA[tt][0] = coords[lid2]; 374 AA[tt][1] = coords[data_stride + lid2]; 375 AA[tt][2] = 1.0; 376 clids[tt] = cid2; /* store for interp */ 377 } 378 for (tt = 0; tt < 3; tt++) alpha[tt] = fcoord[tt]; 379 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */ 380 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&N, &NRHS, (PetscScalar *)AA, &LDA, IPIV, alpha, &LDB, &INFO)); 381 } 382 383 /* put in row of P */ 384 for (idx = 0; idx < 3; idx++) { 385 PetscScalar shp = alpha[idx]; 386 if (PetscAbs(PetscRealPart(shp)) > 1.e-6) { 387 PetscInt cgid = crsGID[clids[idx]]; 388 PetscInt jj = cgid * bs, ii = fgid * bs; /* need to gloalize */ 389 for (tt = 0; tt < bs; tt++, ii++, jj++) { PetscCall(MatSetValues(a_Prol, 1, &ii, 1, &jj, &shp, INSERT_VALUES)); } 390 } 391 } 392 } 393 } /* aggregates iterations */ 394 clid++; 395 } /* a coarse agg */ 396 } /* for all fine nodes */ 397 398 PetscCall(ISRestoreIndices(selected_2, &selected_idx_2)); 399 PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY)); 400 PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY)); 401 402 PetscCall(PetscFree2(node_tri, nTri)); 403 } 404 free(mid.trianglelist); 405 free(mid.neighborlist); 406 free(mid.segmentlist); 407 free(mid.segmentmarkerlist); 408 free(mid.pointlist); 409 free(mid.pointmarkerlist); 410 PetscCall(PetscFree(in.pointlist)); 411 PetscFunctionReturn(0); 412 #else 413 SETERRQ(PetscObjectComm((PetscObject)a_Prol), PETSC_ERR_PLIB, "configure with TRIANGLE to use geometric MG"); 414 #endif 415 } 416 /* -------------------------------------------------------------------------- */ 417 /* 418 getGIDsOnSquareGraph - square graph, get 419 420 Input Parameter: 421 . nselected_1 - selected local indices (includes ghosts in input Gmat1) 422 . clid_lid_1 - [nselected_1] lids of selected nodes 423 . Gmat1 - graph that goes with 'selected_1' 424 Output Parameter: 425 . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2) 426 . a_Gmat_2 - graph that is squared of 'Gmat_1' 427 . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes 428 */ 429 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) { 430 PetscMPIInt size; 431 PetscInt *crsGID, kk, my0, Iend, nloc; 432 MPI_Comm comm; 433 434 PetscFunctionBegin; 435 PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm)); 436 PetscCallMPI(MPI_Comm_size(comm, &size)); 437 PetscCall(MatGetOwnershipRange(Gmat1, &my0, &Iend)); /* AIJ */ 438 nloc = Iend - my0; /* this does not change */ 439 440 if (size == 1) { /* not much to do in serial */ 441 PetscCall(PetscMalloc1(nselected_1, &crsGID)); 442 for (kk = 0; kk < nselected_1; kk++) crsGID[kk] = kk; 443 *a_Gmat_2 = NULL; 444 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nselected_1, clid_lid_1, PETSC_COPY_VALUES, a_selected_2)); 445 } else { 446 PetscInt idx, num_fine_ghosts, num_crs_ghost, myCrs0; 447 Mat_MPIAIJ *mpimat2; 448 Mat Gmat2; 449 Vec locState; 450 PetscScalar *cpcol_state; 451 452 /* scan my coarse zero gid, set 'lid_state' with coarse GID */ 453 kk = nselected_1; 454 PetscCallMPI(MPI_Scan(&kk, &myCrs0, 1, MPIU_INT, MPI_SUM, comm)); 455 myCrs0 -= nselected_1; 456 457 if (a_Gmat_2) { /* output */ 458 /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */ 459 PetscCall(PCGAMGSquareGraph_GAMG(pc, Gmat1, &Gmat2)); 460 *a_Gmat_2 = Gmat2; /* output */ 461 } else Gmat2 = Gmat1; /* use local to get crsGIDs at least */ 462 /* get coarse grid GIDS for selected (locals and ghosts) */ 463 mpimat2 = (Mat_MPIAIJ *)Gmat2->data; 464 PetscCall(MatCreateVecs(Gmat2, &locState, NULL)); 465 PetscCall(VecSet(locState, (PetscScalar)(PetscReal)(-1))); /* set with UNKNOWN state */ 466 for (kk = 0; kk < nselected_1; kk++) { 467 PetscInt fgid = clid_lid_1[kk] + my0; 468 PetscScalar v = (PetscScalar)(kk + myCrs0); 469 PetscCall(VecSetValues(locState, 1, &fgid, &v, INSERT_VALUES)); /* set with PID */ 470 } 471 PetscCall(VecAssemblyBegin(locState)); 472 PetscCall(VecAssemblyEnd(locState)); 473 PetscCall(VecScatterBegin(mpimat2->Mvctx, locState, mpimat2->lvec, INSERT_VALUES, SCATTER_FORWARD)); 474 PetscCall(VecScatterEnd(mpimat2->Mvctx, locState, mpimat2->lvec, INSERT_VALUES, SCATTER_FORWARD)); 475 PetscCall(VecGetLocalSize(mpimat2->lvec, &num_fine_ghosts)); 476 PetscCall(VecGetArray(mpimat2->lvec, &cpcol_state)); 477 for (kk = 0, num_crs_ghost = 0; kk < num_fine_ghosts; kk++) { 478 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) num_crs_ghost++; 479 } 480 PetscCall(PetscMalloc1(nselected_1 + num_crs_ghost, &crsGID)); /* output */ 481 { 482 PetscInt *selected_set; 483 PetscCall(PetscMalloc1(nselected_1 + num_crs_ghost, &selected_set)); 484 /* do ghost of 'crsGID' */ 485 for (kk = 0, idx = nselected_1; kk < num_fine_ghosts; kk++) { 486 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) { 487 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]); 488 selected_set[idx] = nloc + kk; 489 crsGID[idx++] = cgid; 490 } 491 } 492 PetscCheck(idx == (nselected_1 + num_crs_ghost), PETSC_COMM_SELF, PETSC_ERR_PLIB, "idx %" PetscInt_FMT " != (nselected_1 %" PetscInt_FMT " + num_crs_ghost %" PetscInt_FMT ")", idx, nselected_1, num_crs_ghost); 493 PetscCall(VecRestoreArray(mpimat2->lvec, &cpcol_state)); 494 /* do locals in 'crsGID' */ 495 PetscCall(VecGetArray(locState, &cpcol_state)); 496 for (kk = 0, idx = 0; kk < nloc; kk++) { 497 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) { 498 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]); 499 selected_set[idx] = kk; 500 crsGID[idx++] = cgid; 501 } 502 } 503 PetscCheck(idx == nselected_1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "idx %" PetscInt_FMT " != nselected_1 %" PetscInt_FMT, idx, nselected_1); 504 PetscCall(VecRestoreArray(locState, &cpcol_state)); 505 506 if (a_selected_2 != NULL) { /* output */ 507 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (nselected_1 + num_crs_ghost), selected_set, PETSC_OWN_POINTER, a_selected_2)); 508 } else { 509 PetscCall(PetscFree(selected_set)); 510 } 511 } 512 PetscCall(VecDestroy(&locState)); 513 } 514 *a_crsGID = crsGID; /* output */ 515 PetscFunctionReturn(0); 516 } 517 518 /* -------------------------------------------------------------------------- */ 519 /* 520 PCGAMGGraph_GEO 521 522 Input Parameter: 523 . pc - this 524 . Amat - matrix on this fine level 525 Output Parameter: 526 . a_Gmat 527 */ 528 PetscErrorCode PCGAMGGraph_GEO(PC pc, Mat Amat, Mat *a_Gmat) { 529 PC_MG *mg = (PC_MG *)pc->data; 530 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 531 const PetscReal vfilter = pc_gamg->threshold[0]; 532 MPI_Comm comm; 533 Mat Gmat, F = NULL; 534 PetscBool set, flg, symm; 535 536 PetscFunctionBegin; 537 PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm)); 538 539 PetscCall(MatIsSymmetricKnown(Amat, &set, &flg)); 540 symm = (PetscBool) !(set && flg); 541 542 PetscCall(MatCreateGraph(Amat, symm, PETSC_TRUE, &Gmat)); 543 PetscCall(MatFilter(Gmat, vfilter, &F)); 544 if (F) { 545 PetscCall(MatDestroy(&Gmat)); 546 Gmat = F; 547 } 548 *a_Gmat = Gmat; 549 550 PetscFunctionReturn(0); 551 } 552 553 /* -------------------------------------------------------------------------- */ 554 /* 555 PCGAMGCoarsen_GEO 556 557 Input Parameter: 558 . a_pc - this 559 . a_Gmat - graph 560 Output Parameter: 561 . a_llist_parent - linked list from selected indices for data locality only 562 */ 563 PetscErrorCode PCGAMGCoarsen_GEO(PC a_pc, Mat *a_Gmat, PetscCoarsenData **a_llist_parent) { 564 PetscInt Istart, Iend, nloc, kk, Ii, ncols; 565 IS perm; 566 GAMGNode *gnodes; 567 PetscInt *permute; 568 Mat Gmat = *a_Gmat; 569 MPI_Comm comm; 570 MatCoarsen crs; 571 572 PetscFunctionBegin; 573 PetscCall(PetscObjectGetComm((PetscObject)a_pc, &comm)); 574 575 PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend)); 576 nloc = (Iend - Istart); 577 578 /* create random permutation with sort for geo-mg */ 579 PetscCall(PetscMalloc1(nloc, &gnodes)); 580 PetscCall(PetscMalloc1(nloc, &permute)); 581 582 for (Ii = Istart; Ii < Iend; Ii++) { /* locals only? */ 583 PetscCall(MatGetRow(Gmat, Ii, &ncols, NULL, NULL)); 584 { 585 PetscInt lid = Ii - Istart; 586 gnodes[lid].lid = lid; 587 gnodes[lid].degree = ncols; 588 } 589 PetscCall(MatRestoreRow(Gmat, Ii, &ncols, NULL, NULL)); 590 } 591 if (PETSC_TRUE) { 592 PetscRandom rand; 593 PetscBool *bIndexSet; 594 PetscReal rr; 595 PetscInt iSwapIndex; 596 597 PetscCall(PetscRandomCreate(comm, &rand)); 598 PetscCall(PetscCalloc1(nloc, &bIndexSet)); 599 for (Ii = 0; Ii < nloc; Ii++) { 600 PetscCall(PetscRandomGetValueReal(rand, &rr)); 601 iSwapIndex = (PetscInt)(rr * nloc); 602 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) { 603 GAMGNode iTemp = gnodes[iSwapIndex]; 604 gnodes[iSwapIndex] = gnodes[Ii]; 605 gnodes[Ii] = iTemp; 606 bIndexSet[Ii] = PETSC_TRUE; 607 bIndexSet[iSwapIndex] = PETSC_TRUE; 608 } 609 } 610 PetscCall(PetscRandomDestroy(&rand)); 611 PetscCall(PetscFree(bIndexSet)); 612 } 613 /* only sort locals */ 614 qsort(gnodes, nloc, sizeof(GAMGNode), petsc_geo_mg_compare); 615 /* create IS of permutation */ 616 for (kk = 0; kk < nloc; kk++) permute[kk] = gnodes[kk].lid; /* locals only */ 617 PetscCall(PetscFree(gnodes)); 618 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_OWN_POINTER, &perm)); 619 620 /* get MIS aggs */ 621 622 PetscCall(MatCoarsenCreate(comm, &crs)); 623 PetscCall(MatCoarsenSetType(crs, MATCOARSENMIS)); 624 PetscCall(MatCoarsenSetGreedyOrdering(crs, perm)); 625 PetscCall(MatCoarsenSetAdjacency(crs, Gmat)); 626 PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_FALSE)); 627 PetscCall(MatCoarsenApply(crs)); 628 PetscCall(MatCoarsenGetData(crs, a_llist_parent)); 629 PetscCall(MatCoarsenDestroy(&crs)); 630 631 PetscCall(ISDestroy(&perm)); 632 633 PetscFunctionReturn(0); 634 } 635 636 /* -------------------------------------------------------------------------- */ 637 /* 638 PCGAMGProlongator_GEO 639 640 Input Parameter: 641 . pc - this 642 . Amat - matrix on this fine level 643 . Graph - used to get ghost data for nodes in 644 . selected_1 - [nselected] 645 . agg_lists - [nselected] 646 Output Parameter: 647 . a_P_out - prolongation operator to the next level 648 */ 649 PetscErrorCode PCGAMGProlongator_GEO(PC pc, Mat Amat, Mat Gmat, PetscCoarsenData *agg_lists, Mat *a_P_out) { 650 PC_MG *mg = (PC_MG *)pc->data; 651 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 652 const PetscInt dim = pc_gamg->data_cell_cols, data_cols = pc_gamg->data_cell_cols; 653 PetscInt Istart, Iend, nloc, my0, jj, kk, ncols, nLocalSelected, bs, *clid_flid; 654 Mat Prol; 655 PetscMPIInt rank, size; 656 MPI_Comm comm; 657 IS selected_2, selected_1; 658 const PetscInt *selected_idx; 659 MatType mtype; 660 661 PetscFunctionBegin; 662 PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm)); 663 664 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 665 PetscCallMPI(MPI_Comm_size(comm, &size)); 666 PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend)); 667 PetscCall(MatGetBlockSize(Amat, &bs)); 668 nloc = (Iend - Istart) / bs; 669 my0 = Istart / bs; 670 PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") %% bs %" PetscInt_FMT, Iend, Istart, bs); 671 672 /* get 'nLocalSelected' */ 673 PetscCall(PetscCDGetMIS(agg_lists, &selected_1)); 674 PetscCall(ISGetSize(selected_1, &jj)); 675 PetscCall(PetscMalloc1(jj, &clid_flid)); 676 PetscCall(ISGetIndices(selected_1, &selected_idx)); 677 for (kk = 0, nLocalSelected = 0; kk < jj; kk++) { 678 PetscInt lid = selected_idx[kk]; 679 if (lid < nloc) { 680 PetscCall(MatGetRow(Gmat, lid + my0, &ncols, NULL, NULL)); 681 if (ncols > 1) clid_flid[nLocalSelected++] = lid; /* fiter out singletons */ 682 PetscCall(MatRestoreRow(Gmat, lid + my0, &ncols, NULL, NULL)); 683 } 684 } 685 PetscCall(ISRestoreIndices(selected_1, &selected_idx)); 686 PetscCall(ISDestroy(&selected_1)); /* this is selected_1 in serial */ 687 688 /* create prolongator matrix */ 689 PetscCall(MatGetType(Amat, &mtype)); 690 PetscCall(MatCreate(comm, &Prol)); 691 PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * bs, PETSC_DETERMINE, PETSC_DETERMINE)); 692 PetscCall(MatSetBlockSizes(Prol, bs, bs)); 693 PetscCall(MatSetType(Prol, mtype)); 694 PetscCall(MatSeqAIJSetPreallocation(Prol, 3 * data_cols, NULL)); 695 PetscCall(MatMPIAIJSetPreallocation(Prol, 3 * data_cols, NULL, 3 * data_cols, NULL)); 696 697 /* can get all points "removed" - but not on geomg */ 698 PetscCall(MatGetSize(Prol, &kk, &jj)); 699 if (!jj) { 700 PetscCall(PetscInfo(pc, "ERROE: no selected points on coarse grid\n")); 701 PetscCall(PetscFree(clid_flid)); 702 PetscCall(MatDestroy(&Prol)); 703 *a_P_out = NULL; /* out */ 704 PetscFunctionReturn(0); 705 } 706 707 { 708 PetscReal *coords; 709 PetscInt data_stride; 710 PetscInt *crsGID = NULL; 711 Mat Gmat2; 712 713 PetscCheck(dim == data_cols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dim %" PetscInt_FMT " != data_cols %" PetscInt_FMT, dim, data_cols); 714 /* grow ghost data for better coarse grid cover of fine grid */ 715 /* messy method, squares graph and gets some data */ 716 PetscCall(getGIDsOnSquareGraph(pc, nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID)); 717 /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */ 718 /* create global vector of coorindates in 'coords' */ 719 if (size > 1) { 720 PetscCall(PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords)); 721 } else { 722 coords = (PetscReal *)pc_gamg->data; 723 data_stride = pc_gamg->data_sz / pc_gamg->data_cell_cols; 724 } 725 PetscCall(MatDestroy(&Gmat2)); 726 727 /* triangulate */ 728 if (dim == 2) { 729 PetscReal metric, tm; 730 PetscCall(triangulateAndFormProl(selected_2, data_stride, coords, nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric)); 731 PetscCall(PetscFree(crsGID)); 732 733 /* clean up and create coordinates for coarse grid (output) */ 734 if (size > 1) PetscCall(PetscFree(coords)); 735 736 PetscCall(MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm)); 737 if (tm > 1.) { /* needs to be globalized - should not happen */ 738 PetscCall(PetscInfo(pc, " failed metric for coarse grid %e\n", (double)tm)); 739 PetscCall(MatDestroy(&Prol)); 740 } else if (metric > .0) { 741 PetscCall(PetscInfo(pc, "worst metric for coarse grid = %e\n", (double)metric)); 742 } 743 } else SETERRQ(comm, PETSC_ERR_PLIB, "3D not implemented for 'geo' AMG"); 744 { /* create next coords - output */ 745 PetscReal *crs_crds; 746 PetscCall(PetscMalloc1(dim * nLocalSelected, &crs_crds)); 747 for (kk = 0; kk < nLocalSelected; kk++) { /* grab local select nodes to promote - output */ 748 PetscInt lid = clid_flid[kk]; 749 for (jj = 0; jj < dim; jj++) crs_crds[jj * nLocalSelected + kk] = pc_gamg->data[jj * nloc + lid]; 750 } 751 752 PetscCall(PetscFree(pc_gamg->data)); 753 pc_gamg->data = crs_crds; /* out */ 754 pc_gamg->data_sz = dim * nLocalSelected; 755 } 756 PetscCall(ISDestroy(&selected_2)); 757 } 758 759 *a_P_out = Prol; /* out */ 760 PetscCall(PetscFree(clid_flid)); 761 762 PetscFunctionReturn(0); 763 } 764 765 static PetscErrorCode PCDestroy_GAMG_GEO(PC pc) { 766 PetscFunctionBegin; 767 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL)); 768 PetscFunctionReturn(0); 769 } 770 771 /* -------------------------------------------------------------------------- */ 772 /* 773 PCCreateGAMG_GEO 774 775 Input Parameter: 776 . pc - 777 */ 778 PetscErrorCode PCCreateGAMG_GEO(PC pc) { 779 PC_MG *mg = (PC_MG *)pc->data; 780 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 781 782 PetscFunctionBegin; 783 pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO; 784 pc_gamg->ops->destroy = PCDestroy_GAMG_GEO; 785 /* reset does not do anything; setup not virtual */ 786 787 /* set internal function pointers */ 788 pc_gamg->ops->graph = PCGAMGGraph_GEO; 789 pc_gamg->ops->coarsen = PCGAMGCoarsen_GEO; 790 pc_gamg->ops->prolongator = PCGAMGProlongator_GEO; 791 pc_gamg->ops->optprolongator = NULL; 792 pc_gamg->ops->createdefaultdata = PCSetData_GEO; 793 794 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_GEO)); 795 PetscFunctionReturn(0); 796 } 797