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 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); 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,F=NULL; 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(MatCreateGraph(Amat, symm, PETSC_TRUE, &Gmat)); 548 PetscCall(MatFilter(Gmat, vfilter, &F)); 549 if (F) { 550 PetscCall(MatDestroy(&Gmat)); 551 Gmat = F; 552 } 553 *a_Gmat = Gmat; 554 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 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 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 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 PetscCheck((Iend-Istart) % bs == 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") %% bs %" PetscInt_FMT,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 PetscCheck(dim == data_cols,PETSC_COMM_SELF,PETSC_ERR_PLIB,"dim %" PetscInt_FMT " != data_cols %" PetscInt_FMT,dim,data_cols); 720 /* grow ghost data for better coarse grid cover of fine grid */ 721 /* messy method, squares graph and gets some data */ 722 PetscCall(getGIDsOnSquareGraph(pc, nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID)); 723 /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */ 724 /* create global vector of coorindates in 'coords' */ 725 if (size > 1) { 726 PetscCall(PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords)); 727 } else { 728 coords = (PetscReal*)pc_gamg->data; 729 data_stride = pc_gamg->data_sz/pc_gamg->data_cell_cols; 730 } 731 PetscCall(MatDestroy(&Gmat2)); 732 733 /* triangulate */ 734 if (dim == 2) { 735 PetscReal metric,tm; 736 PetscCall(triangulateAndFormProl(selected_2, data_stride, coords,nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric)); 737 PetscCall(PetscFree(crsGID)); 738 739 /* clean up and create coordinates for coarse grid (output) */ 740 if (size > 1) PetscCall(PetscFree(coords)); 741 742 PetscCall(MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm)); 743 if (tm > 1.) { /* needs to be globalized - should not happen */ 744 PetscCall(PetscInfo(pc," failed metric for coarse grid %e\n",(double)tm)); 745 PetscCall(MatDestroy(&Prol)); 746 } else if (metric > .0) { 747 PetscCall(PetscInfo(pc,"worst metric for coarse grid = %e\n",(double)metric)); 748 } 749 } else SETERRQ(comm,PETSC_ERR_PLIB,"3D not implemented for 'geo' AMG"); 750 { /* create next coords - output */ 751 PetscReal *crs_crds; 752 PetscCall(PetscMalloc1(dim*nLocalSelected, &crs_crds)); 753 for (kk=0; kk<nLocalSelected; kk++) { /* grab local select nodes to promote - output */ 754 PetscInt lid = clid_flid[kk]; 755 for (jj=0; jj<dim; jj++) crs_crds[jj*nLocalSelected + kk] = pc_gamg->data[jj*nloc + lid]; 756 } 757 758 PetscCall(PetscFree(pc_gamg->data)); 759 pc_gamg->data = crs_crds; /* out */ 760 pc_gamg->data_sz = dim*nLocalSelected; 761 } 762 PetscCall(ISDestroy(&selected_2)); 763 } 764 765 *a_P_out = Prol; /* out */ 766 PetscCall(PetscFree(clid_flid)); 767 768 PetscFunctionReturn(0); 769 } 770 771 static PetscErrorCode PCDestroy_GAMG_GEO(PC pc) 772 { 773 PetscFunctionBegin; 774 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL)); 775 PetscFunctionReturn(0); 776 } 777 778 /* -------------------------------------------------------------------------- */ 779 /* 780 PCCreateGAMG_GEO 781 782 Input Parameter: 783 . pc - 784 */ 785 PetscErrorCode PCCreateGAMG_GEO(PC pc) 786 { 787 PC_MG *mg = (PC_MG*)pc->data; 788 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 789 790 PetscFunctionBegin; 791 pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO; 792 pc_gamg->ops->destroy = PCDestroy_GAMG_GEO; 793 /* reset does not do anything; setup not virtual */ 794 795 /* set internal function pointers */ 796 pc_gamg->ops->graph = PCGAMGGraph_GEO; 797 pc_gamg->ops->coarsen = PCGAMGCoarsen_GEO; 798 pc_gamg->ops->prolongator = PCGAMGProlongator_GEO; 799 pc_gamg->ops->optprolongator = NULL; 800 pc_gamg->ops->createdefaultdata = PCSetData_GEO; 801 802 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_GEO)); 803 PetscFunctionReturn(0); 804 } 805