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