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