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