Actual source code: agg.c
petsc-3.4.5 2014-06-29
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: PetscBool 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: PetscErrorCode PCGAMGSetNSmooths_GAMG(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: PetscErrorCode PCGAMGSetSymGraph_GAMG(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, PetscBool n)
120: {
125: PetscTryMethod(pc,"PCGAMGSetSquareGraph_C",(PC,PetscBool),(pc,n));
126: return(0);
127: }
131: PetscErrorCode PCGAMGSetSquareGraph_GAMG(PC pc, PetscBool 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: PetscErrorCode PCSetFromOptions_GAMG_AGG(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;
157: PetscBool flag;
161: PetscOptionsHead("GAMG-AGG options");
162: {
163: /* -pc_gamg_agg_nsmooths */
164: pc_gamg_agg->nsmooths = 1;
166: PetscOptionsInt("-pc_gamg_agg_nsmooths",
167: "smoothing steps for smoothed aggregation, usually 1 (0)",
168: "PCGAMGSetNSmooths",
169: pc_gamg_agg->nsmooths,
170: &pc_gamg_agg->nsmooths,
171: &flag);
173: /* -pc_gamg_sym_graph */
174: pc_gamg_agg->sym_graph = PETSC_FALSE;
176: PetscOptionsBool("-pc_gamg_sym_graph",
177: "Set for asymmetric matrices",
178: "PCGAMGSetSymGraph",
179: pc_gamg_agg->sym_graph,
180: &pc_gamg_agg->sym_graph,
181: &flag);
183: /* -pc_gamg_square_graph */
184: pc_gamg_agg->square_graph = PETSC_TRUE;
186: PetscOptionsBool("-pc_gamg_square_graph",
187: "For faster coarsening and lower coarse grid complexity",
188: "PCGAMGSetSquareGraph",
189: pc_gamg_agg->square_graph,
190: &pc_gamg_agg->square_graph,
191: &flag);
192: }
193: PetscOptionsTail();
194: return(0);
195: }
197: /* -------------------------------------------------------------------------- */
198: /*
199: PCDestroy_AGG
201: Input Parameter:
202: . pc -
203: */
206: PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
207: {
209: PC_MG *mg = (PC_MG*)pc->data;
210: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
213: PetscFree(pc_gamg->subctx);
214: return(0);
215: }
217: /* -------------------------------------------------------------------------- */
218: /*
219: PCSetCoordinates_AGG
220: - collective
222: Input Parameter:
223: . pc - the preconditioner context
224: . ndm - dimesion of data (used for dof/vertex for Stokes)
225: . a_nloc - number of vertices local
226: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
227: */
231: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
232: {
233: PC_MG *mg = (PC_MG*)pc->data;
234: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
236: PetscInt arrsz,kk,ii,jj,nloc,ndatarows,ndf;
237: Mat mat = pc->pmat;
242: nloc = a_nloc;
244: /* SA: null space vectors */
245: MatGetBlockSize(mat, &ndf); /* this does not work for Stokes */
246: if (coords && ndf==1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
247: else if (coords) {
248: if (ndm > ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"degrees of motion %d > block size %d",ndm,ndf);
249: pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */
250: if (ndm != ndf) {
251: 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);
252: }
253: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
254: pc_gamg->data_cell_rows = ndatarows = ndf;
255: 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);
256: arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;
258: /* create data - syntactic sugar that should be refactored at some point */
259: if (pc_gamg->data==0 || (pc_gamg->data_sz != arrsz)) {
260: PetscFree(pc_gamg->data);
261: PetscMalloc((arrsz+1)*sizeof(PetscReal), &pc_gamg->data);
262: }
263: /* copy data in - column oriented */
264: for (kk=0; kk<nloc; kk++) {
265: const PetscInt M = nloc*pc_gamg->data_cell_rows; /* stride into data */
266: PetscReal *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */
267: if (pc_gamg->data_cell_cols==1) *data = 1.0;
268: else {
269: /* translational modes */
270: for (ii=0;ii<ndatarows;ii++) {
271: for (jj=0;jj<ndatarows;jj++) {
272: if (ii==jj)data[ii*M + jj] = 1.0;
273: else data[ii*M + jj] = 0.0;
274: }
275: }
277: /* rotational modes */
278: if (coords) {
279: if (ndm == 2) {
280: data += 2*M;
281: data[0] = -coords[2*kk+1];
282: data[1] = coords[2*kk];
283: } else {
284: data += 3*M;
285: data[0] = 0.0; data[M+0] = coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1];
286: data[1] = -coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = coords[3*kk];
287: data[2] = coords[3*kk+1]; data[M+2] = -coords[3*kk]; data[2*M+2] = 0.0;
288: }
289: }
290: }
291: }
293: pc_gamg->data_sz = arrsz;
294: return(0);
295: }
297: typedef PetscInt NState;
298: static const NState NOT_DONE=-2;
299: static const NState DELETED =-1;
300: static const NState REMOVED =-3;
301: #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED)
303: /* -------------------------------------------------------------------------- */
304: /*
305: smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific
306: - AGG-MG specific: clears singletons out of 'selected_2'
308: Input Parameter:
309: . Gmat_2 - glabal matrix of graph (data not defined)
310: . Gmat_1 - base graph to grab with
311: Input/Output Parameter:
312: . aggs_2 - linked list of aggs with gids)
313: */
316: static PetscErrorCode smoothAggs(const Mat Gmat_2, /* base (squared) graph */
317: const Mat Gmat_1, /* base graph */
318: /* const IS selected_2, [nselected local] selected vertices */
319: PetscCoarsenData *aggs_2) /* [nselected local] global ID of aggregate */
320: {
322: PetscBool isMPI;
323: Mat_SeqAIJ *matA_1, *matB_1=0, *matA_2, *matB_2=0;
324: MPI_Comm comm;
325: PetscMPIInt rank,size;
326: PetscInt lid,*ii,*idx,ix,Iend,my0,kk,n,j;
327: Mat_MPIAIJ *mpimat_2 = 0, *mpimat_1=0;
328: const PetscInt nloc = Gmat_2->rmap->n;
329: PetscScalar *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid;
330: PetscInt *lid_cprowID_1;
331: NState *lid_state;
332: Vec ghost_par_orig2;
335: PetscObjectGetComm((PetscObject)Gmat_2,&comm);
336: MPI_Comm_rank(comm, &rank);
337: MPI_Comm_size(comm, &size);
338: MatGetOwnershipRange(Gmat_1,&my0,&Iend);
340: if (PETSC_FALSE) {
341: PetscViewer viewer; char fname[32]; static int llev=0;
342: sprintf(fname,"Gmat2_%d.m",llev++);
343: PetscViewerASCIIOpen(comm,fname,&viewer);
344: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
345: MatView(Gmat_2, viewer);
346: PetscViewerDestroy(&viewer);
347: }
349: /* get submatrices */
350: PetscObjectTypeCompare((PetscObject)Gmat_1, MATMPIAIJ, &isMPI);
351: if (isMPI) {
352: /* grab matrix objects */
353: mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data;
354: mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data;
355: matA_1 = (Mat_SeqAIJ*)mpimat_1->A->data;
356: matB_1 = (Mat_SeqAIJ*)mpimat_1->B->data;
357: matA_2 = (Mat_SeqAIJ*)mpimat_2->A->data;
358: matB_2 = (Mat_SeqAIJ*)mpimat_2->B->data;
360: /* force compressed row storage for B matrix in AuxMat */
361: matB_1->compressedrow.check = PETSC_TRUE;
363: MatCheckCompressedRow(mpimat_1->B,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0);
365: PetscMalloc(nloc*sizeof(PetscInt), &lid_cprowID_1);
366: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
367: for (ix=0; ix<matB_1->compressedrow.nrows; ix++) {
368: PetscInt lid = matB_1->compressedrow.rindex[ix];
369: lid_cprowID_1[lid] = ix;
370: }
371: } else {
372: matA_1 = (Mat_SeqAIJ*)Gmat_1->data;
373: matA_2 = (Mat_SeqAIJ*)Gmat_2->data;
374: lid_cprowID_1 = NULL;
375: }
376: if (!(matA_1 && !matA_1->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matA_1 && !matA_1->compressedrow.use)");
377: if (!(matB_1==0 || matB_1->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matB_1==0 || matB_1->compressedrow.use)");
378: if (!(matA_2 && !matA_2->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matA_2 && !matA_2->compressedrow.use)");
379: if (!(matB_2==0 || matB_2->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matB_2==0 || matB_2->compressedrow.use)");
381: /* get state of locals and selected gid for deleted */
382: PetscMalloc(nloc*sizeof(NState), &lid_state);
383: PetscMalloc(nloc*sizeof(PetscScalar), &lid_parent_gid);
384: for (lid = 0; lid < nloc; lid++) {
385: lid_parent_gid[lid] = -1.0;
386: lid_state[lid] = DELETED;
387: }
389: /* set lid_state */
390: for (lid = 0; lid < nloc; lid++) {
391: PetscCDPos pos;
392: PetscCDGetHeadPos(aggs_2,lid,&pos);
393: if (pos) {
394: PetscInt gid1;
396: PetscLLNGetID(pos, &gid1);
397: if (gid1 != lid+my0) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"gid1 %D != lid %D + my0 %D",gid1,lid,my0);
398: lid_state[lid] = gid1;
399: }
400: }
402: /* map local to selected local, DELETED means a ghost owns it */
403: for (lid=kk=0; lid<nloc; lid++) {
404: NState state = lid_state[lid];
405: if (IS_SELECTED(state)) {
406: PetscCDPos pos;
407: PetscCDGetHeadPos(aggs_2,lid,&pos);
408: while (pos) {
409: PetscInt gid1;
410: PetscLLNGetID(pos, &gid1);
411: PetscCDGetNextPos(aggs_2,lid,&pos);
413: if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0);
414: }
415: }
416: }
417: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
418: if (isMPI) {
419: Vec tempVec;
420: /* get 'cpcol_1_state' */
421: MatGetVecs(Gmat_1, &tempVec, 0);
422: for (kk=0,j=my0; kk<nloc; kk++,j++) {
423: PetscScalar v = (PetscScalar)lid_state[kk];
424: VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
425: }
426: VecAssemblyBegin(tempVec);
427: VecAssemblyEnd(tempVec);
428: VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
429: VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
430: VecGetArray(mpimat_1->lvec, &cpcol_1_state);
431: /* get 'cpcol_2_state' */
432: VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
433: VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
434: VecGetArray(mpimat_2->lvec, &cpcol_2_state);
435: /* get 'cpcol_2_par_orig' */
436: for (kk=0,j=my0; kk<nloc; kk++,j++) {
437: PetscScalar v = (PetscScalar)lid_parent_gid[kk];
438: VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
439: }
440: VecAssemblyBegin(tempVec);
441: VecAssemblyEnd(tempVec);
442: VecDuplicate(mpimat_2->lvec, &ghost_par_orig2);
443: VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
444: VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
445: VecGetArray(ghost_par_orig2, &cpcol_2_par_orig);
447: VecDestroy(&tempVec);
448: } /* ismpi */
450: /* doit */
451: for (lid=0; lid<nloc; lid++) {
452: NState state = lid_state[lid];
453: if (IS_SELECTED(state)) {
454: /* steal locals */
455: ii = matA_1->i; n = ii[lid+1] - ii[lid];
456: idx = matA_1->j + ii[lid];
457: for (j=0; j<n; j++) {
458: PetscInt lidj = idx[j], sgid;
459: NState statej = lid_state[lidj];
460: if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */
461: lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */
462: if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */
463: PetscInt hav=0,slid=sgid-my0,gidj=lidj+my0;
464: PetscCDPos pos,last=NULL;
465: /* looking for local from local so id_llist_2 works */
466: PetscCDGetHeadPos(aggs_2,slid,&pos);
467: while (pos) {
468: PetscInt gid;
469: PetscLLNGetID(pos, &gid);
470: if (gid == gidj) {
471: if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
472: PetscCDRemoveNextNode(aggs_2, slid, last);
473: PetscCDAppendNode(aggs_2, lid, pos);
474: hav = 1;
475: break;
476: } else last = pos;
478: PetscCDGetNextPos(aggs_2,slid,&pos);
479: }
480: if (hav!=1) {
481: if (!hav) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
482: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
483: }
484: } else { /* I'm stealing this local, owned by a ghost */
485: 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 0.0' if the matrix is structurally symmetric.");
486: PetscCDAppendID(aggs_2, lid, lidj+my0);
487: }
488: }
489: } /* local neighbors */
490: } else if (state == DELETED && lid_cprowID_1) {
491: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
492: /* see if I have a selected ghost neighbor that will steal me */
493: if ((ix=lid_cprowID_1[lid]) != -1) {
494: ii = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix];
495: idx = matB_1->j + ii[ix];
496: for (j=0; j<n; j++) {
497: PetscInt cpid = idx[j];
498: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
499: if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */
500: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
501: if (sgidold>=my0 && sgidold<Iend) { /* this was mine */
502: PetscInt hav=0,oldslidj=sgidold-my0;
503: PetscCDPos pos,last=NULL;
504: /* remove from 'oldslidj' list */
505: PetscCDGetHeadPos(aggs_2,oldslidj,&pos);
506: while (pos) {
507: PetscInt gid;
508: PetscLLNGetID(pos, &gid);
509: if (lid+my0 == gid) {
510: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
511: if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
512: PetscCDRemoveNextNode(aggs_2, oldslidj, last);
513: /* ghost (PetscScalar)statej will add this later */
514: hav = 1;
515: break;
516: } else last = pos;
518: PetscCDGetNextPos(aggs_2,oldslidj,&pos);
519: }
520: if (hav!=1) {
521: if (hav==0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
522: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
523: }
524: } else {
525: /* ghosts remove this later */
526: }
527: }
528: }
529: }
530: } /* selected/deleted */
531: } /* node loop */
533: if (isMPI) {
534: PetscScalar *cpcol_2_parent,*cpcol_2_gid;
535: Vec tempVec,ghostgids2,ghostparents2;
536: PetscInt cpid,nghost_2;
537: GAMGHashTable gid_cpid;
539: VecGetSize(mpimat_2->lvec, &nghost_2);
540: MatGetVecs(Gmat_2, &tempVec, 0);
542: /* get 'cpcol_2_parent' */
543: for (kk=0,j=my0; kk<nloc; kk++,j++) {
544: VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES);
545: }
546: VecAssemblyBegin(tempVec);
547: VecAssemblyEnd(tempVec);
548: VecDuplicate(mpimat_2->lvec, &ghostparents2);
549: VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
550: VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
551: VecGetArray(ghostparents2, &cpcol_2_parent);
553: /* get 'cpcol_2_gid' */
554: for (kk=0,j=my0; kk<nloc; kk++,j++) {
555: PetscScalar v = (PetscScalar)j;
556: VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
557: }
558: VecAssemblyBegin(tempVec);
559: VecAssemblyEnd(tempVec);
560: VecDuplicate(mpimat_2->lvec, &ghostgids2);
561: VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
562: VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
563: VecGetArray(ghostgids2, &cpcol_2_gid);
565: VecDestroy(&tempVec);
567: /* look for deleted ghosts and add to table */
568: GAMGTableCreate(2*nghost_2+1, &gid_cpid);
569: for (cpid = 0; cpid < nghost_2; cpid++) {
570: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
571: if (state==DELETED) {
572: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
573: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
574: if (sgid_old == -1 && sgid_new != -1) {
575: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
576: GAMGTableAdd(&gid_cpid, gid, cpid);
577: }
578: }
579: }
581: /* look for deleted ghosts and see if they moved - remove it */
582: for (lid=0; lid<nloc; lid++) {
583: NState state = lid_state[lid];
584: if (IS_SELECTED(state)) {
585: PetscCDPos pos,last=NULL;
586: /* look for deleted ghosts and see if they moved */
587: PetscCDGetHeadPos(aggs_2,lid,&pos);
588: while (pos) {
589: PetscInt gid;
590: PetscLLNGetID(pos, &gid);
592: if (gid < my0 || gid >= Iend) {
593: GAMGTableFind(&gid_cpid, gid, &cpid);
594: if (cpid != -1) {
595: /* a moved ghost - */
596: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
597: PetscCDRemoveNextNode(aggs_2, lid, last);
598: } else last = pos;
599: } else last = pos;
601: PetscCDGetNextPos(aggs_2,lid,&pos);
602: } /* loop over list of deleted */
603: } /* selected */
604: }
605: GAMGTableDestroy(&gid_cpid);
607: /* look at ghosts, see if they changed - and it */
608: for (cpid = 0; cpid < nghost_2; cpid++) {
609: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
610: if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
611: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
612: PetscInt slid_new=sgid_new-my0,hav=0;
613: PetscCDPos pos;
614: /* search for this gid to see if I have it */
615: PetscCDGetHeadPos(aggs_2,slid_new,&pos);
616: while (pos) {
617: PetscInt gidj;
618: PetscLLNGetID(pos, &gidj);
619: PetscCDGetNextPos(aggs_2,slid_new,&pos);
621: if (gidj == gid) { hav = 1; break; }
622: }
623: if (hav != 1) {
624: /* insert 'flidj' into head of llist */
625: PetscCDAppendID(aggs_2, slid_new, gid);
626: }
627: }
628: }
630: VecRestoreArray(mpimat_1->lvec, &cpcol_1_state);
631: VecRestoreArray(mpimat_2->lvec, &cpcol_2_state);
632: VecRestoreArray(ghostparents2, &cpcol_2_parent);
633: VecRestoreArray(ghostgids2, &cpcol_2_gid);
634: PetscFree(lid_cprowID_1);
635: VecDestroy(&ghostgids2);
636: VecDestroy(&ghostparents2);
637: VecDestroy(&ghost_par_orig2);
638: }
640: PetscFree(lid_parent_gid);
641: PetscFree(lid_state);
642: return(0);
643: }
645: /* -------------------------------------------------------------------------- */
646: /*
647: PCSetData_AGG - called if data is not set with PCSetCoordinates.
648: Looks in Mat for near null space.
649: Does not work for Stokes
651: Input Parameter:
652: . pc -
653: . a_A - matrix to get (near) null space out of.
654: */
657: PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
658: {
660: PC_MG *mg = (PC_MG*)pc->data;
661: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
662: MatNullSpace mnull;
665: MatGetNearNullSpace(a_A, &mnull);
666: if (!mnull) {
667: PetscInt bs,NN,MM;
668: MatGetBlockSize(a_A, &bs);
669: MatGetLocalSize(a_A, &MM, &NN);
670: if (MM % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MM %D must be divisible by bs %D",MM,bs);
671: PCSetCoordinates_AGG(pc, bs, MM/bs, NULL);
672: } else {
673: PetscReal *nullvec;
674: PetscBool has_const;
675: PetscInt i,j,mlocal,nvec,bs;
676: const Vec *vecs; const PetscScalar *v;
677: MatGetLocalSize(a_A,&mlocal,NULL);
678: MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
679: pc_gamg->data_sz = (nvec+!!has_const)*mlocal;
680: PetscMalloc((nvec+!!has_const)*mlocal*sizeof(*nullvec),&nullvec);
681: if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
682: for (i=0; i<nvec; i++) {
683: VecGetArrayRead(vecs[i],&v);
684: for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]);
685: VecRestoreArrayRead(vecs[i],&v);
686: }
687: pc_gamg->data = nullvec;
688: pc_gamg->data_cell_cols = (nvec+!!has_const);
690: MatGetBlockSize(a_A, &bs);
692: pc_gamg->data_cell_rows = bs;
693: }
694: return(0);
695: }
697: /* -------------------------------------------------------------------------- */
698: /*
699: formProl0
701: Input Parameter:
702: . agg_llists - list of arrays with aggregates
703: . bs - block size
704: . nSAvec - column bs of new P
705: . my0crs - global index of start of locals
706: . data_stride - bs*(nloc nodes + ghost nodes)
707: . data_in[data_stride*nSAvec] - local data on fine grid
708: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
709: Output Parameter:
710: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
711: . a_Prol - prolongation operator
712: */
715: static PetscErrorCode formProl0(const PetscCoarsenData *agg_llists, /* list from selected vertices of aggregate unselected vertices */
716: const PetscInt bs, /* (row) block size */
717: const PetscInt nSAvec, /* column bs */
718: const PetscInt my0crs, /* global index of start of locals */
719: const PetscInt data_stride, /* (nloc+nghost)*bs */
720: PetscReal data_in[], /* [data_stride][nSAvec] */
721: const PetscInt flid_fgid[], /* [data_stride/bs] */
722: PetscReal **a_data_out,
723: Mat a_Prol) /* prolongation operator (output)*/
724: {
726: PetscInt Istart,my0,Iend,nloc,clid,flid,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride;
727: MPI_Comm comm;
728: PetscMPIInt rank, size;
729: PetscReal *out_data;
730: PetscCDPos pos;
731: GAMGHashTable fgid_flid;
733: /* #define OUT_AGGS */
734: #if defined(OUT_AGGS)
735: static PetscInt llev = 0; char fname[32]; FILE *file = NULL; PetscInt pM;
736: #endif
739: PetscObjectGetComm((PetscObject)a_Prol,&comm);
740: MPI_Comm_rank(comm,&rank);
741: MPI_Comm_size(comm,&size);
742: MatGetOwnershipRange(a_Prol, &Istart, &Iend);
743: nloc = (Iend-Istart)/bs; my0 = Istart/bs;
744: if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Iend %D - Istart %d must be divisible by bs %D",Iend,Istart,bs);
745: Iend /= bs;
746: nghosts = data_stride/bs - nloc;
748: GAMGTableCreate(2*nghosts+1, &fgid_flid);
749: for (kk=0; kk<nghosts; kk++) {
750: GAMGTableAdd(&fgid_flid, flid_fgid[nloc+kk], nloc+kk);
751: }
753: #if defined(OUT_AGGS)
754: sprintf(fname,"aggs_%d_%d.m",llev++,rank);
755: if (llev==1) file = fopen(fname,"w");
756: MatGetSize(a_Prol, &pM, &jj);
757: #endif
759: /* count selected -- same as number of cols of P */
760: for (nSelected=mm=0; mm<nloc; mm++) {
761: PetscBool ise;
762: PetscCDEmptyAt(agg_llists, mm, &ise);
763: if (!ise) nSelected++;
764: }
765: MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);
766: if ((ii/nSAvec) != my0crs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D /nSAvec %D != my0crs %D",ii,nSAvec,my0crs);
767: if (nSelected != (jj-ii)/nSAvec) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nSelected %D != (jj %D - ii %D)/nSAvec %D",nSelected,jj,ii,nSAvec);
769: /* aloc space for coarse point data (output) */
770: out_data_stride = nSelected*nSAvec;
772: PetscMalloc(out_data_stride*nSAvec*sizeof(PetscReal), &out_data);
773: for (ii=0;ii<out_data_stride*nSAvec;ii++) out_data[ii]=PETSC_MAX_REAL;
774: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
776: /* find points and set prolongation */
777: minsz = 100;
778: ndone = 0;
779: for (mm = clid = 0; mm < nloc; mm++) {
780: PetscCDSizeAt(agg_llists, mm, &jj);
781: if (jj > 0) {
782: const PetscInt lid = mm, cgid = my0crs + clid;
783: PetscInt cids[100]; /* max bs */
784: PetscBLASInt asz =jj,M=asz*bs,N=nSAvec,INFO;
785: PetscBLASInt Mdata=M+((N-M>0) ? N-M : 0),LDA=Mdata,LWORK=N*bs;
786: PetscScalar *qqc,*qqr,*TAU,*WORK;
787: PetscInt *fids;
788: PetscReal *data;
789: /* count agg */
790: if (asz<minsz) minsz = asz;
792: /* get block */
793: PetscMalloc((Mdata*N)*sizeof(PetscScalar), &qqc);
794: PetscMalloc((M*N)*sizeof(PetscScalar), &qqr);
795: PetscMalloc(N*sizeof(PetscScalar), &TAU);
796: PetscMalloc(LWORK*sizeof(PetscScalar), &WORK);
797: PetscMalloc(M*sizeof(PetscInt), &fids);
799: aggID = 0;
800: PetscCDGetHeadPos(agg_llists,lid,&pos);
801: while (pos) {
802: PetscInt gid1;
803: PetscLLNGetID(pos, &gid1);
804: PetscCDGetNextPos(agg_llists,lid,&pos);
806: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
807: else {
808: GAMGTableFind(&fgid_flid, gid1, &flid);
809: if (flid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot find gid1 in table");
810: }
811: /* copy in B_i matrix - column oriented */
812: data = &data_in[flid*bs];
813: for (kk = ii = 0; ii < bs; ii++) {
814: for (jj = 0; jj < N; jj++) {
815: PetscReal d = data[jj*data_stride + ii];
816: qqc[jj*Mdata + aggID*bs + ii] = d;
817: }
818: }
819: #if defined(OUT_AGGS)
820: if (llev==1) {
821: char str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^";
822: PetscInt MM,pi,pj;
823: str[12] = col[(clid+17*rank)%6]; str[13] = sim[(clid+17*rank)%11];
824: MM = (PetscInt)(PetscSqrtReal((PetscReal)pM));
825: pj = gid1/MM; pi = gid1%MM;
826: fprintf(file,str,(double)pi,(double)pj);
827: /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */
828: }
829: #endif
830: /* set fine IDs */
831: for (kk=0; kk<bs; kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk;
833: aggID++;
834: }
836: /* pad with zeros */
837: for (ii = asz*bs; ii < Mdata; ii++) {
838: for (jj = 0; jj < N; jj++, kk++) {
839: qqc[jj*Mdata + ii] = .0;
840: }
841: }
843: ndone += aggID;
844: /* QR */
845: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
846: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
847: PetscFPTrapPop();
848: if (INFO != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xGEQRF error");
849: /* get R - column oriented - output B_{i+1} */
850: {
851: PetscReal *data = &out_data[clid*nSAvec];
852: for (jj = 0; jj < nSAvec; jj++) {
853: for (ii = 0; ii < nSAvec; ii++) {
854: 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);
855: if (ii <= jj) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]);
856: else data[jj*out_data_stride + ii] = 0.;
857: }
858: }
859: }
861: /* get Q - row oriented */
862: PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
863: if (INFO != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xORGQR error arg %d",-INFO);
865: for (ii = 0; ii < M; ii++) {
866: for (jj = 0; jj < N; jj++) {
867: qqr[N*ii + jj] = qqc[jj*Mdata + ii];
868: }
869: }
871: /* add diagonal block of P0 */
872: for (kk=0; kk<N; kk++) {
873: cids[kk] = N*cgid + kk; /* global col IDs in P0 */
874: }
875: MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);
877: PetscFree(qqc);
878: PetscFree(qqr);
879: PetscFree(TAU);
880: PetscFree(WORK);
881: PetscFree(fids);
882: clid++;
883: } /* coarse agg */
884: } /* for all fine nodes */
885: MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
886: MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);
888: /* MPI_Allreduce(&ndone, &ii, 1, MPIU_INT, MPIU_SUM, comm); */
889: /* MatGetSize(a_Prol, &kk, &jj); */
890: /* MPI_Allreduce(&minsz, &jj, 1, MPIU_INT, MPIU_MIN, comm); */
891: /* PetscPrintf(comm," **** [%d]%s %d total done, %d nodes (%d local done), min agg. size = %d\n",rank,__FUNCT__,ii,kk/bs,ndone,jj); */
893: #if defined(OUT_AGGS)
894: if (llev==1) fclose(file);
895: #endif
896: GAMGTableDestroy(&fgid_flid);
897: return(0);
898: }
900: /* -------------------------------------------------------------------------- */
901: /*
902: PCGAMGgraph_AGG
904: Input Parameter:
905: . pc - this
906: . Amat - matrix on this fine level
907: Output Parameter:
908: . a_Gmat -
909: */
912: PetscErrorCode PCGAMGgraph_AGG(PC pc,const Mat Amat,Mat *a_Gmat)
913: {
914: PetscErrorCode ierr;
915: PC_MG *mg = (PC_MG*)pc->data;
916: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
917: const PetscInt verbose = pc_gamg->verbose;
918: const PetscReal vfilter = pc_gamg->threshold;
919: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
920: PetscMPIInt rank,size;
921: Mat Gmat;
922: MPI_Comm comm;
923: PetscBool /* set,flg , */ symm;
926: PetscObjectGetComm((PetscObject)Amat,&comm);
927: #if defined PETSC_USE_LOG
928: PetscLogEventBegin(PC_GAMGGgraph_AGG,0,0,0,0);
929: #endif
930: MPI_Comm_rank(comm, &rank);
931: MPI_Comm_size(comm, &size);
933: /* MatIsSymmetricKnown(Amat, &set, &flg); || !(set && flg) -- this causes lot of symm calls */
934: symm = (PetscBool)(pc_gamg_agg->sym_graph); /* && !pc_gamg_agg->square_graph; */
936: PCGAMGCreateGraph(Amat, &Gmat);
937: PCGAMGFilterGraph(&Gmat, vfilter, symm, verbose);
939: *a_Gmat = Gmat;
941: #if defined PETSC_USE_LOG
942: PetscLogEventEnd(PC_GAMGGgraph_AGG,0,0,0,0);
943: #endif
944: return(0);
945: }
947: /* -------------------------------------------------------------------------- */
948: /*
949: PCGAMGCoarsen_AGG
951: Input Parameter:
952: . a_pc - this
953: Input/Output Parameter:
954: . a_Gmat1 - graph on this fine level - coarsening can change this (squares it)
955: Output Parameter:
956: . agg_lists - list of aggregates
957: */
960: PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc,Mat *a_Gmat1,PetscCoarsenData **agg_lists)
961: {
963: PC_MG *mg = (PC_MG*)a_pc->data;
964: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
965: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
966: Mat mat,Gmat2, Gmat1 = *a_Gmat1; /* squared graph */
967: const PetscInt verbose = pc_gamg->verbose;
968: IS perm;
969: PetscInt Ii,nloc,bs,n,m;
970: PetscInt *permute;
971: PetscBool *bIndexSet;
972: MatCoarsen crs;
973: MPI_Comm comm;
974: PetscMPIInt rank,size;
977: #if defined PETSC_USE_LOG
978: PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);
979: #endif
980: PetscObjectGetComm((PetscObject)Gmat1,&comm);
981: MPI_Comm_rank(comm, &rank);
982: MPI_Comm_size(comm, &size);
983: MatGetLocalSize(Gmat1, &n, &m);
984: MatGetBlockSize(Gmat1, &bs);
985: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bs %D must be 1",bs);
986: nloc = n/bs;
988: if (pc_gamg_agg->square_graph) {
989: if (verbose > 1) PetscPrintf(comm,"[%d]%s square graph\n",rank,__FUNCT__);
990: /* MatMatTransposeMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2); */
991: MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);
992: if (verbose > 2) {
993: PetscPrintf(comm,"[%d]%s square graph done\n",rank,__FUNCT__);
994: /* check for symetry */
995: if (verbose > 4) {
997: }
998: }
999: } else Gmat2 = Gmat1;
1001: /* get MIS aggs */
1002: /* randomize */
1003: PetscMalloc(nloc*sizeof(PetscInt), &permute);
1004: PetscMalloc(nloc*sizeof(PetscBool), &bIndexSet);
1005: for (Ii = 0; Ii < nloc; Ii++) {
1006: bIndexSet[Ii] = PETSC_FALSE;
1007: permute[Ii] = Ii;
1008: }
1009: srand(1); /* make deterministic */
1010: for (Ii = 0; Ii < nloc; Ii++) {
1011: PetscInt iSwapIndex = rand()%nloc;
1012: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1013: PetscInt iTemp = permute[iSwapIndex];
1014: permute[iSwapIndex] = permute[Ii];
1015: permute[Ii] = iTemp;
1016: bIndexSet[iSwapIndex] = PETSC_TRUE;
1017: }
1018: }
1019: PetscFree(bIndexSet);
1021: if (verbose > 1) PetscPrintf(comm,"[%d]%s coarsen graph\n",rank,__FUNCT__);
1023: ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);
1024: #if defined PETSC_GAMG_USE_LOG
1025: PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);
1026: #endif
1027: MatCoarsenCreate(comm, &crs);
1028: /* MatCoarsenSetType(crs, MATCOARSENMIS); */
1029: MatCoarsenSetFromOptions(crs);
1030: MatCoarsenSetGreedyOrdering(crs, perm);
1031: MatCoarsenSetAdjacency(crs, Gmat2);
1032: MatCoarsenSetVerbose(crs, pc_gamg->verbose);
1033: MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
1034: MatCoarsenApply(crs);
1035: MatCoarsenGetData(crs, agg_lists); /* output */
1036: MatCoarsenDestroy(&crs);
1038: ISDestroy(&perm);
1039: PetscFree(permute);
1040: #if defined PETSC_GAMG_USE_LOG
1041: PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);
1042: #endif
1043: if (verbose > 2) PetscPrintf(comm,"[%d]%s coarsen graph done\n",rank,__FUNCT__);
1045: /* smooth aggs */
1046: if (Gmat2 != Gmat1) {
1047: const PetscCoarsenData *llist = *agg_lists;
1048: smoothAggs(Gmat2, Gmat1, *agg_lists);
1049: MatDestroy(&Gmat1);
1050: *a_Gmat1 = Gmat2; /* output */
1051: PetscCDGetMat(llist, &mat);
1052: if (mat) SETERRQ(comm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????");
1053: } else {
1054: const PetscCoarsenData *llist = *agg_lists;
1055: /* see if we have a matrix that takes precedence (returned from MatCoarsenAppply) */
1056: PetscCDGetMat(llist, &mat);
1057: if (mat) {
1058: MatDestroy(&Gmat1);
1059: *a_Gmat1 = mat; /* output */
1060: }
1061: }
1062: #if defined PETSC_USE_LOG
1063: PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);
1064: #endif
1065: if (verbose > 2) PetscPrintf(comm,"[%d]%s PCGAMGCoarsen_AGG done\n",rank,__FUNCT__);
1066: return(0);
1067: }
1069: /* -------------------------------------------------------------------------- */
1070: /*
1071: PCGAMGProlongator_AGG
1073: Input Parameter:
1074: . pc - this
1075: . Amat - matrix on this fine level
1076: . Graph - used to get ghost data for nodes in
1077: . agg_lists - list of aggregates
1078: Output Parameter:
1079: . a_P_out - prolongation operator to the next level
1080: */
1083: PetscErrorCode PCGAMGProlongator_AGG(PC pc,const Mat Amat,const Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
1084: {
1085: PC_MG *mg = (PC_MG*)pc->data;
1086: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1087: const PetscInt verbose = pc_gamg->verbose;
1088: const PetscInt data_cols = pc_gamg->data_cell_cols;
1090: PetscInt Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs;
1091: Mat Prol;
1092: PetscMPIInt rank, size;
1093: MPI_Comm comm;
1094: const PetscInt col_bs = data_cols;
1095: PetscReal *data_w_ghost;
1096: PetscInt myCrs0, nbnodes=0, *flid_fgid;
1099: PetscObjectGetComm((PetscObject)Amat,&comm);
1100: #if defined PETSC_USE_LOG
1101: PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);
1102: #endif
1103: MPI_Comm_rank(comm, &rank);
1104: MPI_Comm_size(comm, &size);
1105: MatGetOwnershipRange(Amat, &Istart, &Iend);
1106: MatGetBlockSize(Amat, &bs);
1107: nloc = (Iend-Istart)/bs; my0 = Istart/bs;
1108: if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) not divisible by bs %D",Iend,Istart,bs);
1110: /* get 'nLocalSelected' */
1111: for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1112: PetscBool ise;
1113: /* filter out singletons 0 or 1? */
1114: PetscCDEmptyAt(agg_lists, ii, &ise);
1115: if (!ise) nLocalSelected++;
1116: }
1118: /* create prolongator, create P matrix */
1119: MatCreate(comm, &Prol);
1120: MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);
1121: MatSetBlockSizes(Prol, bs, col_bs);
1122: MatSetType(Prol, MATAIJ);
1123: MatSeqAIJSetPreallocation(Prol, data_cols, NULL);
1124: MatMPIAIJSetPreallocation(Prol,data_cols, NULL,data_cols, NULL);
1125: /* nloc*bs, nLocalSelected*col_bs, */
1126: /* PETSC_DETERMINE, PETSC_DETERMINE, */
1127: /* data_cols, NULL, data_cols, NULL, */
1128: /* &Prol); */
1130: /* can get all points "removed" */
1131: MatGetSize(Prol, &kk, &ii);
1132: if (ii==0) {
1133: if (verbose > 0) PetscPrintf(comm,"[%d]%s no selected points on coarse grid\n",rank,__FUNCT__);
1134: MatDestroy(&Prol);
1135: *a_P_out = NULL; /* out */
1136: return(0);
1137: }
1138: if (verbose > 0) PetscPrintf(comm,"\t\t[%d]%s New grid %d nodes\n",rank,__FUNCT__,ii/col_bs);
1139: MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);
1141: 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);
1142: myCrs0 = myCrs0/col_bs;
1143: 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);
1145: /* create global vector of data in 'data_w_ghost' */
1146: #if defined PETSC_GAMG_USE_LOG
1147: PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);
1148: #endif
1149: if (size > 1) { /* */
1150: PetscReal *tmp_gdata,*tmp_ldata,*tp2;
1151: PetscMalloc(nloc*sizeof(PetscReal), &tmp_ldata);
1152: for (jj = 0; jj < data_cols; jj++) {
1153: for (kk = 0; kk < bs; kk++) {
1154: PetscInt ii,stride;
1155: const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk;
1156: for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1158: PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata);
1160: if (jj==0 && kk==0) { /* now I know how many todal nodes - allocate */
1161: PetscMalloc(stride*bs*data_cols*sizeof(PetscReal), &data_w_ghost);
1162: nbnodes = bs*stride;
1163: }
1164: tp2 = data_w_ghost + jj*bs*stride + kk;
1165: for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1166: PetscFree(tmp_gdata);
1167: }
1168: }
1169: PetscFree(tmp_ldata);
1170: } else {
1171: nbnodes = bs*nloc;
1172: data_w_ghost = (PetscReal*)pc_gamg->data;
1173: }
1175: /* get P0 */
1176: if (size > 1) {
1177: PetscReal *fid_glid_loc,*fiddata;
1178: PetscInt stride;
1180: PetscMalloc(nloc*sizeof(PetscReal), &fid_glid_loc);
1181: for (kk=0; kk<nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk);
1182: PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);
1183: PetscMalloc(stride*sizeof(PetscInt), &flid_fgid);
1184: for (kk=0; kk<stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1185: PetscFree(fiddata);
1187: if (stride != nbnodes/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"stride %D != nbnodes %D/bs %D",stride,nbnodes,bs);
1188: PetscFree(fid_glid_loc);
1189: } else {
1190: PetscMalloc(nloc*sizeof(PetscInt), &flid_fgid);
1191: for (kk=0; kk<nloc; kk++) flid_fgid[kk] = my0 + kk;
1192: }
1193: #if defined PETSC_GAMG_USE_LOG
1194: PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);
1195: PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);
1196: #endif
1197: {
1198: PetscReal *data_out = NULL;
1199: formProl0(agg_lists, bs, data_cols, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);
1200: PetscFree(pc_gamg->data);
1202: pc_gamg->data = data_out;
1203: pc_gamg->data_cell_rows = data_cols;
1204: pc_gamg->data_sz = data_cols*data_cols*nLocalSelected;
1205: }
1206: #if defined PETSC_GAMG_USE_LOG
1207: PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);
1208: #endif
1209: if (size > 1) PetscFree(data_w_ghost);
1210: PetscFree(flid_fgid);
1212: *a_P_out = Prol; /* out */
1213: #if defined PETSC_USE_LOG
1214: PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1215: #endif
1216: return(0);
1217: }
1219: /* -------------------------------------------------------------------------- */
1220: /*
1221: PCGAMGOptprol_AGG
1223: Input Parameter:
1224: . pc - this
1225: . Amat - matrix on this fine level
1226: In/Output Parameter:
1227: . a_P_out - prolongation operator to the next level
1228: */
1231: PetscErrorCode PCGAMGOptprol_AGG(PC pc,const Mat Amat,Mat *a_P)
1232: {
1234: PC_MG *mg = (PC_MG*)pc->data;
1235: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1236: const PetscInt verbose = pc_gamg->verbose;
1237: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1238: PetscInt jj;
1239: PetscMPIInt rank,size;
1240: Mat Prol = *a_P;
1241: MPI_Comm comm;
1244: PetscObjectGetComm((PetscObject)Amat,&comm);
1245: #if defined PETSC_USE_LOG
1246: PetscLogEventBegin(PC_GAMGOptprol_AGG,0,0,0,0);
1247: #endif
1248: MPI_Comm_rank(comm, &rank);
1249: MPI_Comm_size(comm, &size);
1251: /* smooth P0 */
1252: for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1253: Mat tMat;
1254: Vec diag;
1255: PetscReal alpha, emax, emin;
1256: #if defined PETSC_GAMG_USE_LOG
1257: PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);
1258: #endif
1259: if (jj == 0) {
1260: KSP eksp;
1261: Vec bb, xx;
1262: PC epc;
1263: MatGetVecs(Amat, &bb, 0);
1264: MatGetVecs(Amat, &xx, 0);
1265: {
1266: PetscRandom rctx;
1267: PetscRandomCreate(comm,&rctx);
1268: PetscRandomSetFromOptions(rctx);
1269: VecSetRandom(bb,rctx);
1270: PetscRandomDestroy(&rctx);
1271: }
1273: /* zeroing out BC rows -- needed for crazy matrices */
1274: {
1275: PetscInt Istart,Iend,ncols,jj,Ii;
1276: PetscScalar zero = 0.0;
1277: MatGetOwnershipRange(Amat, &Istart, &Iend);
1278: for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
1279: MatGetRow(Amat,Ii,&ncols,0,0);
1280: if (ncols <= 1) {
1281: VecSetValues(bb, 1, &Ii, &zero, INSERT_VALUES);
1282: }
1283: MatRestoreRow(Amat,Ii,&ncols,0,0);
1284: }
1285: VecAssemblyBegin(bb);
1286: VecAssemblyEnd(bb);
1287: }
1289: KSPCreate(comm,&eksp);
1290: KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
1291: KSPSetNormType(eksp, KSP_NORM_NONE);
1292: KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);
1293: KSPAppendOptionsPrefix(eksp, "gamg_est_");
1294: KSPSetFromOptions(eksp);
1296: KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
1297: KSPSetOperators(eksp, Amat, Amat, SAME_NONZERO_PATTERN);
1298: KSPSetComputeSingularValues(eksp,PETSC_TRUE);
1300: KSPGetPC(eksp, &epc);
1301: PCSetType(epc, PCJACOBI); /* smoother in smoothed agg. */
1303: /* solve - keep stuff out of logging */
1304: PetscLogEventDeactivate(KSP_Solve);
1305: PetscLogEventDeactivate(PC_Apply);
1306: KSPSolve(eksp, bb, xx);
1307: PetscLogEventActivate(KSP_Solve);
1308: PetscLogEventActivate(PC_Apply);
1310: KSPComputeExtremeSingularValues(eksp, &emax, &emin);
1311: if (verbose > 0) {
1312: PetscPrintf(comm,"\t\t\t%s smooth P0: max eigen=%e min=%e PC=%s\n",
1313: __FUNCT__,emax,emin,PCJACOBI);
1314: }
1315: VecDestroy(&xx);
1316: VecDestroy(&bb);
1317: KSPDestroy(&eksp);
1319: if (pc_gamg->emax_id == -1) {
1320: PetscObjectComposedDataRegister(&pc_gamg->emax_id);
1321: if (pc_gamg->emax_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->emax_id == -1");
1322: }
1323: PetscObjectComposedDataSetScalar((PetscObject)Amat, pc_gamg->emax_id, emax);
1324: }
1326: /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1327: MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);
1328: MatGetVecs(Amat, &diag, 0);
1329: MatGetDiagonal(Amat, diag); /* effectively PCJACOBI */
1330: VecReciprocal(diag);
1331: MatDiagonalScale(tMat, diag, 0);
1332: VecDestroy(&diag);
1333: alpha = -1.4/emax;
1334: MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);
1335: MatDestroy(&Prol);
1336: Prol = tMat;
1337: #if defined PETSC_GAMG_USE_LOG
1338: PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);
1339: #endif
1340: }
1341: #if defined PETSC_USE_LOG
1342: PetscLogEventEnd(PC_GAMGOptprol_AGG,0,0,0,0);
1343: #endif
1344: *a_P = Prol;
1345: return(0);
1346: }
1348: /* -------------------------------------------------------------------------- */
1349: /*
1350: PCGAMGKKTProl_AGG
1352: Input Parameter:
1353: . pc - this
1354: . Prol11 - matrix on this fine level
1355: . A21 - matrix on this fine level
1356: In/Output Parameter:
1357: . a_P22 - prolongation operator to the next level
1358: */
1361: PetscErrorCode PCGAMGKKTProl_AGG(PC pc,const Mat Prol11,const Mat A21,Mat *a_P22)
1362: {
1363: PetscErrorCode ierr;
1364: PC_MG *mg = (PC_MG*)pc->data;
1365: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1366: const PetscInt verbose = pc_gamg->verbose;
1367: /* PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; */
1368: PetscMPIInt rank,size;
1369: Mat Prol22,Tmat,Gmat;
1370: MPI_Comm comm;
1371: PetscCoarsenData *agg_lists;
1374: PetscObjectGetComm((PetscObject)pc,&comm);
1375: #if defined PETSC_USE_LOG
1376: PetscLogEventBegin(PC_GAMGKKTProl_AGG,0,0,0,0);
1377: #endif
1378: MPI_Comm_rank(comm, &rank);
1379: MPI_Comm_size(comm, &size);
1381: /* form C graph */
1382: MatMatMult(A21, Prol11, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Tmat);
1383: MatMatTransposeMult(Tmat, Tmat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat);
1384: MatDestroy(&Tmat);
1385: PCGAMGFilterGraph(&Gmat, 0.0, PETSC_FALSE, verbose);
1387: /* coarsen constraints */
1388: {
1389: MatCoarsen crs;
1390: MatCoarsenCreate(comm, &crs);
1391: MatCoarsenSetType(crs, MATCOARSENMIS);
1392: MatCoarsenSetAdjacency(crs, Gmat);
1393: MatCoarsenSetVerbose(crs, verbose);
1394: MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
1395: MatCoarsenApply(crs);
1396: MatCoarsenGetData(crs, &agg_lists);
1397: MatCoarsenDestroy(&crs);
1398: }
1400: /* form simple prolongation 'Prol22' */
1401: {
1402: PetscInt ii,mm,clid,my0,nloc,nLocalSelected;
1403: PetscScalar val = 1.0;
1404: /* get 'nLocalSelected' */
1405: MatGetLocalSize(Gmat, &nloc, &ii);
1406: for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1407: PetscBool ise;
1408: /* filter out singletons 0 or 1? */
1409: PetscCDEmptyAt(agg_lists, ii, &ise);
1410: if (!ise) nLocalSelected++;
1411: }
1413: MatCreate(comm,&Prol22);
1414: MatSetSizes(Prol22,nloc, nLocalSelected,PETSC_DETERMINE, PETSC_DETERMINE);
1415: MatSetType(Prol22, MATAIJ);
1416: MatSeqAIJSetPreallocation(Prol22,1,NULL);
1417: MatMPIAIJSetPreallocation(Prol22,1,NULL,1,NULL);
1418: /* MatCreateAIJ(comm, */
1419: /* nloc, nLocalSelected, */
1420: /* PETSC_DETERMINE, PETSC_DETERMINE, */
1421: /* 1, NULL, 1, NULL, */
1422: /* &Prol22); */
1424: MatGetOwnershipRange(Prol22, &my0, &ii);
1425: nloc = ii - my0;
1427: /* make aggregates */
1428: for (mm = clid = 0; mm < nloc; mm++) {
1429: PetscCDSizeAt(agg_lists, mm, &ii);
1430: if (ii > 0) {
1431: PetscInt asz=ii,cgid=my0+clid,rids[1000];
1432: PetscCDPos pos;
1433: if (asz>1000) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Very large aggregate: %d",asz);
1434: ii = 0;
1435: PetscCDGetHeadPos(agg_lists,mm,&pos);
1436: while (pos) {
1437: PetscInt gid1;
1438: PetscLLNGetID(pos, &gid1);
1439: PetscCDGetNextPos(agg_lists,mm,&pos);
1441: rids[ii++] = gid1;
1442: }
1443: if (ii != asz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D != asz %D",ii,asz);
1444: /* add diagonal block of P0 */
1445: MatSetValues(Prol22,asz,rids,1,&cgid,&val,INSERT_VALUES);
1447: clid++;
1448: } /* coarse agg */
1449: } /* for all fine nodes */
1450: MatAssemblyBegin(Prol22,MAT_FINAL_ASSEMBLY);
1451: MatAssemblyEnd(Prol22,MAT_FINAL_ASSEMBLY);
1452: }
1454: /* clean up */
1455: MatDestroy(&Gmat);
1456: PetscCDDestroy(agg_lists);
1457: #if defined PETSC_USE_LOG
1458: PetscLogEventEnd(PC_GAMGKKTProl_AGG,0,0,0,0);
1459: #endif
1460: *a_P22 = Prol22;
1461: return(0);
1462: }
1464: /* -------------------------------------------------------------------------- */
1465: /*
1466: PCCreateGAMG_AGG
1468: Input Parameter:
1469: . pc -
1470: */
1473: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1474: {
1476: PC_MG *mg = (PC_MG*)pc->data;
1477: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1478: PC_GAMG_AGG *pc_gamg_agg;
1481: /* create sub context for SA */
1482: PetscNewLog(pc, PC_GAMG_AGG, &pc_gamg_agg);
1483: pc_gamg->subctx = pc_gamg_agg;
1485: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1486: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
1487: /* reset does not do anything; setup not virtual */
1489: /* set internal function pointers */
1490: pc_gamg->ops->graph = PCGAMGgraph_AGG;
1491: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
1492: pc_gamg->ops->prolongator = PCGAMGProlongator_AGG;
1493: pc_gamg->ops->optprol = PCGAMGOptprol_AGG;
1494: pc_gamg->ops->formkktprol = PCGAMGKKTProl_AGG;
1496: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1498: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);
1499: return(0);
1500: }