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