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