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