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