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