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