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