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