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