Actual source code: agg.c
petsc-3.5.4 2015-05-23
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: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
215: return(0);
216: }
218: /* -------------------------------------------------------------------------- */
219: /*
220: PCSetCoordinates_AGG
221: - collective
223: Input Parameter:
224: . pc - the preconditioner context
225: . ndm - dimesion of data (used for dof/vertex for Stokes)
226: . a_nloc - number of vertices local
227: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
228: */
232: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
233: {
234: PC_MG *mg = (PC_MG*)pc->data;
235: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
237: PetscInt arrsz,kk,ii,jj,nloc,ndatarows,ndf;
238: Mat mat = pc->pmat;
243: nloc = a_nloc;
245: /* SA: null space vectors */
246: MatGetBlockSize(mat, &ndf); /* this does not work for Stokes */
247: if (coords && ndf==1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
248: else if (coords) {
249: if (ndm > ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"degrees of motion %d > block size %d",ndm,ndf);
250: pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */
251: if (ndm != ndf) {
252: 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);
253: }
254: } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
255: pc_gamg->data_cell_rows = ndatarows = ndf;
256: 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);
257: arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;
259: /* create data - syntactic sugar that should be refactored at some point */
260: if (pc_gamg->data==0 || (pc_gamg->data_sz != arrsz)) {
261: PetscFree(pc_gamg->data);
262: PetscMalloc1((arrsz+1), &pc_gamg->data);
263: }
264: /* copy data in - column oriented */
265: for (kk=0; kk<nloc; kk++) {
266: const PetscInt M = nloc*pc_gamg->data_cell_rows; /* stride into data */
267: PetscReal *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */
268: if (pc_gamg->data_cell_cols==1) *data = 1.0;
269: else {
270: /* translational modes */
271: for (ii=0;ii<ndatarows;ii++) {
272: for (jj=0;jj<ndatarows;jj++) {
273: if (ii==jj)data[ii*M + jj] = 1.0;
274: else data[ii*M + jj] = 0.0;
275: }
276: }
278: /* rotational modes */
279: if (coords) {
280: if (ndm == 2) {
281: data += 2*M;
282: data[0] = -coords[2*kk+1];
283: data[1] = coords[2*kk];
284: } else {
285: data += 3*M;
286: data[0] = 0.0; data[M+0] = coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1];
287: data[1] = -coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = coords[3*kk];
288: data[2] = coords[3*kk+1]; data[M+2] = -coords[3*kk]; data[2*M+2] = 0.0;
289: }
290: }
291: }
292: }
294: pc_gamg->data_sz = arrsz;
295: return(0);
296: }
298: typedef PetscInt NState;
299: static const NState NOT_DONE=-2;
300: static const NState DELETED =-1;
301: static const NState REMOVED =-3;
302: #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED)
304: /* -------------------------------------------------------------------------- */
305: /*
306: smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific
307: - AGG-MG specific: clears singletons out of 'selected_2'
309: Input Parameter:
310: . Gmat_2 - glabal matrix of graph (data not defined)
311: . Gmat_1 - base graph to grab with
312: Input/Output Parameter:
313: . aggs_2 - linked list of aggs with gids)
314: */
317: static PetscErrorCode smoothAggs(const Mat Gmat_2, /* base (squared) graph */
318: const Mat Gmat_1, /* base graph */
319: /* const IS selected_2, [nselected local] selected vertices */
320: PetscCoarsenData *aggs_2) /* [nselected local] global ID of aggregate */
321: {
323: PetscBool isMPI;
324: Mat_SeqAIJ *matA_1, *matB_1=0, *matA_2, *matB_2=0;
325: MPI_Comm comm;
326: PetscMPIInt rank,size;
327: PetscInt lid,*ii,*idx,ix,Iend,my0,kk,n,j;
328: Mat_MPIAIJ *mpimat_2 = 0, *mpimat_1=0;
329: const PetscInt nloc = Gmat_2->rmap->n;
330: PetscScalar *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid;
331: PetscInt *lid_cprowID_1;
332: NState *lid_state;
333: Vec ghost_par_orig2;
336: PetscObjectGetComm((PetscObject)Gmat_2,&comm);
337: MPI_Comm_rank(comm, &rank);
338: MPI_Comm_size(comm, &size);
339: MatGetOwnershipRange(Gmat_1,&my0,&Iend);
341: if (PETSC_FALSE) {
342: PetscViewer viewer; char fname[32]; static int llev=0;
343: sprintf(fname,"Gmat2_%d.m",llev++);
344: PetscViewerASCIIOpen(comm,fname,&viewer);
345: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
346: MatView(Gmat_2, viewer);
347: PetscViewerDestroy(&viewer);
348: }
350: /* get submatrices */
351: PetscObjectTypeCompare((PetscObject)Gmat_1, MATMPIAIJ, &isMPI);
352: if (isMPI) {
353: /* grab matrix objects */
354: mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data;
355: mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data;
356: matA_1 = (Mat_SeqAIJ*)mpimat_1->A->data;
357: matB_1 = (Mat_SeqAIJ*)mpimat_1->B->data;
358: matA_2 = (Mat_SeqAIJ*)mpimat_2->A->data;
359: matB_2 = (Mat_SeqAIJ*)mpimat_2->B->data;
361: /* force compressed row storage for B matrix in AuxMat */
362: MatCheckCompressedRow(mpimat_1->B,matB_1->nonzerorowcnt,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0);
364: PetscMalloc1(nloc, &lid_cprowID_1);
365: for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
366: for (ix=0; ix<matB_1->compressedrow.nrows; ix++) {
367: PetscInt lid = matB_1->compressedrow.rindex[ix];
368: lid_cprowID_1[lid] = ix;
369: }
370: } else {
371: matA_1 = (Mat_SeqAIJ*)Gmat_1->data;
372: matA_2 = (Mat_SeqAIJ*)Gmat_2->data;
373: lid_cprowID_1 = NULL;
374: }
375: if (nloc>0) {
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)");
380: }
381: /* get state of locals and selected gid for deleted */
382: PetscMalloc1(nloc, &lid_state);
383: PetscMalloc1(nloc, &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: PetscMalloc1((nvec+!!has_const)*mlocal,&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: PetscMalloc1(out_data_stride*nSAvec, &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: PetscMalloc1((Mdata*N), &qqc);
794: PetscMalloc1((M*N), &qqr);
795: PetscMalloc1(N, &TAU);
796: PetscMalloc1(LWORK, &WORK);
797: PetscMalloc1(M, &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: PetscMalloc1(nloc, &permute);
1004: PetscMalloc1(nloc, &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: PetscMalloc1(nloc, &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: PetscMalloc1(stride*bs*data_cols, &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: PetscMalloc1(nloc, &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: PetscMalloc1(stride, &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: PetscMalloc1(nloc, &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 */
1214: #if defined PETSC_USE_LOG
1215: PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1216: #endif
1217: return(0);
1218: }
1220: /* -------------------------------------------------------------------------- */
1221: /*
1222: PCGAMGOptprol_AGG
1224: Input Parameter:
1225: . pc - this
1226: . Amat - matrix on this fine level
1227: In/Output Parameter:
1228: . a_P_out - prolongation operator to the next level
1229: */
1232: PetscErrorCode PCGAMGOptprol_AGG(PC pc,const Mat Amat,Mat *a_P)
1233: {
1235: PC_MG *mg = (PC_MG*)pc->data;
1236: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1237: const PetscInt verbose = pc_gamg->verbose;
1238: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1239: PetscInt jj;
1240: PetscMPIInt rank,size;
1241: Mat Prol = *a_P;
1242: MPI_Comm comm;
1245: PetscObjectGetComm((PetscObject)Amat,&comm);
1246: #if defined PETSC_USE_LOG
1247: PetscLogEventBegin(PC_GAMGOptprol_AGG,0,0,0,0);
1248: #endif
1249: MPI_Comm_rank(comm, &rank);
1250: MPI_Comm_size(comm, &size);
1252: /* smooth P0 */
1253: for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1254: Mat tMat;
1255: Vec diag;
1256: PetscReal alpha, emax, emin;
1257: #if defined PETSC_GAMG_USE_LOG
1258: PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);
1259: #endif
1260: if (jj == 0) {
1261: KSP eksp;
1262: Vec bb, xx;
1263: PC epc;
1264: MatGetVecs(Amat, &bb, 0);
1265: MatGetVecs(Amat, &xx, 0);
1266: {
1267: PetscRandom rctx;
1268: PetscRandomCreate(comm,&rctx);
1269: PetscRandomSetFromOptions(rctx);
1270: VecSetRandom(bb,rctx);
1271: PetscRandomDestroy(&rctx);
1272: }
1274: /* zeroing out BC rows -- needed for crazy matrices */
1275: {
1276: PetscInt Istart,Iend,ncols,jj,Ii;
1277: PetscScalar zero = 0.0;
1278: MatGetOwnershipRange(Amat, &Istart, &Iend);
1279: for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
1280: MatGetRow(Amat,Ii,&ncols,0,0);
1281: if (ncols <= 1) {
1282: VecSetValues(bb, 1, &Ii, &zero, INSERT_VALUES);
1283: }
1284: MatRestoreRow(Amat,Ii,&ncols,0,0);
1285: }
1286: VecAssemblyBegin(bb);
1287: VecAssemblyEnd(bb);
1288: }
1290: KSPCreate(comm,&eksp);
1291: KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
1292: KSPSetNormType(eksp, KSP_NORM_NONE);
1293: KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);
1294: KSPAppendOptionsPrefix(eksp, "gamg_est_");
1295: KSPSetFromOptions(eksp);
1297: KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
1298: KSPSetOperators(eksp, Amat, Amat);
1299: KSPSetComputeSingularValues(eksp,PETSC_TRUE);
1301: KSPGetPC(eksp, &epc);
1302: PCSetType(epc, PCJACOBI); /* smoother in smoothed agg. */
1304: /* solve - keep stuff out of logging */
1305: PetscLogEventDeactivate(KSP_Solve);
1306: PetscLogEventDeactivate(PC_Apply);
1307: KSPSolve(eksp, bb, xx);
1308: PetscLogEventActivate(KSP_Solve);
1309: PetscLogEventActivate(PC_Apply);
1311: KSPComputeExtremeSingularValues(eksp, &emax, &emin);
1312: if (verbose > 0) {
1313: PetscPrintf(comm,"\t\t\t%s smooth P0: max eigen=%e min=%e PC=%s\n",
1314: __FUNCT__,emax,emin,PCJACOBI);
1315: }
1316: VecDestroy(&xx);
1317: VecDestroy(&bb);
1318: KSPDestroy(&eksp);
1320: if (pc_gamg->emax_id == -1) {
1321: PetscObjectComposedDataRegister(&pc_gamg->emax_id);
1322: if (pc_gamg->emax_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->emax_id == -1");
1323: }
1324: PetscObjectComposedDataSetScalar((PetscObject)Amat, pc_gamg->emax_id, emax);
1325: }
1327: /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1328: MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);
1329: MatGetVecs(Amat, &diag, 0);
1330: MatGetDiagonal(Amat, diag); /* effectively PCJACOBI */
1331: VecReciprocal(diag);
1332: MatDiagonalScale(tMat, diag, 0);
1333: VecDestroy(&diag);
1334: alpha = -1.4/emax;
1335: MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);
1336: MatDestroy(&Prol);
1337: Prol = tMat;
1338: #if defined PETSC_GAMG_USE_LOG
1339: PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);
1340: #endif
1341: }
1342: #if defined PETSC_USE_LOG
1343: PetscLogEventEnd(PC_GAMGOptprol_AGG,0,0,0,0);
1344: #endif
1345: *a_P = Prol;
1346: return(0);
1347: }
1349: /* -------------------------------------------------------------------------- */
1350: /*
1351: PCCreateGAMG_AGG
1353: Input Parameter:
1354: . pc -
1355: */
1358: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1359: {
1361: PC_MG *mg = (PC_MG*)pc->data;
1362: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1363: PC_GAMG_AGG *pc_gamg_agg;
1366: /* create sub context for SA */
1367: PetscNewLog(pc,&pc_gamg_agg);
1368: pc_gamg->subctx = pc_gamg_agg;
1370: pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1371: pc_gamg->ops->destroy = PCDestroy_GAMG_AGG;
1372: /* reset does not do anything; setup not virtual */
1374: /* set internal function pointers */
1375: pc_gamg->ops->graph = PCGAMGgraph_AGG;
1376: pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG;
1377: pc_gamg->ops->prolongator = PCGAMGProlongator_AGG;
1378: pc_gamg->ops->optprol = PCGAMGOptprol_AGG;
1380: pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1382: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);
1383: return(0);
1384: }