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