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 #include <assert.h> 9 #include <petscblaslapack.h> 10 11 typedef struct { 12 PetscInt nsmooths; 13 PetscBool sym_graph; 14 PetscBool square_graph; 15 }PC_GAMG_AGG; 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "PCGAMGSetNSmooths" 19 /*@ 20 PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) 21 22 Not Collective on PC 23 24 Input Parameters: 25 . pc - the preconditioner context 26 27 Options Database Key: 28 . -pc_gamg_agg_nsmooths 29 30 Level: intermediate 31 32 Concepts: Aggregation AMG preconditioner 33 34 .seealso: () 35 @*/ 36 PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n) 37 { 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 42 ierr = PetscTryMethod(pc,"PCGAMGSetNSmooths_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 EXTERN_C_BEGIN 47 #undef __FUNCT__ 48 #define __FUNCT__ "PCGAMGSetNSmooths_GAMG" 49 PetscErrorCode PCGAMGSetNSmooths_GAMG(PC pc, PetscInt n) 50 { 51 PC_MG *mg = (PC_MG*)pc->data; 52 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 53 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 54 55 PetscFunctionBegin; 56 pc_gamg_agg->nsmooths = n; 57 PetscFunctionReturn(0); 58 } 59 EXTERN_C_END 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "PCGAMGSetSymGraph" 63 /*@ 64 PCGAMGSetSymGraph - 65 66 Not Collective on PC 67 68 Input Parameters: 69 . pc - the preconditioner context 70 71 Options Database Key: 72 . -pc_gamg_sym_graph 73 74 Level: intermediate 75 76 Concepts: Aggregation AMG preconditioner 77 78 .seealso: () 79 @*/ 80 PetscErrorCode PCGAMGSetSymGraph(PC pc, PetscBool n) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 86 ierr = PetscTryMethod(pc,"PCGAMGSetSymGraph_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 EXTERN_C_BEGIN 91 #undef __FUNCT__ 92 #define __FUNCT__ "PCGAMGSetSymGraph_GAMG" 93 PetscErrorCode PCGAMGSetSymGraph_GAMG(PC pc, PetscBool n) 94 { 95 PC_MG *mg = (PC_MG*)pc->data; 96 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 97 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 98 99 PetscFunctionBegin; 100 pc_gamg_agg->sym_graph = n; 101 PetscFunctionReturn(0); 102 } 103 EXTERN_C_END 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "PCGAMGSetSquareGraph" 107 /*@ 108 PCGAMGSetSquareGraph - 109 110 Not Collective on PC 111 112 Input Parameters: 113 . pc - the preconditioner context 114 115 Options Database Key: 116 . -pc_gamg_square_graph 117 118 Level: intermediate 119 120 Concepts: Aggregation AMG preconditioner 121 122 .seealso: () 123 @*/ 124 PetscErrorCode PCGAMGSetSquareGraph(PC pc, PetscBool n) 125 { 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 130 ierr = PetscTryMethod(pc,"PCGAMGSetSquareGraph_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 134 EXTERN_C_BEGIN 135 #undef __FUNCT__ 136 #define __FUNCT__ "PCGAMGSetSquareGraph_GAMG" 137 PetscErrorCode PCGAMGSetSquareGraph_GAMG(PC pc, PetscBool n) 138 { 139 PC_MG *mg = (PC_MG*)pc->data; 140 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 141 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 142 143 PetscFunctionBegin; 144 pc_gamg_agg->square_graph = n; 145 PetscFunctionReturn(0); 146 } 147 EXTERN_C_END 148 149 /* -------------------------------------------------------------------------- */ 150 /* 151 PCSetFromOptions_GAMG_AGG 152 153 Input Parameter: 154 . pc - 155 */ 156 #undef __FUNCT__ 157 #define __FUNCT__ "PCSetFromOptions_GAMG_AGG" 158 PetscErrorCode PCSetFromOptions_GAMG_AGG( PC pc ) 159 { 160 PetscErrorCode ierr; 161 PC_MG *mg = (PC_MG*)pc->data; 162 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 163 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 164 PetscBool flag; 165 166 PetscFunctionBegin; 167 /* call base class */ 168 ierr = PCSetFromOptions_GAMG( pc ); CHKERRQ(ierr); 169 170 ierr = PetscOptionsHead("GAMG-AGG options"); CHKERRQ(ierr); 171 { 172 /* -pc_gamg_agg_nsmooths */ 173 pc_gamg_agg->nsmooths = 0; 174 ierr = PetscOptionsInt("-pc_gamg_agg_nsmooths", 175 "smoothing steps for smoothed aggregation, usually 1 (0)", 176 "PCGAMGSetNSmooths", 177 pc_gamg_agg->nsmooths, 178 &pc_gamg_agg->nsmooths, 179 &flag); 180 CHKERRQ(ierr); 181 182 /* -pc_gamg_sym_graph */ 183 pc_gamg_agg->sym_graph = PETSC_FALSE; 184 ierr = PetscOptionsBool("-pc_gamg_sym_graph", 185 "Set for asymmetric matrices", 186 "PCGAMGSetSymGraph", 187 pc_gamg_agg->sym_graph, 188 &pc_gamg_agg->sym_graph, 189 &flag); 190 CHKERRQ(ierr); 191 192 /* -pc_gamg_square_graph */ 193 pc_gamg_agg->square_graph = PETSC_TRUE; 194 ierr = PetscOptionsBool("-pc_gamg_square_graph", 195 "For faster coarsening and lower coarse grid complexity", 196 "PCGAMGSetSquareGraph", 197 pc_gamg_agg->square_graph, 198 &pc_gamg_agg->square_graph, 199 &flag); 200 CHKERRQ(ierr); 201 } 202 ierr = PetscOptionsTail();CHKERRQ(ierr); 203 204 PetscFunctionReturn(0); 205 } 206 207 /* -------------------------------------------------------------------------- */ 208 /* 209 PCDestroy_AGG 210 211 Input Parameter: 212 . pc - 213 */ 214 #undef __FUNCT__ 215 #define __FUNCT__ "PCDestroy_AGG" 216 PetscErrorCode PCDestroy_AGG( PC pc ) 217 { 218 PetscErrorCode ierr; 219 PC_MG *mg = (PC_MG*)pc->data; 220 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 221 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 222 223 PetscFunctionBegin; 224 if( pc_gamg_agg ) { 225 ierr = PetscFree(pc_gamg_agg);CHKERRQ(ierr); 226 pc_gamg_agg = 0; 227 } 228 229 /* call base class */ 230 ierr = PCDestroy_GAMG( pc );CHKERRQ(ierr); 231 232 PetscFunctionReturn(0); 233 } 234 235 /* -------------------------------------------------------------------------- */ 236 /* 237 PCSetCoordinates_AGG 238 - collective 239 240 Input Parameter: 241 . pc - the preconditioner context 242 . ndm - dimesion 243 . a_nloc - number of vertices local 244 . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...} 245 */ 246 EXTERN_C_BEGIN 247 #undef __FUNCT__ 248 #define __FUNCT__ "PCSetCoordinates_AGG" 249 PetscErrorCode PCSetCoordinates_AGG( PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords ) 250 { 251 PC_MG *mg = (PC_MG*)pc->data; 252 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 253 PetscErrorCode ierr; 254 PetscInt arrsz,bs,kk,ii,jj,nloc; 255 Mat Amat = pc->pmat; 256 MPI_Comm wcomm = ((PetscObject)pc)->comm; 257 258 PetscFunctionBegin; 259 PetscValidHeaderSpecific( Amat, MAT_CLASSID, 1 ); 260 /* set 'bs' and 'nloc' */ 261 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ( ierr ); 262 if( a_nloc == -1 ) { 263 PetscInt my0, Iend; 264 /* stokes = PETSC_FALSE; */ 265 ierr = MatGetOwnershipRange( Amat, &my0, &Iend ); CHKERRQ(ierr); 266 nloc = (Iend-my0)/bs; 267 if((Iend-my0)%bs!=0) SETERRQ1(((PetscObject)Amat)->comm,PETSC_ERR_ARG_WRONG, "Bad local size %d.",nloc); 268 } 269 else { 270 ierr = MPI_Allreduce( &a_nloc, &ii, 1, MPIU_INT, MPIU_SUM, wcomm ); CHKERRQ( ierr ); 271 ierr = MatGetSize( Amat, &kk, &jj ); CHKERRQ( ierr ); 272 if( bs==1 && ii!=kk ) { 273 /* stokes = PETSC_TRUE; */ 274 bs = ndm; 275 nloc = a_nloc; 276 } 277 else{ 278 PetscInt my0,Iend; 279 /* stokes = PETSC_FALSE; */ 280 ierr = MatGetOwnershipRange( Amat, &my0, &Iend ); CHKERRQ(ierr); 281 nloc = (Iend-my0)/bs; 282 if((Iend-my0)%bs!=0) SETERRQ1(((PetscObject)Amat)->comm,PETSC_ERR_ARG_WRONG, "Bad local size %d.",nloc); 283 } 284 } 285 286 /* SA: null space vectors */ 287 if( coords && bs==1 ) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */ 288 else if( coords ) pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* elasticity */ 289 else pc_gamg->data_cell_cols = bs; /* no data, force SA with constant null space vectors */ 290 pc_gamg->data_cell_rows = bs; 291 292 arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols; 293 294 /* create data - syntactic sugar that should be refactored at some point */ 295 if (pc_gamg->data==0 || (pc_gamg->data_sz != arrsz)) { 296 ierr = PetscFree( pc_gamg->data ); CHKERRQ(ierr); 297 ierr = PetscMalloc((arrsz+1)*sizeof(PetscReal), &pc_gamg->data ); CHKERRQ(ierr); 298 /* !nul if if nloc==0 */ 299 } 300 /* copy data in - column oriented */ 301 for(kk=0;kk<nloc;kk++){ 302 const PetscInt M = nloc*pc_gamg->data_cell_rows; 303 PetscReal *data = &pc_gamg->data[kk*bs]; 304 if( pc_gamg->data_cell_cols==1 ) *data = 1.0; 305 else { 306 for(ii=0;ii<bs;ii++) 307 for(jj=0;jj<bs;jj++) 308 if(ii==jj)data[ii*M + jj] = 1.0; /* translational modes */ 309 else data[ii*M + jj] = 0.0; 310 if( coords ) { 311 if( ndm == 2 ){ /* rotational modes */ 312 data += 2*M; 313 data[0] = -coords[2*kk+1]; 314 data[1] = coords[2*kk]; 315 } 316 else { 317 data += 3*M; 318 data[0] = 0.0; data[M+0] = coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1]; 319 data[1] = -coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = coords[3*kk]; 320 data[2] = coords[3*kk+1]; data[M+2] = -coords[3*kk]; data[2*M+2] = 0.0; 321 } 322 } 323 } 324 } 325 326 pc_gamg->data_sz = arrsz; 327 328 PetscFunctionReturn(0); 329 } 330 EXTERN_C_END 331 332 typedef PetscInt NState; 333 static const NState NOT_DONE=-2; 334 static const NState DELETED=-1; 335 static const NState REMOVED=-3; 336 #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED) 337 338 /* -------------------------------------------------------------------------- */ 339 /* 340 smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific 341 - AGG-MG specific: clears singletons out of 'selected_2' 342 343 Input Parameter: 344 . Gmat_2 - glabal matrix of graph (data not defined) 345 . Gmat_1 - base graph to grab with 346 Input/Output Parameter: 347 . aggs_2 - linked list of aggs with gids ) 348 */ 349 #undef __FUNCT__ 350 #define __FUNCT__ "smoothAggs" 351 static PetscErrorCode smoothAggs( const Mat Gmat_2, /* base (squared) graph */ 352 const Mat Gmat_1, /* base graph */ 353 /* const IS selected_2, [nselected local] selected vertices */ 354 PetscCoarsenData *aggs_2 /* [nselected local] global ID of aggregate */ 355 ) 356 { 357 PetscErrorCode ierr; 358 PetscBool isMPI; 359 Mat_SeqAIJ *matA_1, *matB_1=0, *matA_2, *matB_2=0; 360 MPI_Comm wcomm = ((PetscObject)Gmat_2)->comm; 361 PetscMPIInt mype,npe; 362 PetscInt lid,*ii,*idx,ix,Iend,my0,kk,n,j; 363 Mat_MPIAIJ *mpimat_2 = 0, *mpimat_1=0; 364 const PetscInt nloc = Gmat_2->rmap->n; 365 PetscScalar *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid; 366 PetscInt *lid_cprowID_1; 367 NState *lid_state; 368 Vec ghost_par_orig2; 369 370 PetscFunctionBegin; 371 ierr = MPI_Comm_rank( wcomm, &mype ); CHKERRQ(ierr); 372 ierr = MPI_Comm_size( wcomm, &npe ); CHKERRQ(ierr); 373 ierr = MatGetOwnershipRange(Gmat_1,&my0,&Iend); CHKERRQ(ierr); 374 375 if( PETSC_FALSE ) { 376 PetscViewer viewer; char fname[32]; static int llev=0; 377 sprintf(fname,"Gmat2_%d.m",llev++); 378 PetscViewerASCIIOpen(wcomm,fname,&viewer); 379 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 380 ierr = MatView(Gmat_2, viewer ); CHKERRQ(ierr); 381 ierr = PetscViewerDestroy( &viewer ); 382 } 383 384 /* get submatrices */ 385 ierr = PetscTypeCompare( (PetscObject)Gmat_1, MATMPIAIJ, &isMPI ); CHKERRQ(ierr); 386 if(isMPI) { 387 /* grab matrix objects */ 388 mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data; 389 mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data; 390 matA_1 = (Mat_SeqAIJ*)mpimat_1->A->data; 391 matB_1 = (Mat_SeqAIJ*)mpimat_1->B->data; 392 matA_2 = (Mat_SeqAIJ*)mpimat_2->A->data; 393 matB_2 = (Mat_SeqAIJ*)mpimat_2->B->data; 394 395 /* force compressed row storage for B matrix in AuxMat */ 396 matB_1->compressedrow.check = PETSC_TRUE; 397 ierr = MatCheckCompressedRow(mpimat_1->B,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0); 398 CHKERRQ(ierr); 399 400 ierr = PetscMalloc( nloc*sizeof(PetscInt), &lid_cprowID_1 ); CHKERRQ(ierr); 401 for( lid = 0 ; lid < nloc ; lid++ ) lid_cprowID_1[lid] = -1; 402 for (ix=0; ix<matB_1->compressedrow.nrows; ix++) { 403 PetscInt lid = matB_1->compressedrow.rindex[ix]; 404 lid_cprowID_1[lid] = ix; 405 } 406 } 407 else { 408 matA_1 = (Mat_SeqAIJ*)Gmat_1->data; 409 matA_2 = (Mat_SeqAIJ*)Gmat_2->data; 410 lid_cprowID_1 = PETSC_NULL; 411 } 412 assert( matA_1 && !matA_1->compressedrow.use ); 413 assert( matB_1==0 || matB_1->compressedrow.use ); 414 assert( matA_2 && !matA_2->compressedrow.use ); 415 assert( matB_2==0 || matB_2->compressedrow.use ); 416 417 /* get state of locals and selected gid for deleted */ 418 ierr = PetscMalloc( nloc*sizeof(NState), &lid_state ); CHKERRQ(ierr); 419 ierr = PetscMalloc( nloc*sizeof(PetscScalar), &lid_parent_gid ); CHKERRQ(ierr); 420 for( lid = 0 ; lid < nloc ; lid++ ) { 421 lid_parent_gid[lid] = -1.0; 422 lid_state[lid] = DELETED; 423 } 424 425 /* set lid_state */ 426 for( lid = 0 ; lid < nloc ; lid++ ) { 427 PetscCDPos pos; 428 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos); CHKERRQ(ierr); 429 if( pos ) { 430 PetscInt gid1; 431 ierr = PetscLLNGetID( pos, &gid1 ); CHKERRQ(ierr); assert(gid1==lid+my0); 432 lid_state[lid] = gid1; 433 } 434 } 435 436 /* map local to selected local, DELETED means a ghost owns it */ 437 for(lid=kk=0;lid<nloc;lid++){ 438 NState state = lid_state[lid]; 439 if( IS_SELECTED(state) ){ 440 PetscCDPos pos; 441 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos); CHKERRQ(ierr); 442 while(pos){ 443 PetscInt gid1; 444 ierr = PetscLLNGetID( pos, &gid1 ); CHKERRQ(ierr); 445 ierr = PetscCDGetNextPos(aggs_2,lid,&pos); CHKERRQ(ierr); 446 447 if( gid1 >= my0 && gid1 < Iend ){ 448 lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0); 449 } 450 } 451 } 452 } 453 /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */ 454 if (isMPI) { 455 Vec tempVec; 456 /* get 'cpcol_1_state' */ 457 ierr = MatGetVecs( Gmat_1, &tempVec, 0 ); CHKERRQ(ierr); 458 for(kk=0,j=my0;kk<nloc;kk++,j++){ 459 PetscScalar v = (PetscScalar)lid_state[kk]; 460 ierr = VecSetValues( tempVec, 1, &j, &v, INSERT_VALUES ); CHKERRQ(ierr); 461 } 462 ierr = VecAssemblyBegin( tempVec ); CHKERRQ(ierr); 463 ierr = VecAssemblyEnd( tempVec ); CHKERRQ(ierr); 464 ierr = VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD); 465 CHKERRQ(ierr); 466 ierr = VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD); 467 CHKERRQ(ierr); 468 ierr = VecGetArray( mpimat_1->lvec, &cpcol_1_state ); CHKERRQ(ierr); 469 /* get 'cpcol_2_state' */ 470 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD); 471 CHKERRQ(ierr); 472 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD); 473 CHKERRQ(ierr); 474 ierr = VecGetArray( mpimat_2->lvec, &cpcol_2_state ); CHKERRQ(ierr); 475 /* get 'cpcol_2_par_orig' */ 476 for(kk=0,j=my0;kk<nloc;kk++,j++){ 477 PetscScalar v = (PetscScalar)lid_parent_gid[kk]; 478 ierr = VecSetValues( tempVec, 1, &j, &v, INSERT_VALUES ); CHKERRQ(ierr); 479 } 480 ierr = VecAssemblyBegin( tempVec ); CHKERRQ(ierr); 481 ierr = VecAssemblyEnd( tempVec ); CHKERRQ(ierr); 482 ierr = VecDuplicate( mpimat_2->lvec, &ghost_par_orig2 ); CHKERRQ(ierr); 483 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD); 484 CHKERRQ(ierr); 485 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD); 486 CHKERRQ(ierr); 487 ierr = VecGetArray( ghost_par_orig2, &cpcol_2_par_orig ); CHKERRQ(ierr); 488 489 ierr = VecDestroy( &tempVec ); CHKERRQ(ierr); 490 } /* ismpi */ 491 492 /* doit */ 493 for(lid=0;lid<nloc;lid++){ 494 NState state = lid_state[lid]; 495 if( IS_SELECTED(state) ) { 496 /* steal locals */ 497 ii = matA_1->i; n = ii[lid+1] - ii[lid]; 498 idx = matA_1->j + ii[lid]; 499 for (j=0; j<n; j++) { 500 PetscInt lidj = idx[j], sgid; 501 NState statej = lid_state[lidj]; 502 if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */ 503 lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */ 504 if( sgid >= my0 && sgid < Iend ){ /* I'm stealing this local from a local sgid */ 505 PetscInt hav=0,slid=sgid-my0,gidj=lidj+my0; 506 PetscCDPos pos,last=PETSC_NULL; 507 /* looking for local from local so id_llist_2 works */ 508 ierr = PetscCDGetHeadPos(aggs_2,slid,&pos); CHKERRQ(ierr); 509 while(pos){ 510 PetscInt gid; 511 ierr = PetscLLNGetID( pos, &gid ); CHKERRQ(ierr); 512 if( gid == gidj ) { 513 assert(last); 514 ierr = PetscCDRemoveNextNode( aggs_2, slid, last ); CHKERRQ(ierr); 515 ierr = PetscCDAppendNode( aggs_2, lid, pos ); CHKERRQ(ierr); 516 hav = 1; 517 break; 518 } 519 else last = pos; 520 521 ierr = PetscCDGetNextPos(aggs_2,slid,&pos); CHKERRQ(ierr); 522 } 523 if(hav!=1){ 524 if(hav==0)SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix"); 525 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"found node %d times???",hav); 526 } 527 } 528 else{ /* I'm stealing this local, owned by a ghost */ 529 assert(sgid==-1); 530 ierr = PetscCDAppendID( aggs_2, lid, lidj+my0 ); CHKERRQ(ierr); 531 } 532 } 533 } /* local neighbors */ 534 } 535 else if( state == DELETED && lid_cprowID_1 ) { 536 PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]); 537 /* see if I have a selected ghost neighbor that will steal me */ 538 if( (ix=lid_cprowID_1[lid]) != -1 ){ 539 ii = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix]; 540 idx = matB_1->j + ii[ix]; 541 for( j=0 ; j<n ; j++ ) { 542 PetscInt cpid = idx[j]; 543 NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]); 544 if( IS_SELECTED(statej) && sgidold != (PetscInt)statej ) { /* ghost will steal this, remove from my list */ 545 lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */ 546 if( sgidold>=my0 && sgidold<Iend ) { /* this was mine */ 547 PetscInt hav=0,oldslidj=sgidold-my0; 548 PetscCDPos pos,last=PETSC_NULL; 549 /* remove from 'oldslidj' list */ 550 ierr = PetscCDGetHeadPos(aggs_2,oldslidj,&pos); CHKERRQ(ierr); 551 while( pos ) { 552 PetscInt gid; 553 ierr = PetscLLNGetID( pos, &gid ); CHKERRQ(ierr); 554 if( lid+my0 == gid ) { 555 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */ 556 assert(last); 557 ierr = PetscCDRemoveNextNode( aggs_2, oldslidj, last ); CHKERRQ(ierr); 558 /* ghost (PetscScalar)statej will add this later */ 559 hav = 1; 560 break; 561 } 562 else last = pos; 563 564 ierr = PetscCDGetNextPos(aggs_2,oldslidj,&pos); CHKERRQ(ierr); 565 } 566 if(hav!=1){ 567 if(hav==0)SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix"); 568 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"found node %d times???",hav); 569 } 570 } 571 else { 572 /* ghosts remove this later */ 573 } 574 } 575 } 576 } 577 } /* selected/deleted */ 578 } /* node loop */ 579 580 if( isMPI ) { 581 PetscScalar *cpcol_2_parent,*cpcol_2_gid; 582 Vec tempVec,ghostgids2,ghostparents2; 583 PetscInt cpid,nghost_2; 584 GAMGHashTable gid_cpid; 585 586 ierr = VecGetSize( mpimat_2->lvec, &nghost_2 ); CHKERRQ(ierr); 587 ierr = MatGetVecs( Gmat_2, &tempVec, 0 ); CHKERRQ(ierr); 588 589 /* get 'cpcol_2_parent' */ 590 for(kk=0,j=my0;kk<nloc;kk++,j++){ 591 ierr = VecSetValues( tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES ); CHKERRQ(ierr); 592 } 593 ierr = VecAssemblyBegin( tempVec ); CHKERRQ(ierr); 594 ierr = VecAssemblyEnd( tempVec ); CHKERRQ(ierr); 595 ierr = VecDuplicate( mpimat_2->lvec, &ghostparents2 ); CHKERRQ(ierr); 596 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD); 597 CHKERRQ(ierr); 598 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD); 599 CHKERRQ(ierr); 600 ierr = VecGetArray( ghostparents2, &cpcol_2_parent ); CHKERRQ(ierr); 601 602 /* get 'cpcol_2_gid' */ 603 for(kk=0,j=my0;kk<nloc;kk++,j++){ 604 PetscScalar v = (PetscScalar)j; 605 ierr = VecSetValues( tempVec, 1, &j, &v, INSERT_VALUES ); CHKERRQ(ierr); 606 } 607 ierr = VecAssemblyBegin( tempVec ); CHKERRQ(ierr); 608 ierr = VecAssemblyEnd( tempVec ); CHKERRQ(ierr); 609 ierr = VecDuplicate( mpimat_2->lvec, &ghostgids2 ); CHKERRQ(ierr); 610 ierr = VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD); 611 CHKERRQ(ierr); 612 ierr = VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD); 613 CHKERRQ(ierr); 614 ierr = VecGetArray( ghostgids2, &cpcol_2_gid ); CHKERRQ(ierr); 615 616 ierr = VecDestroy( &tempVec ); CHKERRQ(ierr); 617 618 /* look for deleted ghosts and add to table */ 619 ierr = GAMGTableCreate( 2*nghost_2, &gid_cpid ); CHKERRQ(ierr); 620 for( cpid = 0 ; cpid < nghost_2 ; cpid++ ) { 621 NState state = (NState)PetscRealPart(cpcol_2_state[cpid]); 622 if( state==DELETED ) { 623 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 624 PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]); 625 if( sgid_old == -1 && sgid_new != -1 ) { 626 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 627 ierr = GAMGTableAdd( &gid_cpid, gid, cpid ); CHKERRQ(ierr); 628 } 629 } 630 } 631 632 /* look for deleted ghosts and see if they moved - remove it */ 633 for(lid=0;lid<nloc;lid++){ 634 NState state = lid_state[lid]; 635 if( IS_SELECTED(state) ){ 636 PetscCDPos pos,last=PETSC_NULL; 637 /* look for deleted ghosts and see if they moved */ 638 ierr = PetscCDGetHeadPos(aggs_2,lid,&pos); CHKERRQ(ierr); 639 while(pos){ 640 PetscInt gid; 641 ierr = PetscLLNGetID( pos, &gid ); CHKERRQ(ierr); 642 643 if( gid < my0 || gid >= Iend ) { 644 ierr = GAMGTableFind( &gid_cpid, gid, &cpid ); CHKERRQ(ierr); 645 if( cpid != -1 ) { 646 /* a moved ghost - */ 647 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */ 648 ierr = PetscCDRemoveNextNode( aggs_2, lid, last ); CHKERRQ(ierr); 649 } 650 else last = pos; 651 } 652 else last = pos; 653 654 ierr = PetscCDGetNextPos(aggs_2,lid,&pos); CHKERRQ(ierr); 655 } /* loop over list of deleted */ 656 } /* selected */ 657 } 658 ierr = GAMGTableDestroy( &gid_cpid ); CHKERRQ(ierr); 659 660 /* look at ghosts, see if they changed - and it */ 661 for( cpid = 0 ; cpid < nghost_2 ; cpid++ ){ 662 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 663 if( sgid_new >= my0 && sgid_new < Iend ) { /* this is mine */ 664 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 665 PetscInt slid_new=sgid_new-my0,hav=0; 666 PetscCDPos pos; 667 /* search for this gid to see if I have it */ 668 ierr = PetscCDGetHeadPos(aggs_2,slid_new,&pos); CHKERRQ(ierr); 669 while(pos){ 670 PetscInt gidj; 671 ierr = PetscLLNGetID( pos, &gidj ); CHKERRQ(ierr); 672 ierr = PetscCDGetNextPos(aggs_2,slid_new,&pos); CHKERRQ(ierr); 673 674 if( gidj == gid ) { hav = 1; break; } 675 } 676 if( hav != 1 ){ 677 /* insert 'flidj' into head of llist */ 678 ierr = PetscCDAppendID( aggs_2, slid_new, gid ); CHKERRQ(ierr); 679 } 680 } 681 } 682 683 ierr = VecRestoreArray( mpimat_1->lvec, &cpcol_1_state ); CHKERRQ(ierr); 684 ierr = VecRestoreArray( mpimat_2->lvec, &cpcol_2_state ); CHKERRQ(ierr); 685 ierr = VecRestoreArray( ghostparents2, &cpcol_2_parent ); CHKERRQ(ierr); 686 ierr = VecRestoreArray( ghostgids2, &cpcol_2_gid ); CHKERRQ(ierr); 687 ierr = PetscFree( lid_cprowID_1 ); CHKERRQ(ierr); 688 ierr = VecDestroy( &ghostgids2 ); CHKERRQ(ierr); 689 ierr = VecDestroy( &ghostparents2 ); CHKERRQ(ierr); 690 ierr = VecDestroy( &ghost_par_orig2 ); CHKERRQ(ierr); 691 } 692 693 ierr = PetscFree( lid_parent_gid ); CHKERRQ(ierr); 694 ierr = PetscFree( lid_state ); CHKERRQ(ierr); 695 696 PetscFunctionReturn(0); 697 } 698 699 /* -------------------------------------------------------------------------- */ 700 /* 701 PCSetData_AGG 702 703 Input Parameter: 704 . pc - 705 */ 706 #undef __FUNCT__ 707 #define __FUNCT__ "PCSetData_AGG" 708 PetscErrorCode PCSetData_AGG( PC pc, Mat a_A ) 709 { 710 PetscErrorCode ierr; 711 PC_MG *mg = (PC_MG*)pc->data; 712 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 713 MatNullSpace mnull; 714 715 PetscFunctionBegin; 716 ierr = MatGetNearNullSpace( a_A, &mnull ); CHKERRQ(ierr); 717 if( !mnull ) { 718 ierr = PCSetCoordinates_AGG( pc, -1, -1, PETSC_NULL ); CHKERRQ(ierr); 719 } 720 else { 721 PetscReal *nullvec; 722 PetscBool has_const; 723 PetscInt i,j,mlocal,nvec,bs; 724 const Vec *vecs; const PetscScalar *v; 725 ierr = MatGetLocalSize(a_A,&mlocal,PETSC_NULL);CHKERRQ(ierr); 726 ierr = MatNullSpaceGetVecs(mnull,&has_const,&nvec,&vecs);CHKERRQ(ierr); 727 ierr = MatGetBlockSize( a_A, &bs ); CHKERRQ( ierr ); 728 ierr = PetscMalloc((nvec+!!has_const)*mlocal*sizeof *nullvec,&nullvec);CHKERRQ(ierr); 729 if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0; 730 for (i=0; i<nvec; i++) { 731 ierr = VecGetArrayRead(vecs[i],&v);CHKERRQ(ierr); 732 for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]); 733 ierr = VecRestoreArrayRead(vecs[i],&v);CHKERRQ(ierr); 734 } 735 pc_gamg->data = nullvec; 736 pc_gamg->data_cell_cols = (nvec+!!has_const); 737 pc_gamg->data_cell_rows = bs; 738 } 739 PetscFunctionReturn(0); 740 } 741 742 /* -------------------------------------------------------------------------- */ 743 /* 744 formProl0 745 746 Input Parameter: 747 . agg_llists - list of arrays with aggregates 748 . bs - block size 749 . nSAvec - column bs of new P 750 . my0crs - global index of start of locals 751 . data_stride - bs*(nloc nodes + ghost nodes) 752 . data_in[data_stride*nSAvec] - local data on fine grid 753 . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist' 754 Output Parameter: 755 . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data 756 . a_Prol - prolongation operator 757 */ 758 #undef __FUNCT__ 759 #define __FUNCT__ "formProl0" 760 static PetscErrorCode formProl0(const PetscCoarsenData *agg_llists,/* list from selected vertices of aggregate unselected vertices */ 761 const PetscInt bs, /* (row) block size */ 762 const PetscInt nSAvec, /* column bs */ 763 const PetscInt my0crs, /* global index of start of locals */ 764 const PetscInt data_stride, /* (nloc+nghost)*bs */ 765 PetscReal data_in[], /* [data_stride][nSAvec] */ 766 const PetscInt flid_fgid[], /* [data_stride/bs] */ 767 PetscReal **a_data_out, 768 Mat a_Prol /* prolongation operator (output)*/ 769 ) 770 { 771 PetscErrorCode ierr; 772 PetscInt Istart,my0,Iend,nloc,clid,flid,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride; 773 MPI_Comm wcomm = ((PetscObject)a_Prol)->comm; 774 PetscMPIInt mype, npe; 775 PetscReal *out_data; 776 PetscCDPos pos; 777 GAMGHashTable fgid_flid; 778 779 /* #define OUT_AGGS */ 780 #ifdef OUT_AGGS 781 static PetscInt llev = 0; char fname[32]; FILE *file = PETSC_NULL; PetscInt pM; 782 #endif 783 784 PetscFunctionBegin; 785 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 786 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 787 ierr = MatGetOwnershipRange( a_Prol, &Istart, &Iend ); CHKERRQ(ierr); 788 nloc = (Iend-Istart)/bs; my0 = Istart/bs; assert((Iend-Istart)%bs==0); 789 Iend /= bs; 790 nghosts = data_stride/bs - nloc; 791 792 ierr = GAMGTableCreate( 2*nghosts, &fgid_flid ); CHKERRQ(ierr); 793 for(kk=0;kk<nghosts;kk++) { 794 ierr = GAMGTableAdd( &fgid_flid, flid_fgid[nloc+kk], nloc+kk ); CHKERRQ(ierr); 795 } 796 797 #ifdef OUT_AGGS 798 sprintf(fname,"aggs_%d_%d.m",llev++,mype); 799 if(llev==1) { 800 file = fopen(fname,"w"); 801 } 802 MatGetSize( a_Prol, &pM, &jj ); 803 #endif 804 805 /* count selected -- same as number of cols of P */ 806 for(nSelected=mm=0;mm<nloc;mm++) { 807 PetscBool ise; 808 ierr = PetscCDEmptyAt( agg_llists, mm, &ise ); CHKERRQ(ierr); 809 if( !ise ) nSelected++; 810 } 811 ierr = MatGetOwnershipRangeColumn( a_Prol, &ii, &jj ); CHKERRQ(ierr); 812 assert((ii/nSAvec)==my0crs); assert(nSelected==(jj-ii)/nSAvec); 813 814 /* aloc space for coarse point data (output) */ 815 out_data_stride = nSelected*nSAvec; 816 ierr = PetscMalloc( out_data_stride*nSAvec*sizeof(PetscReal), &out_data ); CHKERRQ(ierr); 817 for(ii=0;ii<out_data_stride*nSAvec;ii++) { 818 out_data[ii]=1.e300; 819 } 820 *a_data_out = out_data; /* output - stride nSelected*nSAvec */ 821 822 /* find points and set prolongation */ 823 minsz = 100; 824 ndone = 0; 825 for( mm = clid = 0 ; mm < nloc ; mm++ ){ 826 ierr = PetscCDSizeAt( agg_llists, mm, &jj ); CHKERRQ(ierr); 827 if( jj > 0 ) { 828 const PetscInt lid = mm, cgid = my0crs + clid; 829 PetscInt cids[100]; /* max bs */ 830 PetscBLASInt asz=jj,M=asz*bs,N=nSAvec,INFO; 831 PetscBLASInt Mdata=M+((N-M>0)?N-M:0),LDA=Mdata,LWORK=N*bs; 832 PetscScalar *qqc,*qqr,*TAU,*WORK; 833 PetscInt *fids; 834 PetscReal *data; 835 /* count agg */ 836 if( asz<minsz ) minsz = asz; 837 838 /* get block */ 839 ierr = PetscMalloc( (Mdata*N)*sizeof(PetscScalar), &qqc ); CHKERRQ(ierr); 840 ierr = PetscMalloc( (M*N)*sizeof(PetscScalar), &qqr ); CHKERRQ(ierr); 841 ierr = PetscMalloc( N*sizeof(PetscScalar), &TAU ); CHKERRQ(ierr); 842 ierr = PetscMalloc( LWORK*sizeof(PetscScalar), &WORK ); CHKERRQ(ierr); 843 ierr = PetscMalloc( M*sizeof(PetscInt), &fids ); CHKERRQ(ierr); 844 845 aggID = 0; 846 ierr = PetscCDGetHeadPos(agg_llists,lid,&pos); CHKERRQ(ierr); 847 while(pos){ 848 PetscInt gid1; 849 ierr = PetscLLNGetID( pos, &gid1 ); CHKERRQ(ierr); 850 ierr = PetscCDGetNextPos(agg_llists,lid,&pos); CHKERRQ(ierr); 851 852 if( gid1 >= my0 && gid1 < Iend ) flid = gid1 - my0; 853 else { 854 ierr = GAMGTableFind( &fgid_flid, gid1, &flid ); CHKERRQ(ierr); 855 assert(flid>=0); 856 } 857 /* copy in B_i matrix - column oriented */ 858 data = &data_in[flid*bs]; 859 for( kk = ii = 0; ii < bs ; ii++ ) { 860 for( jj = 0; jj < N ; jj++ ) { 861 PetscReal d = data[jj*data_stride + ii]; 862 qqc[jj*Mdata + aggID*bs + ii] = d; 863 } 864 } 865 #ifdef OUT_AGGS 866 if(llev==1) { 867 char str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^"; 868 PetscInt MM,pi,pj; 869 str[12] = col[(clid+17*mype)%6]; str[13] = sim[(clid+17*mype)%11]; 870 MM = (PetscInt)(PetscSqrtReal((PetscReal)pM)); 871 pj = gid1/MM; pi = gid1%MM; 872 fprintf(file,str,(double)pi,(double)pj); 873 /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */ 874 } 875 #endif 876 /* set fine IDs */ 877 for(kk=0;kk<bs;kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk; 878 879 aggID++; 880 } 881 882 /* pad with zeros */ 883 for( ii = asz*bs; ii < Mdata ; ii++ ) { 884 for( jj = 0; jj < N ; jj++, kk++ ) { 885 qqc[jj*Mdata + ii] = .0; 886 } 887 } 888 889 ndone += aggID; 890 /* QR */ 891 LAPACKgeqrf_( &Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO ); 892 if( INFO != 0 ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRS error"); 893 /* get R - column oriented - output B_{i+1} */ 894 { 895 PetscReal *data = &out_data[clid*nSAvec]; 896 for( jj = 0; jj < nSAvec ; jj++ ) { 897 for( ii = 0; ii < nSAvec ; ii++ ) { 898 assert(data[jj*out_data_stride + ii] == 1.e300); 899 if( ii <= jj ) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]); 900 else data[jj*out_data_stride + ii] = 0.; 901 } 902 } 903 } 904 905 /* get Q - row oriented */ 906 LAPACKungqr_( &Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO ); 907 if( INFO != 0 ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR error arg %d",-INFO); 908 909 for( ii = 0 ; ii < M ; ii++ ){ 910 for( jj = 0 ; jj < N ; jj++ ) { 911 qqr[N*ii + jj] = qqc[jj*Mdata + ii]; 912 } 913 } 914 915 /* add diagonal block of P0 */ 916 for(kk=0;kk<N;kk++) { 917 cids[kk] = N*cgid + kk; /* global col IDs in P0 */ 918 } 919 ierr = MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES); CHKERRQ(ierr); 920 921 ierr = PetscFree( qqc ); CHKERRQ(ierr); 922 ierr = PetscFree( qqr ); CHKERRQ(ierr); 923 ierr = PetscFree( TAU ); CHKERRQ(ierr); 924 ierr = PetscFree( WORK ); CHKERRQ(ierr); 925 ierr = PetscFree( fids ); CHKERRQ(ierr); 926 clid++; 927 } /* coarse agg */ 928 } /* for all fine nodes */ 929 ierr = MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 930 ierr = MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 931 932 /* ierr = MPI_Allreduce( &ndone, &ii, 1, MPIU_INT, MPIU_SUM, wcomm ); */ 933 /* MatGetSize( a_Prol, &kk, &jj ); */ 934 /* ierr = MPI_Allreduce( &minsz, &jj, 1, MPIU_INT, MPIU_MIN, wcomm ); */ 935 /* PetscPrintf(wcomm," **** [%d]%s %d total done, %d nodes (%d local done), min agg. size = %d\n",mype,__FUNCT__,ii,kk/bs,ndone,jj); */ 936 937 #ifdef OUT_AGGS 938 if(llev==1) fclose(file); 939 #endif 940 ierr = GAMGTableDestroy( &fgid_flid ); CHKERRQ(ierr); 941 942 PetscFunctionReturn(0); 943 } 944 945 /* -------------------------------------------------------------------------- */ 946 /* 947 PCGAMGgraph_AGG 948 949 Input Parameter: 950 . pc - this 951 . Amat - matrix on this fine level 952 Output Parameter: 953 . a_Gmat - 954 */ 955 #undef __FUNCT__ 956 #define __FUNCT__ "PCGAMGgraph_AGG" 957 PetscErrorCode PCGAMGgraph_AGG( PC pc, 958 const Mat Amat, 959 Mat *a_Gmat 960 ) 961 { 962 PetscErrorCode ierr; 963 PC_MG *mg = (PC_MG*)pc->data; 964 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 965 const PetscInt verbose = pc_gamg->verbose; 966 const PetscReal vfilter = pc_gamg->threshold; 967 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 968 PetscMPIInt mype,npe; 969 Mat Gmat; 970 MPI_Comm wcomm = ((PetscObject)Amat)->comm; 971 PetscBool set,flg,symm; 972 973 PetscFunctionBegin; 974 #if defined PETSC_USE_LOG 975 ierr = PetscLogEventBegin(PC_GAMGGgraph_AGG,0,0,0,0);CHKERRQ(ierr); 976 #endif 977 ierr = MPI_Comm_rank( wcomm, &mype); CHKERRQ(ierr); 978 ierr = MPI_Comm_size( wcomm, &npe); CHKERRQ(ierr); 979 980 ierr = MatIsSymmetricKnown(Amat, &set, &flg); CHKERRQ(ierr); 981 symm = (PetscBool)(pc_gamg_agg->sym_graph || !(set && flg)); 982 983 ierr = PCGAMGCreateGraph( Amat, &Gmat ); CHKERRQ( ierr ); 984 ierr = PCGAMGFilterGraph( &Gmat, vfilter, symm, verbose ); CHKERRQ( ierr ); 985 986 *a_Gmat = Gmat; 987 988 #if defined PETSC_USE_LOG 989 ierr = PetscLogEventEnd(PC_GAMGGgraph_AGG,0,0,0,0);CHKERRQ(ierr); 990 #endif 991 PetscFunctionReturn(0); 992 } 993 994 /* -------------------------------------------------------------------------- */ 995 /* 996 PCGAMGCoarsen_AGG 997 998 Input Parameter: 999 . a_pc - this 1000 Input/Output Parameter: 1001 . a_Gmat1 - graph on this fine level - coarsening can change this (squares it) 1002 Output Parameter: 1003 . agg_lists - list of aggregates 1004 */ 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "PCGAMGCoarsen_AGG" 1007 PetscErrorCode PCGAMGCoarsen_AGG( PC a_pc, 1008 Mat *a_Gmat1, 1009 PetscCoarsenData **agg_lists 1010 ) 1011 { 1012 PetscErrorCode ierr; 1013 PC_MG *mg = (PC_MG*)a_pc->data; 1014 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1015 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 1016 Mat mat,Gmat2, Gmat1 = *a_Gmat1; /* squared graph */ 1017 IS perm; 1018 PetscInt Ii,nloc,bs,n,m; 1019 PetscInt *permute; 1020 PetscBool *bIndexSet; 1021 MatCoarsen crs; 1022 MPI_Comm wcomm = ((PetscObject)Gmat1)->comm; 1023 PetscMPIInt mype,npe; 1024 1025 PetscFunctionBegin; 1026 #if defined PETSC_USE_LOG 1027 ierr = PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 1028 #endif 1029 ierr = MPI_Comm_rank( wcomm, &mype); CHKERRQ(ierr); 1030 ierr = MPI_Comm_size( wcomm, &npe); CHKERRQ(ierr); 1031 ierr = MatGetLocalSize( Gmat1, &n, &m ); CHKERRQ(ierr); 1032 ierr = MatGetBlockSize( Gmat1, &bs ); CHKERRQ(ierr); assert(bs==1); 1033 nloc = n/bs; 1034 1035 if( pc_gamg_agg->square_graph ) { 1036 ierr = MatTransposeMatMult( Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2 ); 1037 CHKERRQ(ierr); 1038 } 1039 else Gmat2 = Gmat1; 1040 1041 /* get MIS aggs */ 1042 /* randomize */ 1043 ierr = PetscMalloc( nloc*sizeof(PetscInt), &permute ); CHKERRQ(ierr); 1044 ierr = PetscMalloc( nloc*sizeof(PetscBool), &bIndexSet ); CHKERRQ(ierr); 1045 for ( Ii = 0; Ii < nloc ; Ii++ ){ 1046 bIndexSet[Ii] = PETSC_FALSE; 1047 permute[Ii] = Ii; 1048 } 1049 srand(1); /* make deterministic */ 1050 for ( Ii = 0; Ii < nloc ; Ii++ ) { 1051 PetscInt iSwapIndex = rand()%nloc; 1052 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) { 1053 PetscInt iTemp = permute[iSwapIndex]; 1054 permute[iSwapIndex] = permute[Ii]; 1055 permute[Ii] = iTemp; 1056 bIndexSet[iSwapIndex] = PETSC_TRUE; 1057 } 1058 } 1059 ierr = PetscFree( bIndexSet ); CHKERRQ(ierr); 1060 1061 ierr = ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm); 1062 CHKERRQ(ierr); 1063 #if defined PETSC_GAMG_USE_LOG 1064 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 1065 #endif 1066 ierr = MatCoarsenCreate( wcomm, &crs ); CHKERRQ(ierr); 1067 /* ierr = MatCoarsenSetType( crs, MATCOARSENMIS ); CHKERRQ(ierr); */ 1068 ierr = MatCoarsenSetFromOptions( crs ); CHKERRQ(ierr); 1069 ierr = MatCoarsenSetGreedyOrdering( crs, perm ); CHKERRQ(ierr); 1070 ierr = MatCoarsenSetAdjacency( crs, Gmat2 ); CHKERRQ(ierr); 1071 ierr = MatCoarsenSetVerbose( crs, pc_gamg->verbose ); CHKERRQ(ierr); 1072 ierr = MatCoarsenSetStrictAggs( crs, PETSC_TRUE ); CHKERRQ(ierr); 1073 ierr = MatCoarsenApply( crs ); CHKERRQ(ierr); 1074 ierr = MatCoarsenGetData( crs, agg_lists ); CHKERRQ(ierr); /* output */ 1075 ierr = MatCoarsenDestroy( &crs ); CHKERRQ(ierr); 1076 1077 ierr = ISDestroy( &perm ); CHKERRQ(ierr); 1078 ierr = PetscFree( permute ); CHKERRQ(ierr); 1079 #if defined PETSC_GAMG_USE_LOG 1080 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);CHKERRQ(ierr); 1081 #endif 1082 /* smooth aggs */ 1083 if( Gmat2 != Gmat1 ) { 1084 const PetscCoarsenData *llist = *agg_lists; 1085 ierr = smoothAggs( Gmat2, Gmat1, *agg_lists ); CHKERRQ(ierr); 1086 ierr = MatDestroy( &Gmat1 ); CHKERRQ(ierr); 1087 *a_Gmat1 = Gmat2; /* output */ 1088 ierr = PetscCDGetMat( llist, &mat ); CHKERRQ(ierr); 1089 if(mat) SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????"); 1090 } 1091 else { 1092 const PetscCoarsenData *llist = *agg_lists; 1093 /* see if we have a matrix that takes pecedence (returned from MatCoarsenAppply) */ 1094 ierr = PetscCDGetMat( llist, &mat ); CHKERRQ(ierr); 1095 if( mat ) { 1096 ierr = MatDestroy( &Gmat1 ); CHKERRQ(ierr); 1097 *a_Gmat1 = mat; /* output */ 1098 } 1099 } 1100 #if defined PETSC_USE_LOG 1101 ierr = PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);CHKERRQ(ierr); 1102 #endif 1103 PetscFunctionReturn(0); 1104 } 1105 1106 /* -------------------------------------------------------------------------- */ 1107 /* 1108 PCGAMGProlongator_AGG 1109 1110 Input Parameter: 1111 . pc - this 1112 . Amat - matrix on this fine level 1113 . Graph - used to get ghost data for nodes in 1114 . agg_lists - list of aggregates 1115 Output Parameter: 1116 . a_P_out - prolongation operator to the next level 1117 */ 1118 #undef __FUNCT__ 1119 #define __FUNCT__ "PCGAMGProlongator_AGG" 1120 PetscErrorCode PCGAMGProlongator_AGG( PC pc, 1121 const Mat Amat, 1122 const Mat Gmat, 1123 PetscCoarsenData *agg_lists, 1124 Mat *a_P_out 1125 ) 1126 { 1127 PC_MG *mg = (PC_MG*)pc->data; 1128 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1129 const PetscInt verbose = pc_gamg->verbose; 1130 const PetscInt data_cols = pc_gamg->data_cell_cols; 1131 PetscErrorCode ierr; 1132 PetscInt Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs; 1133 Mat Prol; 1134 PetscMPIInt mype, npe; 1135 MPI_Comm wcomm = ((PetscObject)Amat)->comm; 1136 const PetscInt col_bs=data_cols; 1137 PetscReal *data_w_ghost; 1138 PetscInt myCrs0, nbnodes=0, *flid_fgid; 1139 1140 PetscFunctionBegin; 1141 #if defined PETSC_USE_LOG 1142 ierr = PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1143 #endif 1144 ierr = MPI_Comm_rank( wcomm, &mype); CHKERRQ(ierr); 1145 ierr = MPI_Comm_size( wcomm, &npe); CHKERRQ(ierr); 1146 ierr = MatGetOwnershipRange( Amat, &Istart, &Iend ); CHKERRQ(ierr); 1147 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ( ierr ); 1148 nloc = (Iend-Istart)/bs; my0 = Istart/bs; assert((Iend-Istart)%bs==0); 1149 1150 /* get 'nLocalSelected' */ 1151 for( ii=0, nLocalSelected = 0 ; ii < nloc ; ii++ ){ 1152 PetscBool ise; 1153 /* filter out singletons 0 or 1? */ 1154 ierr = PetscCDEmptyAt( agg_lists, ii, &ise ); CHKERRQ(ierr); 1155 if( !ise ) nLocalSelected++; 1156 } 1157 1158 /* create prolongator, create P matrix */ 1159 ierr = MatCreateAIJ( wcomm, 1160 nloc*bs, nLocalSelected*col_bs, 1161 PETSC_DETERMINE, PETSC_DETERMINE, 1162 data_cols, PETSC_NULL, data_cols, PETSC_NULL, 1163 &Prol ); 1164 CHKERRQ(ierr); 1165 1166 /* can get all points "removed" */ 1167 ierr = MatGetSize( Prol, &kk, &ii ); CHKERRQ(ierr); 1168 if( ii==0 ) { 1169 if( verbose ) { 1170 PetscPrintf(wcomm,"[%d]%s no selected points on coarse grid\n",mype,__FUNCT__); 1171 } 1172 ierr = MatDestroy( &Prol ); CHKERRQ(ierr); 1173 *a_P_out = PETSC_NULL; /* out */ 1174 PetscFunctionReturn(0); 1175 } 1176 if( verbose ) { 1177 PetscPrintf(wcomm,"\t\t[%d]%s New grid %d nodes\n",mype,__FUNCT__,ii/col_bs); 1178 } 1179 ierr = MatGetOwnershipRangeColumn( Prol, &myCrs0, &kk ); CHKERRQ(ierr); 1180 1181 assert((kk-myCrs0)%col_bs==0); 1182 myCrs0 = myCrs0/col_bs; 1183 assert((kk/col_bs-myCrs0)==nLocalSelected); 1184 1185 /* create global vector of data in 'data_w_ghost' */ 1186 #if defined PETSC_GAMG_USE_LOG 1187 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 1188 #endif 1189 if (npe > 1) { /* */ 1190 PetscReal *tmp_gdata,*tmp_ldata,*tp2; 1191 ierr = PetscMalloc( nloc*sizeof(PetscReal), &tmp_ldata ); CHKERRQ(ierr); 1192 for( jj = 0 ; jj < data_cols ; jj++ ){ 1193 for( kk = 0 ; kk < bs ; kk++) { 1194 PetscInt ii,nnodes; 1195 const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk; 1196 for( ii = 0 ; ii < nloc ; ii++, tp += bs ){ 1197 tmp_ldata[ii] = *tp; 1198 } 1199 ierr = PCGAMGGetDataWithGhosts( Gmat, 1, tmp_ldata, &nnodes, &tmp_gdata ); 1200 CHKERRQ(ierr); 1201 if(jj==0 && kk==0) { /* now I know how many todal nodes - allocate */ 1202 ierr = PetscMalloc( nnodes*bs*data_cols*sizeof(PetscReal), &data_w_ghost ); CHKERRQ(ierr); 1203 nbnodes = bs*nnodes; 1204 } 1205 tp2 = data_w_ghost + jj*bs*nnodes + kk; 1206 for( ii = 0 ; ii < nnodes ; ii++, tp2 += bs ){ 1207 *tp2 = tmp_gdata[ii]; 1208 } 1209 ierr = PetscFree( tmp_gdata ); CHKERRQ(ierr); 1210 } 1211 } 1212 ierr = PetscFree( tmp_ldata ); CHKERRQ(ierr); 1213 } 1214 else { 1215 nbnodes = bs*nloc; 1216 data_w_ghost = (PetscReal*)pc_gamg->data; 1217 } 1218 1219 /* get P0 */ 1220 if( npe > 1 ){ 1221 PetscReal *fid_glid_loc,*fiddata; 1222 PetscInt nnodes; 1223 1224 ierr = PetscMalloc( nloc*sizeof(PetscReal), &fid_glid_loc ); CHKERRQ(ierr); 1225 for(kk=0;kk<nloc;kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk); 1226 ierr = PCGAMGGetDataWithGhosts( Gmat, 1, fid_glid_loc, &nnodes, &fiddata ); 1227 CHKERRQ(ierr); 1228 ierr = PetscMalloc( nnodes*sizeof(PetscInt), &flid_fgid ); CHKERRQ(ierr); 1229 for(kk=0;kk<nnodes;kk++) flid_fgid[kk] = (PetscInt)fiddata[kk]; 1230 ierr = PetscFree( fiddata ); CHKERRQ(ierr); 1231 assert(nnodes==nbnodes/bs); 1232 ierr = PetscFree( fid_glid_loc ); CHKERRQ(ierr); 1233 } 1234 else { 1235 ierr = PetscMalloc( nloc*sizeof(PetscInt), &flid_fgid ); CHKERRQ(ierr); 1236 for(kk=0;kk<nloc;kk++) flid_fgid[kk] = my0 + kk; 1237 } 1238 #if defined PETSC_GAMG_USE_LOG 1239 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);CHKERRQ(ierr); 1240 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1241 #endif 1242 { 1243 PetscReal *data_out = PETSC_NULL; 1244 ierr = formProl0( agg_lists, bs, data_cols, myCrs0, nbnodes, 1245 data_w_ghost, flid_fgid, &data_out, Prol ); 1246 CHKERRQ(ierr); 1247 ierr = PetscFree( pc_gamg->data ); CHKERRQ( ierr ); 1248 pc_gamg->data = data_out; 1249 pc_gamg->data_cell_rows = data_cols; 1250 pc_gamg->data_sz = data_cols*data_cols*nLocalSelected; 1251 } 1252 #if defined PETSC_GAMG_USE_LOG 1253 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);CHKERRQ(ierr); 1254 #endif 1255 if (npe > 1) ierr = PetscFree( data_w_ghost ); CHKERRQ(ierr); 1256 ierr = PetscFree( flid_fgid ); CHKERRQ(ierr); 1257 1258 /* attach block size of columns */ 1259 if( pc_gamg->col_bs_id == -1 ) { 1260 ierr = PetscObjectComposedDataRegister( &pc_gamg->col_bs_id ); assert(pc_gamg->col_bs_id != -1 ); 1261 } 1262 ierr = PetscObjectComposedDataSetInt( (PetscObject)Prol, pc_gamg->col_bs_id, data_cols ); CHKERRQ(ierr); 1263 1264 *a_P_out = Prol; /* out */ 1265 #if defined PETSC_USE_LOG 1266 ierr = PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);CHKERRQ(ierr); 1267 #endif 1268 PetscFunctionReturn(0); 1269 } 1270 1271 /* -------------------------------------------------------------------------- */ 1272 /* 1273 PCGAMGOptprol_AGG 1274 1275 Input Parameter: 1276 . pc - this 1277 . Amat - matrix on this fine level 1278 In/Output Parameter: 1279 . a_P_out - prolongation operator to the next level 1280 */ 1281 #undef __FUNCT__ 1282 #define __FUNCT__ "PCGAMGOptprol_AGG" 1283 PetscErrorCode PCGAMGOptprol_AGG( PC pc, 1284 const Mat Amat, 1285 Mat *a_P 1286 ) 1287 { 1288 PetscErrorCode ierr; 1289 PC_MG *mg = (PC_MG*)pc->data; 1290 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1291 const PetscInt verbose = pc_gamg->verbose; 1292 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; 1293 PetscInt jj; 1294 PetscMPIInt mype,npe; 1295 Mat Prol = *a_P; 1296 MPI_Comm wcomm = ((PetscObject)Amat)->comm; 1297 1298 PetscFunctionBegin; 1299 #if defined PETSC_USE_LOG 1300 ierr = PetscLogEventBegin(PC_GAMGOptprol_AGG,0,0,0,0);CHKERRQ(ierr); 1301 #endif 1302 ierr = MPI_Comm_rank( wcomm, &mype); CHKERRQ(ierr); 1303 ierr = MPI_Comm_size( wcomm, &npe); CHKERRQ(ierr); 1304 1305 /* smooth P0 */ 1306 for( jj = 0 ; jj < pc_gamg_agg->nsmooths ; jj++ ){ 1307 Mat tMat; 1308 Vec diag; 1309 PetscReal alpha, emax, emin; 1310 #if defined PETSC_GAMG_USE_LOG 1311 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1312 #endif 1313 if( jj == 0 ) { 1314 KSP eksp; 1315 Vec bb, xx; 1316 PC pc; 1317 ierr = MatGetVecs( Amat, &bb, 0 ); CHKERRQ(ierr); 1318 ierr = MatGetVecs( Amat, &xx, 0 ); CHKERRQ(ierr); 1319 { 1320 PetscRandom rctx; 1321 ierr = PetscRandomCreate(wcomm,&rctx);CHKERRQ(ierr); 1322 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 1323 ierr = VecSetRandom(bb,rctx);CHKERRQ(ierr); 1324 ierr = PetscRandomDestroy( &rctx ); CHKERRQ(ierr); 1325 } 1326 ierr = KSPCreate(wcomm,&eksp); CHKERRQ(ierr); 1327 ierr = KSPAppendOptionsPrefix( eksp, "gamg_est_"); CHKERRQ(ierr); 1328 ierr = KSPSetFromOptions( eksp ); CHKERRQ(ierr); 1329 ierr = KSPSetInitialGuessNonzero( eksp, PETSC_FALSE ); CHKERRQ(ierr); 1330 ierr = KSPSetOperators( eksp, Amat, Amat, SAME_NONZERO_PATTERN ); 1331 CHKERRQ( ierr ); 1332 ierr = KSPGetPC( eksp, &pc ); CHKERRQ( ierr ); 1333 ierr = PCSetType( pc, PCJACOBI ); CHKERRQ(ierr); /* smoother */ 1334 ierr = KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10); 1335 CHKERRQ(ierr); 1336 ierr = KSPSetNormType( eksp, KSP_NORM_NONE ); CHKERRQ(ierr); 1337 ierr = KSPSetComputeSingularValues( eksp,PETSC_TRUE ); CHKERRQ(ierr); 1338 1339 /* solve - keep stuff out of logging */ 1340 ierr = PetscLogEventDeactivate(KSP_Solve);CHKERRQ(ierr); 1341 ierr = PetscLogEventDeactivate(PC_Apply);CHKERRQ(ierr); 1342 ierr = KSPSolve( eksp, bb, xx ); CHKERRQ(ierr); 1343 ierr = PetscLogEventActivate(KSP_Solve);CHKERRQ(ierr); 1344 ierr = PetscLogEventActivate(PC_Apply);CHKERRQ(ierr); 1345 1346 ierr = KSPComputeExtremeSingularValues( eksp, &emax, &emin ); CHKERRQ(ierr); 1347 if( verbose ) { 1348 PetscPrintf(wcomm,"\t\t\t%s smooth P0: max eigen=%e min=%e PC=%s\n", 1349 __FUNCT__,emax,emin,PCJACOBI); 1350 } 1351 ierr = VecDestroy( &xx ); CHKERRQ(ierr); 1352 ierr = VecDestroy( &bb ); CHKERRQ(ierr); 1353 ierr = KSPDestroy( &eksp ); CHKERRQ(ierr); 1354 1355 if( pc_gamg->emax_id == -1 ) { 1356 ierr = PetscObjectComposedDataRegister( &pc_gamg->emax_id ); 1357 assert(pc_gamg->emax_id != -1 ); 1358 } 1359 ierr = PetscObjectComposedDataSetScalar( (PetscObject)Amat, pc_gamg->emax_id, emax ); CHKERRQ(ierr); 1360 } 1361 1362 /* smooth P1 := (I - omega/lam D^{-1}A)P0 */ 1363 ierr = MatMatMult( Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat ); CHKERRQ(ierr); 1364 ierr = MatGetVecs( Amat, &diag, 0 ); CHKERRQ(ierr); 1365 ierr = MatGetDiagonal( Amat, diag ); CHKERRQ(ierr); /* effectively PCJACOBI */ 1366 ierr = VecReciprocal( diag ); CHKERRQ(ierr); 1367 ierr = MatDiagonalScale( tMat, diag, 0 ); CHKERRQ(ierr); 1368 ierr = VecDestroy( &diag ); CHKERRQ(ierr); 1369 alpha = -1.5/emax; 1370 ierr = MatAYPX( tMat, alpha, Prol, SUBSET_NONZERO_PATTERN ); CHKERRQ(ierr); 1371 ierr = MatDestroy( &Prol ); CHKERRQ(ierr); 1372 Prol = tMat; 1373 #if defined PETSC_GAMG_USE_LOG 1374 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);CHKERRQ(ierr); 1375 #endif 1376 } 1377 #if defined PETSC_USE_LOG 1378 ierr = PetscLogEventEnd(PC_GAMGOptprol_AGG,0,0,0,0);CHKERRQ(ierr); 1379 #endif 1380 *a_P = Prol; 1381 1382 PetscFunctionReturn(0); 1383 } 1384 1385 /* -------------------------------------------------------------------------- */ 1386 /* 1387 PCCreateGAMG_AGG 1388 1389 Input Parameter: 1390 . pc - 1391 */ 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "PCCreateGAMG_AGG" 1394 PetscErrorCode PCCreateGAMG_AGG( PC pc ) 1395 { 1396 PetscErrorCode ierr; 1397 PC_MG *mg = (PC_MG*)pc->data; 1398 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1399 PC_GAMG_AGG *pc_gamg_agg; 1400 1401 PetscFunctionBegin; 1402 /* create sub context for SA */ 1403 ierr = PetscNewLog( pc, PC_GAMG_AGG, &pc_gamg_agg ); CHKERRQ(ierr); 1404 assert(!pc_gamg->subctx); 1405 pc_gamg->subctx = pc_gamg_agg; 1406 1407 pc->ops->setfromoptions = PCSetFromOptions_GAMG_AGG; 1408 pc->ops->destroy = PCDestroy_AGG; 1409 /* reset does not do anything; setup not virtual */ 1410 1411 /* set internal function pointers */ 1412 pc_gamg->graph = PCGAMGgraph_AGG; 1413 pc_gamg->coarsen = PCGAMGCoarsen_AGG; 1414 pc_gamg->prolongator = PCGAMGProlongator_AGG; 1415 pc_gamg->optprol = PCGAMGOptprol_AGG; 1416 1417 pc_gamg->createdefaultdata = PCSetData_AGG; 1418 1419 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 1420 "PCSetCoordinates_C", 1421 "PCSetCoordinates_AGG", 1422 PCSetCoordinates_AGG); 1423 PetscFunctionReturn(0); 1424 } 1425