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