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