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