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