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