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