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