Actual source code: agg.c
petsc-3.4.2 2013-07-02
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: PetscMalloc((nvec+!!has_const)*mlocal*sizeof(*nullvec),&nullvec);
680: if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
681: for (i=0; i<nvec; i++) {
682: VecGetArrayRead(vecs[i],&v);
683: for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]);
684: VecRestoreArrayRead(vecs[i],&v);
685: }
686: pc_gamg->data = nullvec;
687: pc_gamg->data_cell_cols = (nvec+!!has_const);
689: MatGetBlockSize(a_A, &bs); /* this does not work for Stokes */
691: pc_gamg->data_cell_rows = bs;
692: }
693: return(0);
694: }
696: /* -------------------------------------------------------------------------- */
697: /*
698: formProl0
700: Input Parameter:
701: . agg_llists - list of arrays with aggregates
702: . bs - block size
703: . nSAvec - column bs of new P
704: . my0crs - global index of start of locals
705: . data_stride - bs*(nloc nodes + ghost nodes)
706: . data_in[data_stride*nSAvec] - local data on fine grid
707: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
708: Output Parameter:
709: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
710: . a_Prol - prolongation operator
711: */
714: static PetscErrorCode formProl0(const PetscCoarsenData *agg_llists, /* list from selected vertices of aggregate unselected vertices */
715: const PetscInt bs, /* (row) block size */
716: const PetscInt nSAvec, /* column bs */
717: const PetscInt my0crs, /* global index of start of locals */
718: const PetscInt data_stride, /* (nloc+nghost)*bs */
719: PetscReal data_in[], /* [data_stride][nSAvec] */
720: const PetscInt flid_fgid[], /* [data_stride/bs] */
721: PetscReal **a_data_out,
722: Mat a_Prol) /* prolongation operator (output)*/
723: {
725: PetscInt Istart,my0,Iend,nloc,clid,flid,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride;
726: MPI_Comm comm;
727: PetscMPIInt rank, size;
728: PetscReal *out_data;
729: PetscCDPos pos;
730: GAMGHashTable fgid_flid;
732: /* #define OUT_AGGS */
733: #if defined(OUT_AGGS)
734: static PetscInt llev = 0; char fname[32]; FILE *file = NULL; PetscInt pM;
735: #endif
738: PetscObjectGetComm((PetscObject)a_Prol,&comm);
739: MPI_Comm_rank(comm,&rank);
740: MPI_Comm_size(comm,&size);
741: MatGetOwnershipRange(a_Prol, &Istart, &Iend);
742: nloc = (Iend-Istart)/bs; my0 = Istart/bs;
743: if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Iend %D - Istart %d must be divisible by bs %D",Iend,Istart,bs);
744: Iend /= bs;
745: nghosts = data_stride/bs - nloc;
747: GAMGTableCreate(2*nghosts+1, &fgid_flid);
748: for (kk=0; kk<nghosts; kk++) {
749: GAMGTableAdd(&fgid_flid, flid_fgid[nloc+kk], nloc+kk);
750: }
752: #if defined(OUT_AGGS)
753: sprintf(fname,"aggs_%d_%d.m",llev++,rank);
754: if (llev==1) file = fopen(fname,"w");
755: MatGetSize(a_Prol, &pM, &jj);
756: #endif
758: /* count selected -- same as number of cols of P */
759: for (nSelected=mm=0; mm<nloc; mm++) {
760: PetscBool ise;
761: PetscCDEmptyAt(agg_llists, mm, &ise);
762: if (!ise) nSelected++;
763: }
764: MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);
765: if ((ii/nSAvec) != my0crs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D /nSAvec %D != my0crs %D",ii,nSAvec,my0crs);
766: if (nSelected != (jj-ii)/nSAvec) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nSelected %D != (jj %D - ii %D)/nSAvec %D",nSelected,jj,ii,nSAvec);
768: /* aloc space for coarse point data (output) */
769: out_data_stride = nSelected*nSAvec;
771: PetscMalloc(out_data_stride*nSAvec*sizeof(PetscReal), &out_data);
772: for (ii=0;ii<out_data_stride*nSAvec;ii++) out_data[ii]=PETSC_MAX_REAL;
773: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
775: /* find points and set prolongation */
776: minsz = 100;
777: ndone = 0;
778: for (mm = clid = 0; mm < nloc; mm++) {
779: PetscCDSizeAt(agg_llists, mm, &jj);
780: if (jj > 0) {
781: const PetscInt lid = mm, cgid = my0crs + clid;
782: PetscInt cids[100]; /* max bs */
783: PetscBLASInt asz =jj,M=asz*bs,N=nSAvec,INFO;
784: PetscBLASInt Mdata=M+((N-M>0) ? N-M : 0),LDA=Mdata,LWORK=N*bs;
785: PetscScalar *qqc,*qqr,*TAU,*WORK;
786: PetscInt *fids;
787: PetscReal *data;
788: /* count agg */
789: if (asz<minsz) minsz = asz;
791: /* get block */
792: PetscMalloc((Mdata*N)*sizeof(PetscScalar), &qqc);
793: PetscMalloc((M*N)*sizeof(PetscScalar), &qqr);
794: PetscMalloc(N*sizeof(PetscScalar), &TAU);
795: PetscMalloc(LWORK*sizeof(PetscScalar), &WORK);
796: PetscMalloc(M*sizeof(PetscInt), &fids);
798: aggID = 0;
799: PetscCDGetHeadPos(agg_llists,lid,&pos);
800: while (pos) {
801: PetscInt gid1;
802: PetscLLNGetID(pos, &gid1);
803: PetscCDGetNextPos(agg_llists,lid,&pos);
805: if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
806: else {
807: GAMGTableFind(&fgid_flid, gid1, &flid);
808: if (flid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot find gid1 in table");
809: }
810: /* copy in B_i matrix - column oriented */
811: data = &data_in[flid*bs];
812: for (kk = ii = 0; ii < bs; ii++) {
813: for (jj = 0; jj < N; jj++) {
814: PetscReal d = data[jj*data_stride + ii];
815: qqc[jj*Mdata + aggID*bs + ii] = d;
816: }
817: }
818: #if defined(OUT_AGGS)
819: if (llev==1) {
820: char str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^";
821: PetscInt MM,pi,pj;
822: str[12] = col[(clid+17*rank)%6]; str[13] = sim[(clid+17*rank)%11];
823: MM = (PetscInt)(PetscSqrtReal((PetscReal)pM));
824: pj = gid1/MM; pi = gid1%MM;
825: fprintf(file,str,(double)pi,(double)pj);
826: /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */
827: }
828: #endif
829: /* set fine IDs */
830: for (kk=0; kk<bs; kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk;
832: aggID++;
833: }
835: /* pad with zeros */
836: for (ii = asz*bs; ii < Mdata; ii++) {
837: for (jj = 0; jj < N; jj++, kk++) {
838: qqc[jj*Mdata + ii] = .0;
839: }
840: }
842: ndone += aggID;
843: /* QR */
844: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
845: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
846: PetscFPTrapPop();
847: if (INFO != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xGEQRF error");
848: /* get R - column oriented - output B_{i+1} */
849: {
850: PetscReal *data = &out_data[clid*nSAvec];
851: for (jj = 0; jj < nSAvec; jj++) {
852: for (ii = 0; ii < nSAvec; ii++) {
853: 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);
854: if (ii <= jj) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]);
855: else data[jj*out_data_stride + ii] = 0.;
856: }
857: }
858: }
860: /* get Q - row oriented */
861: PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
862: if (INFO != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xORGQR error arg %d",-INFO);
864: for (ii = 0; ii < M; ii++) {
865: for (jj = 0; jj < N; jj++) {
866: qqr[N*ii + jj] = qqc[jj*Mdata + ii];
867: }
868: }
870: /* add diagonal block of P0 */
871: for (kk=0; kk<N; kk++) {
872: cids[kk] = N*cgid + kk; /* global col IDs in P0 */
873: }
874: MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);
876: PetscFree(qqc);
877: PetscFree(qqr);
878: PetscFree(TAU);
879: PetscFree(WORK);
880: PetscFree(fids);
881: clid++;
882: } /* coarse agg */
883: } /* for all fine nodes */
884: MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
885: MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);
887: /* MPI_Allreduce(&ndone, &ii, 1, MPIU_INT, MPIU_SUM, comm); */
888: /* MatGetSize(a_Prol, &kk, &jj); */
889: /* MPI_Allreduce(&minsz, &jj, 1, MPIU_INT, MPIU_MIN, comm); */
890: /* PetscPrintf(comm," **** [%d]%s %d total done, %d nodes (%d local done), min agg. size = %d\n",rank,__FUNCT__,ii,kk/bs,ndone,jj); */
892: #if defined(OUT_AGGS)
893: if (llev==1) fclose(file);
894: #endif
895: GAMGTableDestroy(&fgid_flid);
896: return(0);
897: }
899: /* -------------------------------------------------------------------------- */
900: /*
901: PCGAMGgraph_AGG
903: Input Parameter:
904: . pc - this
905: . Amat - matrix on this fine level
906: Output Parameter:
907: . a_Gmat -
908: */
911: PetscErrorCode PCGAMGgraph_AGG(PC pc,const Mat Amat,Mat *a_Gmat)
912: {
913: PetscErrorCode ierr;
914: PC_MG *mg = (PC_MG*)pc->data;
915: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
916: const PetscInt verbose = pc_gamg->verbose;
917: const PetscReal vfilter = pc_gamg->threshold;
918: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
919: PetscMPIInt rank,size;
920: Mat Gmat;
921: MPI_Comm comm;
922: PetscBool /* set,flg , */ symm;
925: PetscObjectGetComm((PetscObject)Amat,&comm);
926: #if defined PETSC_USE_LOG
927: PetscLogEventBegin(PC_GAMGGgraph_AGG,0,0,0,0);
928: #endif
929: MPI_Comm_rank(comm, &rank);
930: MPI_Comm_size(comm, &size);
932: /* MatIsSymmetricKnown(Amat, &set, &flg); || !(set && flg) -- this causes lot of symm calls */
933: symm = (PetscBool)(pc_gamg_agg->sym_graph); /* && !pc_gamg_agg->square_graph; */
935: PCGAMGCreateGraph(Amat, &Gmat);
936: PCGAMGFilterGraph(&Gmat, vfilter, symm, verbose);
938: *a_Gmat = Gmat;
940: #if defined PETSC_USE_LOG
941: PetscLogEventEnd(PC_GAMGGgraph_AGG,0,0,0,0);
942: #endif
943: return(0);
944: }
946: /* -------------------------------------------------------------------------- */
947: /*
948: PCGAMGCoarsen_AGG
950: Input Parameter:
951: . a_pc - this
952: Input/Output Parameter:
953: . a_Gmat1 - graph on this fine level - coarsening can change this (squares it)
954: Output Parameter:
955: . agg_lists - list of aggregates
956: */
959: PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc,Mat *a_Gmat1,PetscCoarsenData **agg_lists)
960: {
962: PC_MG *mg = (PC_MG*)a_pc->data;
963: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
964: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
965: Mat mat,Gmat2, Gmat1 = *a_Gmat1; /* squared graph */
966: const PetscInt verbose = pc_gamg->verbose;
967: IS perm;
968: PetscInt Ii,nloc,bs,n,m;
969: PetscInt *permute;
970: PetscBool *bIndexSet;
971: MatCoarsen crs;
972: MPI_Comm comm;
973: PetscMPIInt rank,size;
976: #if defined PETSC_USE_LOG
977: PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);
978: #endif
979: PetscObjectGetComm((PetscObject)Gmat1,&comm);
980: MPI_Comm_rank(comm, &rank);
981: MPI_Comm_size(comm, &size);
982: MatGetLocalSize(Gmat1, &n, &m);
983: MatGetBlockSize(Gmat1, &bs);
984: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bs %D must be 1",bs);
985: nloc = n/bs;
987: if (pc_gamg_agg->square_graph) {
988: if (verbose > 1) PetscPrintf(comm,"[%d]%s square graph\n",rank,__FUNCT__);
989: /* MatMatTransposeMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2); */
990: MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);
991: if (verbose > 2) {
992: PetscPrintf(comm,"[%d]%s square graph done\n",rank,__FUNCT__);
993: /* check for symetry */
994: if (verbose > 4) {
996: }
997: }
998: } else Gmat2 = Gmat1;
1000: /* get MIS aggs */
1001: /* randomize */
1002: PetscMalloc(nloc*sizeof(PetscInt), &permute);
1003: PetscMalloc(nloc*sizeof(PetscBool), &bIndexSet);
1004: for (Ii = 0; Ii < nloc; Ii++) {
1005: bIndexSet[Ii] = PETSC_FALSE;
1006: permute[Ii] = Ii;
1007: }
1008: srand(1); /* make deterministic */
1009: for (Ii = 0; Ii < nloc; Ii++) {
1010: PetscInt iSwapIndex = rand()%nloc;
1011: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1012: PetscInt iTemp = permute[iSwapIndex];
1013: permute[iSwapIndex] = permute[Ii];
1014: permute[Ii] = iTemp;
1015: bIndexSet[iSwapIndex] = PETSC_TRUE;
1016: }
1017: }
1018: PetscFree(bIndexSet);
1020: if (verbose > 1) PetscPrintf(comm,"[%d]%s coarsen graph\n",rank,__FUNCT__);
1022: ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);
1023: #if defined PETSC_GAMG_USE_LOG
1024: PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);
1025: #endif
1026: MatCoarsenCreate(comm, &crs);
1027: /* MatCoarsenSetType(crs, MATCOARSENMIS); */
1028: MatCoarsenSetFromOptions(crs);
1029: MatCoarsenSetGreedyOrdering(crs, perm);
1030: MatCoarsenSetAdjacency(crs, Gmat2);
1031: MatCoarsenSetVerbose(crs, pc_gamg->verbose);
1032: MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
1033: MatCoarsenApply(crs);
1034: MatCoarsenGetData(crs, agg_lists); /* output */
1035: MatCoarsenDestroy(&crs);
1037: ISDestroy(&perm);
1038: PetscFree(permute);
1039: #if defined PETSC_GAMG_USE_LOG
1040: PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);
1041: #endif
1042: if (verbose > 2) PetscPrintf(comm,"[%d]%s coarsen graph done\n",rank,__FUNCT__);
1044: /* smooth aggs */
1045: if (Gmat2 != Gmat1) {
1046: const PetscCoarsenData *llist = *agg_lists;
1047: smoothAggs(Gmat2, Gmat1, *agg_lists);
1048: MatDestroy(&Gmat1);
1049: *a_Gmat1 = Gmat2; /* output */
1050: PetscCDGetMat(llist, &mat);
1051: if (mat) SETERRQ(comm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????");
1052: } else {
1053: const PetscCoarsenData *llist = *agg_lists;
1054: /* see if we have a matrix that takes precedence (returned from MatCoarsenAppply) */
1055: PetscCDGetMat(llist, &mat);
1056: if (mat) {
1057: MatDestroy(&Gmat1);
1058: *a_Gmat1 = mat; /* output */
1059: }
1060: }
1061: #if defined PETSC_USE_LOG
1062: PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);
1063: #endif
1064: if (verbose > 2) PetscPrintf(comm,"[%d]%s PCGAMGCoarsen_AGG done\n",rank,__FUNCT__);
1065: return(0);
1066: }
1068: /* -------------------------------------------------------------------------- */
1069: /*
1070: PCGAMGProlongator_AGG
1072: Input Parameter:
1073: . pc - this
1074: . Amat - matrix on this fine level
1075: . Graph - used to get ghost data for nodes in
1076: . agg_lists - list of aggregates
1077: Output Parameter:
1078: . a_P_out - prolongation operator to the next level
1079: */
1082: PetscErrorCode PCGAMGProlongator_AGG(PC pc,const Mat Amat,const Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
1083: {
1084: PC_MG *mg = (PC_MG*)pc->data;
1085: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1086: const PetscInt verbose = pc_gamg->verbose;
1087: const PetscInt data_cols = pc_gamg->data_cell_cols;
1089: PetscInt Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs;
1090: Mat Prol;
1091: PetscMPIInt rank, size;
1092: MPI_Comm comm;
1093: const PetscInt col_bs = data_cols;
1094: PetscReal *data_w_ghost;
1095: PetscInt myCrs0, nbnodes=0, *flid_fgid;
1098: PetscObjectGetComm((PetscObject)Amat,&comm);
1099: #if defined PETSC_USE_LOG
1100: PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);
1101: #endif
1102: MPI_Comm_rank(comm, &rank);
1103: MPI_Comm_size(comm, &size);
1104: MatGetOwnershipRange(Amat, &Istart, &Iend);
1105: MatGetBlockSize(Amat, &bs);
1106: nloc = (Iend-Istart)/bs; my0 = Istart/bs;
1107: if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) not divisible by bs %D",Iend,Istart,bs);
1109: /* get 'nLocalSelected' */
1110: for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1111: PetscBool ise;
1112: /* filter out singletons 0 or 1? */
1113: PetscCDEmptyAt(agg_lists, ii, &ise);
1114: if (!ise) nLocalSelected++;
1115: }
1117: /* create prolongator, create P matrix */
1118: MatCreate(comm, &Prol);
1119: MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);
1120: MatSetBlockSizes(Prol, bs, col_bs);
1121: MatSetType(Prol, MATAIJ);
1122: MatSeqAIJSetPreallocation(Prol, data_cols, NULL);
1123: MatMPIAIJSetPreallocation(Prol,data_cols, NULL,data_cols, NULL);
1124: /* nloc*bs, nLocalSelected*col_bs, */
1125: /* PETSC_DETERMINE, PETSC_DETERMINE, */
1126: /* data_cols, NULL, data_cols, NULL, */
1127: /* &Prol); */
1129: /* can get all points "removed" */
1130: MatGetSize(Prol, &kk, &ii);
1131: if (ii==0) {
1132: if (verbose > 0) PetscPrintf(comm,"[%d]%s no selected points on coarse grid\n",rank,__FUNCT__);
1133: MatDestroy(&Prol);
1134: *a_P_out = NULL; /* out */
1135: return(0);
1136: }
1137: if (verbose > 0) PetscPrintf(comm,"\t\t[%d]%s New grid %d nodes\n",rank,__FUNCT__,ii/col_bs);
1138: MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);
1140: 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);
1141: myCrs0 = myCrs0/col_bs;
1142: 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);
1144: /* create global vector of data in 'data_w_ghost' */
1145: #if defined PETSC_GAMG_USE_LOG
1146: PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);
1147: #endif
1148: if (size > 1) { /* */
1149: PetscReal *tmp_gdata,*tmp_ldata,*tp2;
1150: PetscMalloc(nloc*sizeof(PetscReal), &tmp_ldata);
1151: for (jj = 0; jj < data_cols; jj++) {
1152: for (kk = 0; kk < bs; kk++) {
1153: PetscInt ii,stride;
1154: const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk;
1155: for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;
1157: PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata);
1159: if (jj==0 && kk==0) { /* now I know how many todal nodes - allocate */
1160: PetscMalloc(stride*bs*data_cols*sizeof(PetscReal), &data_w_ghost);
1161: nbnodes = bs*stride;
1162: }
1163: tp2 = data_w_ghost + jj*bs*stride + kk;
1164: for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1165: PetscFree(tmp_gdata);
1166: }
1167: }
1168: PetscFree(tmp_ldata);
1169: } else {
1170: nbnodes = bs*nloc;
1171: data_w_ghost = (PetscReal*)pc_gamg->data;
1172: }
1174: /* get P0 */
1175: if (size > 1) {
1176: PetscReal *fid_glid_loc,*fiddata;
1177: PetscInt stride;
1179: PetscMalloc(nloc*sizeof(PetscReal), &fid_glid_loc);
1180: for (kk=0; kk<nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk);
1181: PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);
1182: PetscMalloc(stride*sizeof(PetscInt), &flid_fgid);
1183: for (kk=0; kk<stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1184: PetscFree(fiddata);
1186: if (stride != nbnodes/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"stride %D != nbnodes %D/bs %D",stride,nbnodes,bs);
1187: PetscFree(fid_glid_loc);
1188: } else {
1189: PetscMalloc(nloc*sizeof(PetscInt), &flid_fgid);
1190: for (kk=0; kk<nloc; kk++) flid_fgid[kk] = my0 + kk;
1191: }
1192: #if defined PETSC_GAMG_USE_LOG
1193: PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);
1194: PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);
1195: #endif
1196: {
1197: PetscReal *data_out = NULL;
1198: formProl0(agg_lists, bs, data_cols, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);
1199: PetscFree(pc_gamg->data);
1201: pc_gamg->data = data_out;
1202: pc_gamg->data_cell_rows = data_cols;
1203: pc_gamg->data_sz = data_cols*data_cols*nLocalSelected;
1204: }
1205: #if defined PETSC_GAMG_USE_LOG
1206: PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);
1207: #endif
1208: if (size > 1) PetscFree(data_w_ghost);
1209: PetscFree(flid_fgid);
1211: *a_P_out = Prol; /* out */
1212: #if defined PETSC_USE_LOG
1213: PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1214: #endif
1215: return(0);
1216: }
1218: /* -------------------------------------------------------------------------- */
1219: /*
1220: PCGAMGOptprol_AGG
1222: Input Parameter:
1223: . pc - this
1224: . Amat - matrix on this fine level
1225: In/Output Parameter:
1226: . a_P_out - prolongation operator to the next level
1227: */
1230: PetscErrorCode PCGAMGOptprol_AGG(PC pc,const Mat Amat,Mat *a_P)
1231: {
1233: PC_MG *mg = (PC_MG*)pc->data;
1234: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1235: const PetscInt verbose = pc_gamg->verbose;
1236: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1237: PetscInt jj;
1238: PetscMPIInt rank,size;
1239: Mat Prol = *a_P;
1240: MPI_Comm comm;
1243: PetscObjectGetComm((PetscObject)Amat,&comm);
1244: #if defined PETSC_USE_LOG
1245: PetscLogEventBegin(PC_GAMGOptprol_AGG,0,0,0,0);
1246: #endif
1247: MPI_Comm_rank(comm, &rank);
1248: MPI_Comm_size(comm, &size);
1250: /* smooth P0 */
1251: for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1252: Mat tMat;
1253: Vec diag;
1254: PetscReal alpha, emax, emin;
1255: #if defined PETSC_GAMG_USE_LOG
1256: PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);
1257: #endif
1258: if (jj == 0) {
1259: KSP eksp;
1260: Vec bb, xx;
1261: PC epc;
1262: MatGetVecs(Amat, &bb, 0);
1263: MatGetVecs(Amat, &xx, 0);
1264: {
1265: PetscRandom rctx;
1266: PetscRandomCreate(comm,&rctx);
1267: PetscRandomSetFromOptions(rctx);
1268: VecSetRandom(bb,rctx);
1269: PetscRandomDestroy(&rctx);
1270: }
1272: /* zeroing out BC rows -- needed for crazy matrices */
1273: {
1274: PetscInt Istart,Iend,ncols,jj,Ii;
1275: PetscScalar zero = 0.0;
1276: MatGetOwnershipRange(Amat, &Istart, &Iend);
1277: for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
1278: MatGetRow(Amat,Ii,&ncols,0,0);
1279: if (ncols <= 1) {
1280: VecSetValues(bb, 1, &Ii, &zero, INSERT_VALUES);
1281: }
1282: MatRestoreRow(Amat,Ii,&ncols,0,0);
1283: }
1284: VecAssemblyBegin(bb);
1285: VecAssemblyEnd(bb);
1286: }
1288: KSPCreate(comm,&eksp);
1289: KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
1290: KSPSetNormType(eksp, KSP_NORM_NONE);
1291: KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);
1292: KSPAppendOptionsPrefix(eksp, "gamg_est_");
1293: KSPSetFromOptions(eksp);
1295: KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
1296: KSPSetOperators(eksp, Amat, Amat, SAME_NONZERO_PATTERN);
1297: KSPSetComputeSingularValues(eksp,PETSC_TRUE);
1299: KSPGetPC(eksp, &epc);
1300: PCSetType(epc, PCJACOBI); /* smoother in smoothed agg. */
1302: /* solve - keep stuff out of logging */
1303: PetscLogEventDeactivate(KSP_Solve);
1304: PetscLogEventDeactivate(PC_Apply);
1305: KSPSolve(eksp, bb, xx);
1306: PetscLogEventActivate(KSP_Solve);
1307: PetscLogEventActivate(PC_Apply);
1309: KSPComputeExtremeSingularValues(eksp, &emax, &emin);
1310: if (verbose > 0) {
1311: PetscPrintf(comm,"\t\t\t%s smooth P0: max eigen=%e min=%e PC=%s\n",
1312: __FUNCT__,emax,emin,PCJACOBI);
1313: }
1314: VecDestroy(&xx);
1315: VecDestroy(&bb);
1316: KSPDestroy(&eksp);
1318: if (pc_gamg->emax_id == -1) {
1319: PetscObjectComposedDataRegister(&pc_gamg->emax_id);
1320: if (pc_gamg->emax_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->emax_id == -1");
1321: }
1322: PetscObjectComposedDataSetScalar((PetscObject)Amat, pc_gamg->emax_id, emax);
1323: }
1325: /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1326: MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);
1327: MatGetVecs(Amat, &diag, 0);
1328: MatGetDiagonal(Amat, diag); /* effectively PCJACOBI */
1329: VecReciprocal(diag);
1330: MatDiagonalScale(tMat, diag, 0);
1331: VecDestroy(&diag);
1332: alpha = -1.4/emax;
1333: MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);
1334: MatDestroy(&Prol);
1335: Prol = tMat;
1336: #if defined PETSC_GAMG_USE_LOG
1337: PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);
1338: #endif
1339: }
1340: #if defined PETSC_USE_LOG
1341: PetscLogEventEnd(PC_GAMGOptprol_AGG,0,0,0,0);
1342: #endif
1343: *a_P = Prol;
1344: return(0);
1345: }
1347: /* -------------------------------------------------------------------------- */
1348: /*
1349: PCGAMGKKTProl_AGG
1351: Input Parameter:
1352: . pc - this
1353: . Prol11 - matrix on this fine level
1354: . A21 - matrix on this fine level
1355: In/Output Parameter:
1356: . a_P22 - prolongation operator to the next level
1357: */
1360: PetscErrorCode PCGAMGKKTProl_AGG(PC pc,const Mat Prol11,const Mat A21,Mat *a_P22)
1361: {
1362: PetscErrorCode ierr;
1363: PC_MG *mg = (PC_MG*)pc->data;
1364: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1365: const PetscInt verbose = pc_gamg->verbose;
1366: /* PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; */
1367: PetscMPIInt rank,size;
1368: Mat Prol22,Tmat,Gmat;
1369: MPI_Comm comm;
1370: PetscCoarsenData *agg_lists;
1373: PetscObjectGetComm((PetscObject)pc,&comm);
1374: #if defined PETSC_USE_LOG
1375: PetscLogEventBegin(PC_GAMGKKTProl_AGG,0,0,0,0);
1376: #endif
1377: MPI_Comm_rank(comm, &rank);
1378: MPI_Comm_size(comm, &size);
1380: /* form C graph */
1381: MatMatMult(A21, Prol11, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Tmat);
1382: MatMatTransposeMult(Tmat, Tmat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat);
1383: MatDestroy(&Tmat);
1384: PCGAMGFilterGraph(&Gmat, 0.0, PETSC_FALSE, verbose);
1386: /* coarsen constraints */
1387: {
1388: MatCoarsen crs;
1389: MatCoarsenCreate(comm, &crs);
1390: MatCoarsenSetType(crs, MATCOARSENMIS);
1391: MatCoarsenSetAdjacency(crs, Gmat);
1392: MatCoarsenSetVerbose(crs, verbose);
1393: MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
1394: MatCoarsenApply(crs);
1395: MatCoarsenGetData(crs, &agg_lists);
1396: MatCoarsenDestroy(&crs);
1397: }
1399: /* form simple prolongation 'Prol22' */
1400: {
1401: PetscInt ii,mm,clid,my0,nloc,nLocalSelected;
1402: PetscScalar val = 1.0;
1403: /* get 'nLocalSelected' */
1404: MatGetLocalSize(Gmat, &nloc, &ii);
1405: for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1406: PetscBool ise;
1407: /* filter out singletons 0 or 1? */
1408: PetscCDEmptyAt(agg_lists, ii, &ise);
1409: if (!ise) nLocalSelected++;
1410: }
1412: MatCreate(comm,&Prol22);
1413: MatSetSizes(Prol22,nloc, nLocalSelected,PETSC_DETERMINE, PETSC_DETERMINE);
1414: MatSetType(Prol22, MATAIJ);
1415: MatSeqAIJSetPreallocation(Prol22,1,NULL);
1416: MatMPIAIJSetPreallocation(Prol22,1,NULL,1,NULL);
1417: /* MatCreateAIJ(comm, */
1418: /* nloc, nLocalSelected, */
1419: /* PETSC_DETERMINE, PETSC_DETERMINE, */
1420: /* 1, NULL, 1, NULL, */
1421: /* &Prol22); */
1423: MatGetOwnershipRange(Prol22, &my0, &ii);
1424: nloc = ii - my0;
1426: /* make aggregates */
1427: for (mm = clid = 0; mm < nloc; mm++) {
1428: PetscCDSizeAt(agg_lists, mm, &ii);
1429: if (ii > 0) {
1430: PetscInt asz=ii,cgid=my0+clid,rids[1000];
1431: PetscCDPos pos;
1432: if (asz>1000) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Very large aggregate: %d",asz);
1433: ii = 0;
1434: PetscCDGetHeadPos(agg_lists,mm,&pos);
1435: while (pos) {
1436: PetscInt gid1;
1437: PetscLLNGetID(pos, &gid1);
1438: PetscCDGetNextPos(agg_lists,mm,&pos);
1440: rids[ii++] = gid1;
1441: }
1442: if (ii != asz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D != asz %D",ii,asz);
1443: /* add diagonal block of P0 */
1444: MatSetValues(Prol22,asz,rids,1,&cgid,&val,INSERT_VALUES);
1446: clid++;
1447: } /* coarse agg */
1448: } /* for all fine nodes */
1449: MatAssemblyBegin(Prol22,MAT_FINAL_ASSEMBLY);
1450: MatAssemblyEnd(Prol22,MAT_FINAL_ASSEMBLY);
1451: }
1453: /* clean up */
1454: MatDestroy(&Gmat);
1455: PetscCDDestroy(agg_lists);
1456: #if defined PETSC_USE_LOG
1457: PetscLogEventEnd(PC_GAMGKKTProl_AGG,0,0,0,0);
1458: #endif
1459: *a_P22 = Prol22;
1460: return(0);
1461: }
1463: /* -------------------------------------------------------------------------- */
1464: /*
1465: PCCreateGAMG_AGG
1467: Input Parameter:
1468: . pc -
1469: */
1472: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1473: {
1475: PC_MG *mg = (PC_MG*)pc->data;
1476: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1477: PC_GAMG_AGG *pc_gamg_agg;
1480: /* create sub context for SA */
1481: PetscNewLog(pc, PC_GAMG_AGG, &pc_gamg_agg);
1482: pc_gamg->subctx = pc_gamg_agg;
1484: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1485: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
1486: /* reset does not do anything; setup not virtual */
1488: /* set internal function pointers */
1489: pc_gamg->ops->graph = PCGAMGgraph_AGG;
1490: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
1491: pc_gamg->ops->prolongator = PCGAMGProlongator_AGG;
1492: pc_gamg->ops->optprol = PCGAMGOptprol_AGG;
1493: pc_gamg->ops->formkktprol = PCGAMGKKTProl_AGG;
1495: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1497: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);
1498: return(0);
1499: }