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