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