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 /* -------------------------------------------------------------------------- */ 17 /* 18 PCSetCoordinates_GEO 19 20 Input Parameter: 21 . pc - the preconditioner context 22 */ 23 EXTERN_C_BEGIN 24 #undef __FUNCT__ 25 #define __FUNCT__ "PCSetCoordinates_GEO" 26 PetscErrorCode PCSetCoordinates_GEO( PC pc, PetscInt ndm, PetscReal *coords ) 27 { 28 PC_MG *mg = (PC_MG*)pc->data; 29 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 30 PetscErrorCode ierr; 31 PetscInt arrsz,bs,my0,kk,ii,nloc,Iend; 32 Mat Amat = pc->pmat; 33 34 PetscFunctionBegin; 35 PetscValidHeaderSpecific( Amat, MAT_CLASSID, 1 ); 36 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ( ierr ); 37 ierr = MatGetOwnershipRange( Amat, &my0, &Iend ); CHKERRQ(ierr); 38 nloc = (Iend-my0)/bs; 39 if((Iend-my0)%bs!=0) SETERRQ1(((PetscObject)Amat)->comm,PETSC_ERR_ARG_WRONG, "Bad local size %d.",nloc); 40 41 pc_gamg->data_rows = 1; 42 if( coords==0 && nloc > 0 ) { 43 SETERRQ(((PetscObject)Amat)->comm,PETSC_ERR_ARG_WRONG, "Need coordinates for pc_gamg_type 'geo'."); 44 } 45 pc_gamg->data_cols = ndm; /* coordinates */ 46 47 arrsz = nloc*pc_gamg->data_rows*pc_gamg->data_cols; 48 49 /* create data - syntactic sugar that should be refactored at some point */ 50 if (pc_gamg->data==0 || (pc_gamg->data_sz != arrsz)) { 51 ierr = PetscFree( pc_gamg->data ); CHKERRQ(ierr); 52 ierr = PetscMalloc((arrsz+1)*sizeof(PetscReal), &pc_gamg->data ); CHKERRQ(ierr); 53 } 54 for(kk=0;kk<arrsz;kk++)pc_gamg->data[kk] = -999.; 55 pc_gamg->data[arrsz] = -99.; 56 /* copy data in - column oriented */ 57 for( kk = 0 ; kk < nloc ; kk++ ){ 58 for( ii = 0 ; ii < ndm ; ii++ ) { 59 pc_gamg->data[ii*nloc + kk] = coords[kk*ndm + ii]; 60 } 61 } 62 assert(pc_gamg->data[arrsz] == -99.); 63 64 pc_gamg->data_sz = arrsz; 65 66 PetscFunctionReturn(0); 67 } 68 EXTERN_C_END 69 70 /* -------------------------------------------------------------------------- */ 71 /* 72 PCSetData_GEO 73 74 Input Parameter: 75 . pc - 76 */ 77 #undef __FUNCT__ 78 #define __FUNCT__ "PCSetData_GEO" 79 PetscErrorCode PCSetData_GEO( PC pc ) 80 { 81 PetscFunctionBegin; 82 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_LIB,"GEO MG needs coordinates"); 83 } 84 85 /* -------------------------------------------------------------------------- */ 86 /* 87 PCSetFromOptions_GEO 88 89 Input Parameter: 90 . pc - 91 */ 92 #undef __FUNCT__ 93 #define __FUNCT__ "PCSetFromOptions_GEO" 94 PetscErrorCode PCSetFromOptions_GEO( PC pc ) 95 { 96 PetscErrorCode ierr; 97 PC_MG *mg = (PC_MG*)pc->data; 98 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 99 100 PetscFunctionBegin; 101 ierr = PetscOptionsHead("GAMG-GEO options"); CHKERRQ(ierr); 102 { 103 /* -pc_gamg_sa_nsmooths */ 104 /* pc_gamg_sa->smooths = 0; */ 105 /* ierr = PetscOptionsInt("-pc_gamg_agg_nsmooths", */ 106 /* "smoothing steps for smoothed aggregation, usually 1 (0)", */ 107 /* "PCGAMGSetNSmooths_AGG", */ 108 /* pc_gamg_sa->smooths, */ 109 /* &pc_gamg_sa->smooths, */ 110 /* &flag); */ 111 /* CHKERRQ(ierr); */ 112 } 113 ierr = PetscOptionsTail();CHKERRQ(ierr); 114 115 /* call base class */ 116 ierr = PCSetFromOptions_GAMG( pc ); CHKERRQ(ierr); 117 118 if( pc_gamg->verbose ) { 119 PetscPrintf(PETSC_COMM_WORLD,"[%d]%s done\n",0,__FUNCT__); 120 } 121 122 PetscFunctionReturn(0); 123 } 124 125 /* -------------------------------------------------------------------------- */ 126 /* 127 PCCreateGAMG_GEO - nothing to do here now 128 129 Input Parameter: 130 . pc - 131 */ 132 #undef __FUNCT__ 133 #define __FUNCT__ "PCCreateGAMG_GEO" 134 PetscErrorCode PCCreateGAMG_GEO( PC pc ) 135 { 136 PetscErrorCode ierr; 137 PC_MG *mg = (PC_MG*)pc->data; 138 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 139 140 PetscFunctionBegin; 141 pc->ops->setfromoptions = PCSetFromOptions_GEO; 142 /* pc->ops->destroy = PCDestroy_GEO; */ 143 /* reset does not do anything; setup not virtual */ 144 145 /* set internal function pointers */ 146 pc_gamg->createprolongator = PCGAMGcreateProl_GEO; 147 pc_gamg->createdefaultdata = PCSetData_GEO; 148 149 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 150 "PCSetCoordinates_C", 151 "PCSetCoordinates_GEO", 152 PCSetCoordinates_GEO);CHKERRQ(ierr); 153 154 PetscFunctionReturn(0); 155 } 156 157 /* -------------------------------------------------------------------------- */ 158 /* 159 triangulateAndFormProl 160 161 Input Parameter: 162 . selected_2 - list of selected local ID, includes selected ghosts 163 . nnodes - 164 . coords[2*nnodes] - column vector of local coordinates w/ ghosts 165 . selected_1 - selected IDs that go with base (1) graph 166 . locals_llist - linked list with (some) locality info of base graph 167 . crsGID[selected.size()] - global index for prolongation operator 168 . bs - block size 169 Output Parameter: 170 . a_Prol - prolongation operator 171 . a_worst_best - measure of worst missed fine vertex, 0 is no misses 172 */ 173 #undef __FUNCT__ 174 #define __FUNCT__ "triangulateAndFormProl" 175 PetscErrorCode triangulateAndFormProl( IS selected_2, /* list of selected local ID, includes selected ghosts */ 176 const PetscInt nnodes, 177 const PetscReal coords[], /* column vector of local coordinates w/ ghosts */ 178 IS selected_1, /* list of selected local ID, includes selected ghosts */ 179 IS locals_llist, /* linked list from selected vertices of aggregate unselected vertices */ 180 const PetscInt crsGID[], 181 const PetscInt bs, 182 Mat a_Prol, /* prolongation operator (output) */ 183 PetscReal *a_worst_best /* measure of worst missed fine vertex, 0 is no misses */ 184 ) 185 { 186 #if defined(PETSC_HAVE_TRIANGLE) 187 PetscErrorCode ierr; 188 PetscInt jj,tid,tt,idx,nselected_1,nselected_2; 189 struct triangulateio in,mid; 190 const PetscInt *selected_idx_1,*selected_idx_2,*llist_idx; 191 PetscMPIInt mype,npe; 192 PetscInt Istart,Iend,nFineLoc,myFine0; 193 int kk,nPlotPts,sid; 194 195 PetscFunctionBegin; 196 *a_worst_best = 0.0; 197 ierr = MPI_Comm_rank(((PetscObject)a_Prol)->comm,&mype); CHKERRQ(ierr); 198 ierr = MPI_Comm_size(((PetscObject)a_Prol)->comm,&npe); CHKERRQ(ierr); 199 ierr = ISGetLocalSize( selected_1, &nselected_1 ); CHKERRQ(ierr); 200 ierr = ISGetLocalSize( selected_2, &nselected_2 ); CHKERRQ(ierr); 201 if(nselected_2 == 1 || nselected_2 == 2 ){ /* 0 happens on idle processors */ 202 /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Not enough points - error in stopping logic",nselected_2); */ 203 *a_worst_best = 100.0; /* this will cause a stop, but not globalized (should not happen) */ 204 #ifdef VERBOSE 205 PetscPrintf(PETSC_COMM_SELF,"[%d]%s %d selected point - bailing out\n",mype,__FUNCT__,nselected_2); 206 #endif 207 PetscFunctionReturn(0); 208 } 209 ierr = MatGetOwnershipRange( a_Prol, &Istart, &Iend ); CHKERRQ(ierr); 210 nFineLoc = (Iend-Istart)/bs; myFine0 = Istart/bs; 211 nPlotPts = nFineLoc; /* locals */ 212 /* traingle */ 213 /* Define input points - in*/ 214 in.numberofpoints = nselected_2; 215 in.numberofpointattributes = 0; 216 /* get nselected points */ 217 ierr = PetscMalloc( 2*(nselected_2)*sizeof(REAL), &in.pointlist ); CHKERRQ(ierr); 218 ierr = ISGetIndices( selected_2, &selected_idx_2 ); CHKERRQ(ierr); 219 220 for(kk=0,sid=0;kk<nselected_2;kk++,sid += 2){ 221 PetscInt lid = selected_idx_2[kk]; 222 in.pointlist[sid] = coords[lid]; 223 in.pointlist[sid+1] = coords[nnodes + lid]; 224 if(lid>=nFineLoc) nPlotPts++; 225 } 226 assert(sid==2*nselected_2); 227 228 in.numberofsegments = 0; 229 in.numberofedges = 0; 230 in.numberofholes = 0; 231 in.numberofregions = 0; 232 in.trianglelist = 0; 233 in.segmentmarkerlist = 0; 234 in.pointattributelist = 0; 235 in.pointmarkerlist = 0; 236 in.triangleattributelist = 0; 237 in.trianglearealist = 0; 238 in.segmentlist = 0; 239 in.holelist = 0; 240 in.regionlist = 0; 241 in.edgelist = 0; 242 in.edgemarkerlist = 0; 243 in.normlist = 0; 244 /* triangulate */ 245 mid.pointlist = 0; /* Not needed if -N switch used. */ 246 /* Not needed if -N switch used or number of point attributes is zero: */ 247 mid.pointattributelist = 0; 248 mid.pointmarkerlist = 0; /* Not needed if -N or -B switch used. */ 249 mid.trianglelist = 0; /* Not needed if -E switch used. */ 250 /* Not needed if -E switch used or number of triangle attributes is zero: */ 251 mid.triangleattributelist = 0; 252 mid.neighborlist = 0; /* Needed only if -n switch used. */ 253 /* Needed only if segments are output (-p or -c) and -P not used: */ 254 mid.segmentlist = 0; 255 /* Needed only if segments are output (-p or -c) and -P and -B not used: */ 256 mid.segmentmarkerlist = 0; 257 mid.edgelist = 0; /* Needed only if -e switch used. */ 258 mid.edgemarkerlist = 0; /* Needed if -e used and -B not used. */ 259 mid.numberoftriangles = 0; 260 261 /* Triangulate the points. Switches are chosen to read and write a */ 262 /* PSLG (p), preserve the convex hull (c), number everything from */ 263 /* zero (z), assign a regional attribute to each element (A), and */ 264 /* produce an edge list (e), a Voronoi diagram (v), and a triangle */ 265 /* neighbor list (n). */ 266 if(nselected_2 != 0){ /* inactive processor */ 267 char args[] = "npczQ"; /* c is needed ? */ 268 triangulate(args, &in, &mid, (struct triangulateio *) NULL ); 269 /* output .poly files for 'showme' */ 270 if( !PETSC_TRUE ) { 271 static int level = 1; 272 FILE *file; char fname[32]; 273 274 sprintf(fname,"C%d_%d.poly",level,mype); file = fopen(fname, "w"); 275 /*First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>*/ 276 fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0); 277 /*Following lines: <vertex #> <x> <y> */ 278 for(kk=0,sid=0;kk<in.numberofpoints;kk++,sid += 2){ 279 fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]); 280 } 281 /*One line: <# of segments> <# of boundary markers (0 or 1)> */ 282 fprintf(file, "%d %d\n",0,0); 283 /*Following lines: <segment #> <endpoint> <endpoint> [boundary marker] */ 284 /* One line: <# of holes> */ 285 fprintf(file, "%d\n",0); 286 /* Following lines: <hole #> <x> <y> */ 287 /* Optional line: <# of regional attributes and/or area constraints> */ 288 /* Optional following lines: <region #> <x> <y> <attribute> <maximum area> */ 289 fclose(file); 290 291 /* elems */ 292 sprintf(fname,"C%d_%d.ele",level,mype); file = fopen(fname, "w"); 293 /* First line: <# of triangles> <nodes per triangle> <# of attributes> */ 294 fprintf(file, "%d %d %d\n",mid.numberoftriangles,3,0); 295 /* Remaining lines: <triangle #> <node> <node> <node> ... [attributes] */ 296 for(kk=0,sid=0;kk<mid.numberoftriangles;kk++,sid += 3){ 297 fprintf(file, "%d %d %d %d\n",kk,mid.trianglelist[sid],mid.trianglelist[sid+1],mid.trianglelist[sid+2]); 298 } 299 fclose(file); 300 301 sprintf(fname,"C%d_%d.node",level,mype); file = fopen(fname, "w"); 302 /* First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> */ 303 /* fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0); */ 304 fprintf(file, "%d %d %d %d\n",nPlotPts,2,0,0); 305 /*Following lines: <vertex #> <x> <y> */ 306 for(kk=0,sid=0;kk<in.numberofpoints;kk++,sid+=2){ 307 fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]); 308 } 309 310 sid /= 2; 311 for(jj=0;jj<nFineLoc;jj++){ 312 PetscBool sel = PETSC_TRUE; 313 for( kk=0 ; kk<nselected_2 && sel ; kk++ ){ 314 PetscInt lid = selected_idx_2[kk]; 315 if( lid == jj ) sel = PETSC_FALSE; 316 } 317 if( sel ) { 318 fprintf(file, "%d %e %e\n",sid++,coords[jj],coords[nnodes + jj]); 319 } 320 } 321 fclose(file); 322 assert(sid==nPlotPts); 323 level++; 324 } 325 } 326 #if defined PETSC_USE_LOG 327 ierr = PetscLogEventBegin(gamg_setup_events[FIND_V],0,0,0,0);CHKERRQ(ierr); 328 #endif 329 { /* form P - setup some maps */ 330 PetscInt clid_iterator; 331 PetscInt *nTri, *node_tri; 332 333 ierr = PetscMalloc( nselected_2*sizeof(PetscInt), &node_tri ); CHKERRQ(ierr); 334 ierr = PetscMalloc( nselected_2*sizeof(PetscInt), &nTri ); CHKERRQ(ierr); 335 336 /* need list of triangles on node */ 337 for(kk=0;kk<nselected_2;kk++) nTri[kk] = 0; 338 for(tid=0,kk=0;tid<mid.numberoftriangles;tid++){ 339 for(jj=0;jj<3;jj++) { 340 PetscInt cid = mid.trianglelist[kk++]; 341 if( nTri[cid] == 0 ) node_tri[cid] = tid; 342 nTri[cid]++; 343 } 344 } 345 #define EPS 1.e-12 346 /* find points and set prolongation */ 347 ierr = ISGetIndices( selected_1, &selected_idx_1 ); CHKERRQ(ierr); 348 ierr = ISGetIndices( locals_llist, &llist_idx ); CHKERRQ(ierr); 349 for( clid_iterator = 0 ; clid_iterator < nselected_1 ; clid_iterator++ ){ 350 PetscInt flid = selected_idx_1[clid_iterator]; assert(flid != -1); 351 PetscScalar AA[3][3]; 352 PetscBLASInt N=3,NRHS=1,LDA=3,IPIV[3],LDB=3,INFO; 353 do{ 354 if( flid < nFineLoc ) { /* could be a ghost */ 355 PetscInt bestTID = -1; PetscReal best_alpha = 1.e10; 356 const PetscInt fgid = flid + myFine0; 357 /* compute shape function for gid */ 358 const PetscReal fcoord[3] = { coords[flid], coords[nnodes + flid], 1.0 }; 359 PetscBool haveit = PETSC_FALSE; PetscScalar alpha[3]; PetscInt clids[3]; 360 /* look for it */ 361 for( tid = node_tri[clid_iterator], jj=0; 362 jj < 5 && !haveit && tid != -1; 363 jj++ ){ 364 for(tt=0;tt<3;tt++){ 365 PetscInt cid2 = mid.trianglelist[3*tid + tt]; 366 PetscInt lid2 = selected_idx_2[cid2]; 367 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[nnodes + lid2]; AA[tt][2] = 1.0; 368 clids[tt] = cid2; /* store for interp */ 369 } 370 371 for(tt=0;tt<3;tt++) alpha[tt] = (PetscScalar)fcoord[tt]; 372 373 /* SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) */ 374 LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO); 375 { 376 PetscBool have=PETSC_TRUE; PetscReal lowest=1.e10; 377 for( tt = 0, idx = 0 ; tt < 3 ; tt++ ) { 378 if( PetscRealPart(alpha[tt]) > (1.0+EPS) || PetscRealPart(alpha[tt]) < -EPS ) have = PETSC_FALSE; 379 if( PetscRealPart(alpha[tt]) < lowest ){ 380 lowest = PetscRealPart(alpha[tt]); 381 idx = tt; 382 } 383 } 384 haveit = have; 385 } 386 tid = mid.neighborlist[3*tid + idx]; 387 } 388 389 if( !haveit ) { 390 /* brute force */ 391 for(tid=0 ; tid<mid.numberoftriangles && !haveit ; tid++ ){ 392 for(tt=0;tt<3;tt++){ 393 PetscInt cid2 = mid.trianglelist[3*tid + tt]; 394 PetscInt lid2 = selected_idx_2[cid2]; 395 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[nnodes + lid2]; AA[tt][2] = 1.0; 396 clids[tt] = cid2; /* store for interp */ 397 } 398 for(tt=0;tt<3;tt++) alpha[tt] = fcoord[tt]; 399 /* SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) */ 400 LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO); 401 { 402 PetscBool have=PETSC_TRUE; PetscReal worst=0.0, v; 403 for(tt=0; tt<3 && have ;tt++) { 404 if( PetscRealPart(alpha[tt]) > 1.0+EPS || PetscRealPart(alpha[tt]) < -EPS ) have=PETSC_FALSE; 405 if( (v=PetscAbs(PetscRealPart(alpha[tt])-0.5)) > worst ) worst = v; 406 } 407 if( worst < best_alpha ) { 408 best_alpha = worst; bestTID = tid; 409 } 410 haveit = have; 411 } 412 } 413 } 414 if( !haveit ) { 415 if( best_alpha > *a_worst_best ) *a_worst_best = best_alpha; 416 /* use best one */ 417 for(tt=0;tt<3;tt++){ 418 PetscInt cid2 = mid.trianglelist[3*bestTID + tt]; 419 PetscInt lid2 = selected_idx_2[cid2]; 420 AA[tt][0] = coords[lid2]; AA[tt][1] = coords[nnodes + lid2]; AA[tt][2] = 1.0; 421 clids[tt] = cid2; /* store for interp */ 422 } 423 for(tt=0;tt<3;tt++) alpha[tt] = fcoord[tt]; 424 /* SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) */ 425 LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO); 426 } 427 428 /* put in row of P */ 429 for(idx=0;idx<3;idx++){ 430 PetscScalar shp = alpha[idx]; 431 if( PetscAbs(PetscRealPart(shp)) > 1.e-6 ) { 432 PetscInt cgid = crsGID[clids[idx]]; 433 PetscInt jj = cgid*bs, ii = fgid*bs; /* need to gloalize */ 434 for(tt=0 ; tt < bs ; tt++, ii++, jj++ ){ 435 ierr = MatSetValues(a_Prol,1,&ii,1,&jj,&shp,INSERT_VALUES); CHKERRQ(ierr); 436 } 437 } 438 } 439 } /* local vertex test */ 440 } while( (flid=llist_idx[flid]) != -1 ); 441 } 442 ierr = ISRestoreIndices( selected_2, &selected_idx_2 ); CHKERRQ(ierr); 443 ierr = ISRestoreIndices( selected_1, &selected_idx_1 ); CHKERRQ(ierr); 444 ierr = ISRestoreIndices( locals_llist, &llist_idx ); CHKERRQ(ierr); 445 ierr = MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 446 ierr = MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 447 448 ierr = PetscFree( node_tri ); CHKERRQ(ierr); 449 ierr = PetscFree( nTri ); CHKERRQ(ierr); 450 } 451 #if defined PETSC_USE_LOG 452 ierr = PetscLogEventEnd(gamg_setup_events[FIND_V],0,0,0,0);CHKERRQ(ierr); 453 #endif 454 free( mid.trianglelist ); 455 free( mid.neighborlist ); 456 ierr = PetscFree( in.pointlist ); CHKERRQ(ierr); 457 458 PetscFunctionReturn(0); 459 #else 460 SETERRQ(((PetscObject)a_Prol)->comm,PETSC_ERR_LIB,"configure with TRIANGLE to use geometric MG"); 461 #endif 462 } 463 /* -------------------------------------------------------------------------- */ 464 /* 465 getGIDsOnSquareGraph - square graph, get 466 467 Input Parameter: 468 . selected_1 - selected local indices (includes ghosts in input Gmat_1) 469 . Gmat1 - graph that goes with 'selected_1' 470 Output Parameter: 471 . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2) 472 . a_Gmat_2 - graph that is squared of 'Gmat_1' 473 . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes 474 */ 475 #undef __FUNCT__ 476 #define __FUNCT__ "getGIDsOnSquareGraph" 477 PetscErrorCode getGIDsOnSquareGraph( const IS selected_1, 478 const Mat Gmat1, 479 IS *a_selected_2, 480 Mat *a_Gmat_2, 481 PetscInt **a_crsGID 482 ) 483 { 484 PetscErrorCode ierr; 485 PetscMPIInt mype,npe; 486 PetscInt *crsGID, kk,my0,Iend,nloc,nSelected_1; 487 const PetscInt *selected_idx; 488 MPI_Comm wcomm = ((PetscObject)Gmat1)->comm; 489 490 PetscFunctionBegin; 491 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 492 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 493 ierr = MatGetOwnershipRange(Gmat1,&my0,&Iend); CHKERRQ(ierr); /* AIJ */ 494 nloc = Iend - my0; /* this does not change */ 495 ierr = ISGetLocalSize( selected_1, &nSelected_1 ); CHKERRQ(ierr); 496 497 if (npe == 1) { /* not much to do in serial */ 498 ierr = PetscMalloc( nSelected_1*sizeof(PetscInt), &crsGID ); CHKERRQ(ierr); 499 for(kk=0;kk<nSelected_1;kk++) crsGID[kk] = kk; 500 *a_Gmat_2 = 0; 501 *a_selected_2 = selected_1; /* needed? */ 502 } 503 else { 504 PetscInt idx,num_fine_ghosts,num_crs_ghost,nLocalSelected,myCrs0; 505 Mat_MPIAIJ *mpimat2; 506 Mat Gmat2; 507 Vec locState; 508 PetscScalar *cpcol_state; 509 510 /* get 'nLocalSelected' */ 511 ierr = ISGetIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 512 for(kk=0,nLocalSelected=0;kk<nSelected_1;kk++){ 513 PetscInt lid = selected_idx[kk]; 514 if(lid<nloc) nLocalSelected++; 515 } 516 ierr = ISRestoreIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 517 /* scan my coarse zero gid, set 'lid_state' with coarse GID */ 518 MPI_Scan( &nLocalSelected, &myCrs0, 1, MPIU_INT, MPIU_SUM, wcomm ); 519 myCrs0 -= nLocalSelected; 520 521 if( a_Gmat_2 != 0 ) { /* output */ 522 /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */ 523 ierr = MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2 ); CHKERRQ(ierr); 524 *a_Gmat_2 = Gmat2; /* output */ 525 } 526 else Gmat2 = Gmat1; /* use local to get crsGIDs at least */ 527 /* get coarse grid GIDS for selected (locals and ghosts) */ 528 mpimat2 = (Mat_MPIAIJ*)Gmat2->data; 529 ierr = MatGetVecs( Gmat2, &locState, 0 ); CHKERRQ(ierr); 530 ierr = VecSet( locState, (PetscScalar)(PetscReal)(-1) ); CHKERRQ(ierr); /* set with UNKNOWN state */ 531 ierr = ISGetIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 532 for(kk=0;kk<nLocalSelected;kk++){ 533 PetscInt fgid = selected_idx[kk] + my0; 534 PetscScalar v = (PetscScalar)(kk+myCrs0); 535 ierr = VecSetValues( locState, 1, &fgid, &v, INSERT_VALUES ); CHKERRQ(ierr); /* set with PID */ 536 } 537 ierr = ISRestoreIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 538 ierr = VecAssemblyBegin( locState ); CHKERRQ(ierr); 539 ierr = VecAssemblyEnd( locState ); CHKERRQ(ierr); 540 ierr = VecScatterBegin(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 541 ierr = VecScatterEnd(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 542 ierr = VecGetLocalSize( mpimat2->lvec, &num_fine_ghosts ); CHKERRQ(ierr); 543 ierr = VecGetArray( mpimat2->lvec, &cpcol_state ); CHKERRQ(ierr); 544 for(kk=0,num_crs_ghost=0;kk<num_fine_ghosts;kk++){ 545 if( (PetscInt)PetscRealPart(cpcol_state[kk]) != -1 ) num_crs_ghost++; 546 } 547 ierr = PetscMalloc( (nLocalSelected+num_crs_ghost)*sizeof(PetscInt), &crsGID ); CHKERRQ(ierr); /* output */ 548 { 549 PetscInt *selected_set; 550 ierr = PetscMalloc( (nLocalSelected+num_crs_ghost)*sizeof(PetscInt), &selected_set ); CHKERRQ(ierr); 551 /* do ghost of 'crsGID' */ 552 for(kk=0,idx=nLocalSelected;kk<num_fine_ghosts;kk++){ 553 if( (PetscInt)PetscRealPart(cpcol_state[kk]) != -1 ){ 554 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]); 555 selected_set[idx] = nloc + kk; 556 crsGID[idx++] = cgid; 557 } 558 } 559 assert(idx==(nLocalSelected+num_crs_ghost)); 560 ierr = VecRestoreArray( mpimat2->lvec, &cpcol_state ); CHKERRQ(ierr); 561 /* do locals in 'crsGID' */ 562 ierr = VecGetArray( locState, &cpcol_state ); CHKERRQ(ierr); 563 for(kk=0,idx=0;kk<nloc;kk++){ 564 if( (PetscInt)PetscRealPart(cpcol_state[kk]) != -1 ){ 565 PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]); 566 selected_set[idx] = kk; 567 crsGID[idx++] = cgid; 568 } 569 } 570 assert(idx==nLocalSelected); 571 ierr = VecRestoreArray( locState, &cpcol_state ); CHKERRQ(ierr); 572 573 if( a_selected_2 != 0 ) { /* output */ 574 ierr = ISCreateGeneral(PETSC_COMM_SELF,(nLocalSelected+num_crs_ghost),selected_set,PETSC_COPY_VALUES,a_selected_2); 575 CHKERRQ(ierr); 576 } 577 ierr = PetscFree( selected_set ); CHKERRQ(ierr); 578 } 579 ierr = VecDestroy( &locState ); CHKERRQ(ierr); 580 } 581 *a_crsGID = crsGID; /* output */ 582 583 PetscFunctionReturn(0); 584 } 585 586 /* -------------------------------------------------------------------------- */ 587 /* 588 PCGAMGcreateProl_GEO 589 590 Input Parameter: 591 . pc - this 592 . Amat - matrix on this fine level 593 . data[nloc*data_sz(in)] 594 Output Parameter: 595 . a_P_out - prolongation operator to the next level 596 . a_data_out - data of coarse grid points (num local columns in 'a_P_out') 597 */ 598 #undef __FUNCT__ 599 #define __FUNCT__ "PCGAMGcreateProl_GEO" 600 PetscErrorCode PCGAMGcreateProl_GEO( PC pc, 601 const Mat Amat, 602 const PetscReal data[], 603 Mat *a_P_out, 604 PetscReal **a_data_out 605 ) 606 { 607 PC_MG *mg = (PC_MG*)pc->data; 608 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 609 const PetscInt verbose = pc_gamg->verbose; 610 const PetscInt dim = pc_gamg->data_cols, data_cols = pc_gamg->data_cols; 611 PetscErrorCode ierr; 612 PetscInt Istart,Iend,nloc,jj,kk,my0,nLocalSelected,NN,MM,bs; 613 Mat Prol, Gmat; 614 PetscMPIInt mype, npe; 615 MPI_Comm wcomm = ((PetscObject)Amat)->comm; 616 IS permIS, llist_1, selected_1, selected_2; 617 const PetscInt *selected_idx; 618 619 PetscFunctionBegin; 620 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 621 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 622 ierr = MatGetOwnershipRange( Amat, &Istart, &Iend ); CHKERRQ(ierr); 623 ierr = MatGetSize( Amat, &MM, &NN ); CHKERRQ(ierr); 624 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ( ierr ); 625 nloc = (Iend-Istart)/bs; my0 = Istart/bs; assert((Iend-Istart)%bs==0); 626 627 ierr = createGraph( pc, Amat, &Gmat, PETSC_NULL, &permIS ); CHKERRQ(ierr); 628 629 /* SELECT COARSE POINTS */ 630 #if defined PETSC_USE_LOG 631 ierr = PetscLogEventBegin(gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 632 #endif 633 634 ierr = maxIndSetAgg( permIS, Gmat, PETSC_NULL, PETSC_FALSE, &selected_1, &llist_1 ); 635 CHKERRQ(ierr); 636 #if defined PETSC_USE_LOG 637 ierr = PetscLogEventEnd(gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 638 #endif 639 ierr = ISDestroy(&permIS); CHKERRQ(ierr); 640 641 /* get 'nLocalSelected' */ 642 ierr = ISGetLocalSize( selected_1, &NN ); CHKERRQ(ierr); 643 ierr = ISGetIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 644 for(kk=0,nLocalSelected=0;kk<NN;kk++){ 645 PetscInt lid = selected_idx[kk]; 646 if(lid<nloc) nLocalSelected++; 647 } 648 ierr = ISRestoreIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 649 650 /* create prolongator, create P matrix */ 651 ierr = MatCreateMPIAIJ(wcomm, 652 nloc*bs, nLocalSelected*bs, 653 PETSC_DETERMINE, PETSC_DETERMINE, 654 3*data_cols, PETSC_NULL, /* don't have a good way to set this!!! */ 655 3*data_cols, PETSC_NULL, 656 &Prol ); 657 CHKERRQ(ierr); 658 659 /* can get all points "removed" */ 660 ierr = MatGetSize( Prol, &kk, &jj ); CHKERRQ(ierr); 661 if( jj==0 ) { 662 if( verbose ) { 663 PetscPrintf(PETSC_COMM_WORLD,"[%d]%s no selected points on coarse grid\n",mype,__FUNCT__); 664 } 665 ierr = MatDestroy( &Prol ); CHKERRQ(ierr); 666 ierr = ISDestroy( &llist_1 ); CHKERRQ(ierr); 667 ierr = ISDestroy( &selected_1 ); CHKERRQ(ierr); 668 ierr = MatDestroy( &Gmat ); CHKERRQ(ierr); 669 *a_P_out = PETSC_NULL; /* out */ 670 PetscFunctionReturn(0); 671 } 672 673 { 674 PetscReal *coords; 675 PetscInt nnodes; 676 PetscInt *crsGID; 677 Mat Gmat2; 678 679 assert(dim==data_cols); 680 /* grow ghost data for better coarse grid cover of fine grid */ 681 #if defined PETSC_USE_LOG 682 ierr = PetscLogEventBegin(gamg_setup_events[SET5],0,0,0,0);CHKERRQ(ierr); 683 #endif 684 ierr = getGIDsOnSquareGraph( selected_1, Gmat, &selected_2, &Gmat2, &crsGID ); 685 CHKERRQ(ierr); 686 ierr = MatDestroy( &Gmat ); CHKERRQ(ierr); 687 /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */ 688 #if defined PETSC_USE_LOG 689 ierr = PetscLogEventEnd(gamg_setup_events[SET5],0,0,0,0);CHKERRQ(ierr); 690 #endif 691 /* create global vector of coorindates in 'coords' */ 692 if (npe > 1) { 693 ierr = getDataWithGhosts( Gmat2, dim, data, &nnodes, &coords ); 694 CHKERRQ(ierr); 695 } 696 else { 697 coords = (PetscReal*)data; 698 nnodes = nloc; 699 } 700 ierr = MatDestroy( &Gmat2 ); CHKERRQ(ierr); 701 702 /* triangulate */ 703 if( dim == 2 ) { 704 PetscReal metric=0.0; 705 #if defined PETSC_USE_LOG 706 ierr = PetscLogEventBegin(gamg_setup_events[SET6],0,0,0,0);CHKERRQ(ierr); 707 #endif 708 ierr = triangulateAndFormProl( selected_2, nnodes, coords, 709 selected_1, llist_1, crsGID, bs, Prol, &metric ); 710 CHKERRQ(ierr); 711 #if defined PETSC_USE_LOG 712 ierr = PetscLogEventEnd(gamg_setup_events[SET6],0,0,0,0); CHKERRQ(ierr); 713 #endif 714 ierr = PetscFree( crsGID ); CHKERRQ(ierr); 715 716 /* clean up and create coordinates for coarse grid (output) */ 717 if (npe > 1) ierr = PetscFree( coords ); CHKERRQ(ierr); 718 719 if( metric > 1. ) { /* needs to be globalized - should not happen */ 720 if( verbose ) { 721 PetscPrintf(PETSC_COMM_SELF,"[%d]%s failed metric for coarse grid %e\n",mype,__FUNCT__,metric); 722 } 723 ierr = MatDestroy( &Prol ); CHKERRQ(ierr); 724 Prol = PETSC_NULL; 725 } 726 else if( metric > .0 ) { 727 if( verbose ) { 728 PetscPrintf(PETSC_COMM_SELF,"[%d]%s metric for coarse grid = %e\n",mype,__FUNCT__,metric); 729 } 730 } 731 } else { 732 SETERRQ(wcomm,PETSC_ERR_LIB,"3D not implemented for 'geo' AMG"); 733 } 734 { /* create next coords - output */ 735 PetscReal *crs_crds; 736 ierr = PetscMalloc( dim*nLocalSelected*sizeof(PetscReal), &crs_crds ); 737 CHKERRQ(ierr); 738 ierr = ISGetIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 739 for(kk=0;kk<nLocalSelected;kk++){/* grab local select nodes to promote - output */ 740 PetscInt lid = selected_idx[kk]; 741 for(jj=0;jj<dim;jj++) crs_crds[jj*nLocalSelected + kk] = data[jj*nloc + lid]; 742 } 743 ierr = ISRestoreIndices( selected_1, &selected_idx ); CHKERRQ(ierr); 744 *a_data_out = crs_crds; /* out */ 745 } 746 if (npe > 1) { 747 ierr = ISDestroy( &selected_2 ); CHKERRQ(ierr); /* this is selected_1 in serial */ 748 } 749 } 750 *a_P_out = Prol; /* out */ 751 752 ierr = ISDestroy( &llist_1 ); CHKERRQ(ierr); 753 ierr = ISDestroy( &selected_1 ); CHKERRQ(ierr); 754 755 PetscFunctionReturn(0); 756 } 757