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