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 %" PetscInt_FMT " must be %" PetscInt_FMT " or %" PetscInt_FMT ".",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 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]); 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 PetscOptionsHeadBegin(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 PetscOptionsHeadEnd(); 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 PetscCheck(sid == 2*nselected_2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %" PetscInt_FMT " != 2*nselected_2 %" PetscInt_FMT,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 PetscCheck(sid == nPlotPts,PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %" PetscInt_FMT " != nPlotPts %" PetscInt_FMT,sid,nPlotPts); 275 level++; 276 } 277 } 278 { /* form P - setup some maps */ 279 PetscInt clid,mm,*nTri,*node_tri; 280 281 PetscCall(PetscMalloc2(nselected_2, &node_tri,nselected_2, &nTri)); 282 283 /* need list of triangles on node */ 284 for (kk=0; kk<nselected_2; kk++) nTri[kk] = 0; 285 for (tid=0,kk=0; tid<mid.numberoftriangles; tid++) { 286 for (jj=0; jj<3; jj++) { 287 PetscInt cid = mid.trianglelist[kk++]; 288 if (nTri[cid] == 0) node_tri[cid] = tid; 289 nTri[cid]++; 290 } 291 } 292 #define EPS 1.e-12 293 /* find points and set prolongation */ 294 for (mm = clid = 0; mm < nFineLoc; mm++) { 295 PetscBool ise; 296 PetscCall(PetscCDEmptyAt(agg_lists_1,mm,&ise)); 297 if (!ise) { 298 const PetscInt lid = mm; 299 PetscScalar AA[3][3]; 300 PetscBLASInt N=3,NRHS=1,LDA=3,IPIV[3],LDB=3,INFO; 301 PetscCDIntNd *pos; 302 303 PetscCall(PetscCDGetHeadPos(agg_lists_1,lid,&pos)); 304 while (pos) { 305 PetscInt flid; 306 PetscCall(PetscCDIntNdGetID(pos, &flid)); 307 PetscCall(PetscCDGetNextPos(agg_lists_1,lid,&pos)); 308 309 if (flid < nFineLoc) { /* could be a ghost */ 310 PetscInt bestTID = -1; PetscReal best_alpha = 1.e10; 311 const PetscInt fgid = flid + myFine0; 312 /* compute shape function for gid */ 313 const PetscReal fcoord[3] = {coords[flid],coords[data_stride+flid],1.0}; 314 PetscBool haveit =PETSC_FALSE; PetscScalar alpha[3]; PetscInt clids[3]; 315 316 /* look for it */ 317 for (tid = node_tri[clid], jj=0; 318 jj < 5 && !haveit && tid != -1; 319 jj++) { 320 for (tt=0; tt<3; tt++) { 321 PetscInt cid2 = mid.trianglelist[3*tid + tt]; 322 PetscInt lid2 = selected_idx_2[cid2]; 323 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0; 324 clids[tt] = cid2; /* store for interp */ 325 } 326 327 for (tt=0; tt<3; tt++) alpha[tt] = (PetscScalar)fcoord[tt]; 328 329 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */ 330 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO)); 331 { 332 PetscBool have=PETSC_TRUE; PetscReal lowest=1.e10; 333 for (tt = 0, idx = 0; tt < 3; tt++) { 334 if (PetscRealPart(alpha[tt]) > (1.0+EPS) || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE; 335 if (PetscRealPart(alpha[tt]) < lowest) { 336 lowest = PetscRealPart(alpha[tt]); 337 idx = tt; 338 } 339 } 340 haveit = have; 341 } 342 tid = mid.neighborlist[3*tid + idx]; 343 } 344 345 if (!haveit) { 346 /* brute force */ 347 for (tid=0; tid<mid.numberoftriangles && !haveit; tid++) { 348 for (tt=0; tt<3; tt++) { 349 PetscInt cid2 = mid.trianglelist[3*tid + tt]; 350 PetscInt lid2 = selected_idx_2[cid2]; 351 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0; 352 clids[tt] = cid2; /* store for interp */ 353 } 354 for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt]; 355 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */ 356 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO)); 357 { 358 PetscBool have=PETSC_TRUE; PetscReal worst=0.0, v; 359 for (tt=0; tt<3 && have; tt++) { 360 if (PetscRealPart(alpha[tt]) > 1.0+EPS || PetscRealPart(alpha[tt]) < -EPS) have=PETSC_FALSE; 361 if ((v=PetscAbs(PetscRealPart(alpha[tt])-0.5)) > worst) worst = v; 362 } 363 if (worst < best_alpha) { 364 best_alpha = worst; bestTID = tid; 365 } 366 haveit = have; 367 } 368 } 369 } 370 if (!haveit) { 371 if (best_alpha > *a_worst_best) *a_worst_best = best_alpha; 372 /* use best one */ 373 for (tt=0; tt<3; tt++) { 374 PetscInt cid2 = mid.trianglelist[3*bestTID + tt]; 375 PetscInt lid2 = selected_idx_2[cid2]; 376 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0; 377 clids[tt] = cid2; /* store for interp */ 378 } 379 for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt]; 380 /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */ 381 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO)); 382 } 383 384 /* put in row of P */ 385 for (idx=0; idx<3; idx++) { 386 PetscScalar shp = alpha[idx]; 387 if (PetscAbs(PetscRealPart(shp)) > 1.e-6) { 388 PetscInt cgid = crsGID[clids[idx]]; 389 PetscInt jj = cgid*bs, ii = fgid*bs; /* need to gloalize */ 390 for (tt=0; tt < bs; tt++, ii++, jj++) { 391 PetscCall(MatSetValues(a_Prol,1,&ii,1,&jj,&shp,INSERT_VALUES)); 392 } 393 } 394 } 395 } 396 } /* aggregates iterations */ 397 clid++; 398 } /* a coarse agg */ 399 } /* for all fine nodes */ 400 401 PetscCall(ISRestoreIndices(selected_2, &selected_idx_2)); 402 PetscCall(MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY)); 403 PetscCall(MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY)); 404 405 PetscCall(PetscFree2(node_tri,nTri)); 406 } 407 free(mid.trianglelist); 408 free(mid.neighborlist); 409 free(mid.segmentlist); 410 free(mid.segmentmarkerlist); 411 free(mid.pointlist); 412 free(mid.pointmarkerlist); 413 PetscCall(PetscFree(in.pointlist)); 414 PetscFunctionReturn(0); 415 #else 416 SETERRQ(PetscObjectComm((PetscObject)a_Prol),PETSC_ERR_PLIB,"configure with TRIANGLE to use geometric MG"); 417 #endif 418 } 419 /* -------------------------------------------------------------------------- */ 420 /* 421 getGIDsOnSquareGraph - square graph, get 422 423 Input Parameter: 424 . nselected_1 - selected local indices (includes ghosts in input Gmat1) 425 . clid_lid_1 - [nselected_1] lids of selected nodes 426 . Gmat1 - graph that goes with 'selected_1' 427 Output Parameter: 428 . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2) 429 . a_Gmat_2 - graph that is squared of 'Gmat_1' 430 . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes 431 */ 432 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) 433 { 434 PetscMPIInt size; 435 PetscInt *crsGID, kk,my0,Iend,nloc; 436 MPI_Comm comm; 437 438 PetscFunctionBegin; 439 PetscCall(PetscObjectGetComm((PetscObject)Gmat1,&comm)); 440 PetscCallMPI(MPI_Comm_size(comm,&size)); 441 PetscCall(MatGetOwnershipRange(Gmat1,&my0,&Iend)); /* AIJ */ 442 nloc = Iend - my0; /* this does not change */ 443 444 if (size == 1) { /* not much to do in serial */ 445 PetscCall(PetscMalloc1(nselected_1, &crsGID)); 446 for (kk=0; kk<nselected_1; kk++) crsGID[kk] = kk; 447 *a_Gmat_2 = NULL; 448 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,nselected_1,clid_lid_1,PETSC_COPY_VALUES,a_selected_2)); 449 } else { 450 PetscInt idx,num_fine_ghosts,num_crs_ghost,myCrs0; 451 Mat_MPIAIJ *mpimat2; 452 Mat Gmat2; 453 Vec locState; 454 PetscScalar *cpcol_state; 455 456 /* scan my coarse zero gid, set 'lid_state' with coarse GID */ 457 kk = nselected_1; 458 PetscCallMPI(MPI_Scan(&kk, &myCrs0, 1, MPIU_INT, MPI_SUM, comm)); 459 myCrs0 -= nselected_1; 460 461 if (a_Gmat_2) { /* output */ 462 /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */ 463 PetscCall(PCGAMGSquareGraph_GAMG(pc,Gmat1,&Gmat2)); 464 *a_Gmat_2 = Gmat2; /* output */ 465 } else Gmat2 = Gmat1; /* use local to get crsGIDs at least */ 466 /* get coarse grid GIDS for selected (locals and ghosts) */ 467 mpimat2 = (Mat_MPIAIJ*)Gmat2->data; 468 PetscCall(MatCreateVecs(Gmat2, &locState, NULL)); 469 PetscCall(VecSet(locState, (PetscScalar)(PetscReal)(-1))); /* set with UNKNOWN state */ 470 for (kk=0; kk<nselected_1; kk++) { 471 PetscInt fgid = clid_lid_1[kk] + my0; 472 PetscScalar v = (PetscScalar)(kk+myCrs0); 473 PetscCall(VecSetValues(locState, 1, &fgid, &v, INSERT_VALUES)); /* set with PID */ 474 } 475 PetscCall(VecAssemblyBegin(locState)); 476 PetscCall(VecAssemblyEnd(locState)); 477 PetscCall(VecScatterBegin(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD)); 478 PetscCall(VecScatterEnd(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD)); 479 PetscCall(VecGetLocalSize(mpimat2->lvec, &num_fine_ghosts)); 480 PetscCall(VecGetArray(mpimat2->lvec, &cpcol_state)); 481 for (kk=0,num_crs_ghost=0; kk<num_fine_ghosts; kk++) { 482 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) num_crs_ghost++; 483 } 484 PetscCall(PetscMalloc1(nselected_1+num_crs_ghost, &crsGID)); /* output */ 485 { 486 PetscInt *selected_set; 487 PetscCall(PetscMalloc1(nselected_1+num_crs_ghost, &selected_set)); 488 /* do ghost of 'crsGID' */ 489 for (kk=0,idx=nselected_1; kk<num_fine_ghosts; kk++) { 490 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) { 491 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]); 492 selected_set[idx] = nloc + kk; 493 crsGID[idx++] = cgid; 494 } 495 } 496 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); 497 PetscCall(VecRestoreArray(mpimat2->lvec, &cpcol_state)); 498 /* do locals in 'crsGID' */ 499 PetscCall(VecGetArray(locState, &cpcol_state)); 500 for (kk=0,idx=0; kk<nloc; kk++) { 501 if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) { 502 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]); 503 selected_set[idx] = kk; 504 crsGID[idx++] = cgid; 505 } 506 } 507 PetscCheck(idx == nselected_1,PETSC_COMM_SELF,PETSC_ERR_PLIB,"idx %" PetscInt_FMT " != nselected_1 %" PetscInt_FMT,idx,nselected_1); 508 PetscCall(VecRestoreArray(locState, &cpcol_state)); 509 510 if (a_selected_2 != NULL) { /* output */ 511 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,(nselected_1+num_crs_ghost),selected_set,PETSC_OWN_POINTER,a_selected_2)); 512 } else { 513 PetscCall(PetscFree(selected_set)); 514 } 515 } 516 PetscCall(VecDestroy(&locState)); 517 } 518 *a_crsGID = crsGID; /* output */ 519 PetscFunctionReturn(0); 520 } 521 522 /* -------------------------------------------------------------------------- */ 523 /* 524 PCGAMGGraph_GEO 525 526 Input Parameter: 527 . pc - this 528 . Amat - matrix on this fine level 529 Output Parameter: 530 . a_Gmat 531 */ 532 PetscErrorCode PCGAMGGraph_GEO(PC pc,Mat Amat,Mat *a_Gmat) 533 { 534 PC_MG *mg = (PC_MG*)pc->data; 535 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 536 const PetscReal vfilter = pc_gamg->threshold[0]; 537 MPI_Comm comm; 538 Mat Gmat; 539 PetscBool set,flg,symm; 540 541 PetscFunctionBegin; 542 PetscCall(PetscObjectGetComm((PetscObject)Amat,&comm)); 543 544 PetscCall(MatIsSymmetricKnown(Amat, &set, &flg)); 545 symm = (PetscBool)!(set && flg); 546 547 PetscCall(PCGAMGCreateGraph(Amat, &Gmat)); 548 PetscCall(PCGAMGFilterGraph(&Gmat, vfilter, symm)); 549 550 *a_Gmat = Gmat; 551 552 PetscFunctionReturn(0); 553 } 554 555 /* -------------------------------------------------------------------------- */ 556 /* 557 PCGAMGCoarsen_GEO 558 559 Input Parameter: 560 . a_pc - this 561 . a_Gmat - graph 562 Output Parameter: 563 . a_llist_parent - linked list from selected indices for data locality only 564 */ 565 PetscErrorCode PCGAMGCoarsen_GEO(PC a_pc,Mat *a_Gmat,PetscCoarsenData **a_llist_parent) 566 { 567 PetscInt Istart,Iend,nloc,kk,Ii,ncols; 568 IS perm; 569 GAMGNode *gnodes; 570 PetscInt *permute; 571 Mat Gmat = *a_Gmat; 572 MPI_Comm comm; 573 MatCoarsen crs; 574 575 PetscFunctionBegin; 576 PetscCall(PetscObjectGetComm((PetscObject)a_pc,&comm)); 577 578 PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend)); 579 nloc = (Iend-Istart); 580 581 /* create random permutation with sort for geo-mg */ 582 PetscCall(PetscMalloc1(nloc, &gnodes)); 583 PetscCall(PetscMalloc1(nloc, &permute)); 584 585 for (Ii=Istart; Ii<Iend; Ii++) { /* locals only? */ 586 PetscCall(MatGetRow(Gmat,Ii,&ncols,NULL,NULL)); 587 { 588 PetscInt lid = Ii - Istart; 589 gnodes[lid].lid = lid; 590 gnodes[lid].degree = ncols; 591 } 592 PetscCall(MatRestoreRow(Gmat,Ii,&ncols,NULL,NULL)); 593 } 594 if (PETSC_TRUE) { 595 PetscRandom rand; 596 PetscBool *bIndexSet; 597 PetscReal rr; 598 PetscInt iSwapIndex; 599 600 PetscCall(PetscRandomCreate(comm,&rand)); 601 PetscCall(PetscCalloc1(nloc, &bIndexSet)); 602 for (Ii = 0; Ii < nloc; Ii++) { 603 PetscCall(PetscRandomGetValueReal(rand,&rr)); 604 iSwapIndex = (PetscInt) (rr*nloc); 605 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) { 606 GAMGNode iTemp = gnodes[iSwapIndex]; 607 gnodes[iSwapIndex] = gnodes[Ii]; 608 gnodes[Ii] = iTemp; 609 bIndexSet[Ii] = PETSC_TRUE; 610 bIndexSet[iSwapIndex] = PETSC_TRUE; 611 } 612 } 613 PetscCall(PetscRandomDestroy(&rand)); 614 PetscCall(PetscFree(bIndexSet)); 615 } 616 /* only sort locals */ 617 qsort(gnodes, nloc, sizeof(GAMGNode), petsc_geo_mg_compare); 618 /* create IS of permutation */ 619 for (kk=0; kk<nloc; kk++) permute[kk] = gnodes[kk].lid; /* locals only */ 620 PetscCall(PetscFree(gnodes)); 621 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_OWN_POINTER, &perm)); 622 623 /* get MIS aggs */ 624 625 PetscCall(MatCoarsenCreate(comm, &crs)); 626 PetscCall(MatCoarsenSetType(crs, MATCOARSENMIS)); 627 PetscCall(MatCoarsenSetGreedyOrdering(crs, perm)); 628 PetscCall(MatCoarsenSetAdjacency(crs, Gmat)); 629 PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_FALSE)); 630 PetscCall(MatCoarsenApply(crs)); 631 PetscCall(MatCoarsenGetData(crs, a_llist_parent)); 632 PetscCall(MatCoarsenDestroy(&crs)); 633 634 PetscCall(ISDestroy(&perm)); 635 636 PetscFunctionReturn(0); 637 } 638 639 /* -------------------------------------------------------------------------- */ 640 /* 641 PCGAMGProlongator_GEO 642 643 Input Parameter: 644 . pc - this 645 . Amat - matrix on this fine level 646 . Graph - used to get ghost data for nodes in 647 . selected_1 - [nselected] 648 . agg_lists - [nselected] 649 Output Parameter: 650 . a_P_out - prolongation operator to the next level 651 */ 652 PetscErrorCode PCGAMGProlongator_GEO(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out) 653 { 654 PC_MG *mg = (PC_MG*)pc->data; 655 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 656 const PetscInt dim = pc_gamg->data_cell_cols, data_cols = pc_gamg->data_cell_cols; 657 PetscInt Istart,Iend,nloc,my0,jj,kk,ncols,nLocalSelected,bs,*clid_flid; 658 Mat Prol; 659 PetscMPIInt rank, size; 660 MPI_Comm comm; 661 IS selected_2,selected_1; 662 const PetscInt *selected_idx; 663 MatType mtype; 664 665 PetscFunctionBegin; 666 PetscCall(PetscObjectGetComm((PetscObject)Amat,&comm)); 667 668 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 669 PetscCallMPI(MPI_Comm_size(comm,&size)); 670 PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend)); 671 PetscCall(MatGetBlockSize(Amat, &bs)); 672 nloc = (Iend-Istart)/bs; my0 = Istart/bs; 673 PetscCheck((Iend-Istart) % bs == 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") %% bs %" PetscInt_FMT,Iend,Istart,bs); 674 675 /* get 'nLocalSelected' */ 676 PetscCall(PetscCDGetMIS(agg_lists, &selected_1)); 677 PetscCall(ISGetSize(selected_1, &jj)); 678 PetscCall(PetscMalloc1(jj, &clid_flid)); 679 PetscCall(ISGetIndices(selected_1, &selected_idx)); 680 for (kk=0,nLocalSelected=0; kk<jj; kk++) { 681 PetscInt lid = selected_idx[kk]; 682 if (lid<nloc) { 683 PetscCall(MatGetRow(Gmat,lid+my0,&ncols,NULL,NULL)); 684 if (ncols>1) clid_flid[nLocalSelected++] = lid; /* fiter out singletons */ 685 PetscCall(MatRestoreRow(Gmat,lid+my0,&ncols,NULL,NULL)); 686 } 687 } 688 PetscCall(ISRestoreIndices(selected_1, &selected_idx)); 689 PetscCall(ISDestroy(&selected_1)); /* this is selected_1 in serial */ 690 691 /* create prolongator matrix */ 692 PetscCall(MatGetType(Amat,&mtype)); 693 PetscCall(MatCreate(comm, &Prol)); 694 PetscCall(MatSetSizes(Prol,nloc*bs,nLocalSelected*bs,PETSC_DETERMINE,PETSC_DETERMINE)); 695 PetscCall(MatSetBlockSizes(Prol, bs, bs)); 696 PetscCall(MatSetType(Prol, mtype)); 697 PetscCall(MatSeqAIJSetPreallocation(Prol,3*data_cols,NULL)); 698 PetscCall(MatMPIAIJSetPreallocation(Prol,3*data_cols,NULL,3*data_cols,NULL)); 699 700 /* can get all points "removed" - but not on geomg */ 701 PetscCall(MatGetSize(Prol, &kk, &jj)); 702 if (!jj) { 703 PetscCall(PetscInfo(pc,"ERROE: no selected points on coarse grid\n")); 704 PetscCall(PetscFree(clid_flid)); 705 PetscCall(MatDestroy(&Prol)); 706 *a_P_out = NULL; /* out */ 707 PetscFunctionReturn(0); 708 } 709 710 { 711 PetscReal *coords; 712 PetscInt data_stride; 713 PetscInt *crsGID = NULL; 714 Mat Gmat2; 715 716 PetscCheck(dim == data_cols,PETSC_COMM_SELF,PETSC_ERR_PLIB,"dim %" PetscInt_FMT " != data_cols %" PetscInt_FMT,dim,data_cols); 717 /* grow ghost data for better coarse grid cover of fine grid */ 718 /* messy method, squares graph and gets some data */ 719 PetscCall(getGIDsOnSquareGraph(pc, nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID)); 720 /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */ 721 /* create global vector of coorindates in 'coords' */ 722 if (size > 1) { 723 PetscCall(PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords)); 724 } else { 725 coords = (PetscReal*)pc_gamg->data; 726 data_stride = pc_gamg->data_sz/pc_gamg->data_cell_cols; 727 } 728 PetscCall(MatDestroy(&Gmat2)); 729 730 /* triangulate */ 731 if (dim == 2) { 732 PetscReal metric,tm; 733 PetscCall(triangulateAndFormProl(selected_2, data_stride, coords,nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric)); 734 PetscCall(PetscFree(crsGID)); 735 736 /* clean up and create coordinates for coarse grid (output) */ 737 if (size > 1) PetscCall(PetscFree(coords)); 738 739 PetscCall(MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm)); 740 if (tm > 1.) { /* needs to be globalized - should not happen */ 741 PetscCall(PetscInfo(pc," failed metric for coarse grid %e\n",(double)tm)); 742 PetscCall(MatDestroy(&Prol)); 743 } else if (metric > .0) { 744 PetscCall(PetscInfo(pc,"worst metric for coarse grid = %e\n",(double)metric)); 745 } 746 } else SETERRQ(comm,PETSC_ERR_PLIB,"3D not implemented for 'geo' AMG"); 747 { /* create next coords - output */ 748 PetscReal *crs_crds; 749 PetscCall(PetscMalloc1(dim*nLocalSelected, &crs_crds)); 750 for (kk=0; kk<nLocalSelected; kk++) { /* grab local select nodes to promote - output */ 751 PetscInt lid = clid_flid[kk]; 752 for (jj=0; jj<dim; jj++) crs_crds[jj*nLocalSelected + kk] = pc_gamg->data[jj*nloc + lid]; 753 } 754 755 PetscCall(PetscFree(pc_gamg->data)); 756 pc_gamg->data = crs_crds; /* out */ 757 pc_gamg->data_sz = dim*nLocalSelected; 758 } 759 PetscCall(ISDestroy(&selected_2)); 760 } 761 762 *a_P_out = Prol; /* out */ 763 PetscCall(PetscFree(clid_flid)); 764 765 PetscFunctionReturn(0); 766 } 767 768 static PetscErrorCode PCDestroy_GAMG_GEO(PC pc) 769 { 770 PetscFunctionBegin; 771 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL)); 772 PetscFunctionReturn(0); 773 } 774 775 /* -------------------------------------------------------------------------- */ 776 /* 777 PCCreateGAMG_GEO 778 779 Input Parameter: 780 . pc - 781 */ 782 PetscErrorCode PCCreateGAMG_GEO(PC pc) 783 { 784 PC_MG *mg = (PC_MG*)pc->data; 785 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 786 787 PetscFunctionBegin; 788 pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO; 789 pc_gamg->ops->destroy = PCDestroy_GAMG_GEO; 790 /* reset does not do anything; setup not virtual */ 791 792 /* set internal function pointers */ 793 pc_gamg->ops->graph = PCGAMGGraph_GEO; 794 pc_gamg->ops->coarsen = PCGAMGCoarsen_GEO; 795 pc_gamg->ops->prolongator = PCGAMGProlongator_GEO; 796 pc_gamg->ops->optprolongator = NULL; 797 pc_gamg->ops->createdefaultdata = PCSetData_GEO; 798 799 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_GEO)); 800 PetscFunctionReturn(0); 801 } 802