Actual source code: geo.c
petsc-3.10.5 2019-03-28
1: /*
2: GAMG geometric-algebric multiogrid PC - Mark Adams 2011
3: */
5: #include <../src/ksp/pc/impls/gamg/gamg.h>
7: #if defined(PETSC_HAVE_TRIANGLE)
8: #define REAL PetscReal
9: #include <triangle.h>
10: #endif
12: #include <petscblaslapack.h>
14: /* Private context for the GAMG preconditioner */
15: typedef struct {
16: PetscInt lid; /* local vertex index */
17: PetscInt degree; /* vertex degree */
18: } GAMGNode;
20: PETSC_STATIC_INLINE int petsc_geo_mg_compare(const void *a, const void *b)
21: {
22: return (((GAMGNode*)a)->degree - ((GAMGNode*)b)->degree);
23: }
25: /* -------------------------------------------------------------------------- */
26: /*
27: PCSetCoordinates_GEO
29: Input Parameter:
30: . pc - the preconditioner context
31: */
32: PetscErrorCode PCSetCoordinates_GEO(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
33: {
34: PC_MG *mg = (PC_MG*)pc->data;
35: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
37: PetscInt arrsz,bs,my0,kk,ii,nloc,Iend;
38: Mat Amat = pc->pmat;
42: MatGetBlockSize(Amat, &bs);
43: MatGetOwnershipRange(Amat, &my0, &Iend);
44: nloc = (Iend-my0)/bs;
46: if (nloc!=a_nloc) SETERRQ2(PetscObjectComm((PetscObject)Amat),PETSC_ERR_ARG_WRONG, "Stokes not supported nloc = %D %D.",a_nloc,nloc);
47: if ((Iend-my0)%bs!=0) SETERRQ1(PetscObjectComm((PetscObject)Amat),PETSC_ERR_ARG_WRONG, "Bad local size %D.",nloc);
49: pc_gamg->data_cell_rows = 1;
50: if (!coords && nloc > 0) SETERRQ(PetscObjectComm((PetscObject)Amat),PETSC_ERR_ARG_WRONG, "Need coordinates for pc_gamg_type 'geo'.");
51: pc_gamg->data_cell_cols = ndm; /* coordinates */
53: arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;
55: /* create data - syntactic sugar that should be refactored at some point */
56: if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
57: PetscFree(pc_gamg->data);
58: PetscMalloc1(arrsz+1, &pc_gamg->data);
59: }
60: for (kk=0; kk<arrsz; kk++) pc_gamg->data[kk] = -999.;
61: pc_gamg->data[arrsz] = -99.;
62: /* copy data in - column oriented */
63: for (kk = 0; kk < nloc; kk++) {
64: for (ii = 0; ii < ndm; ii++) {
65: pc_gamg->data[ii*nloc + kk] = coords[kk*ndm + ii];
66: }
67: }
68: if (pc_gamg->data[arrsz] != -99.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->data[arrsz %D] %g != -99.",arrsz,pc_gamg->data[arrsz]);
69: pc_gamg->data_sz = arrsz;
70: return(0);
71: }
73: /* -------------------------------------------------------------------------- */
74: /*
75: PCSetData_GEO
77: Input Parameter:
78: . pc -
79: */
80: PetscErrorCode PCSetData_GEO(PC pc, Mat m)
81: {
83: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"GEO MG needs coordinates");
84: }
86: /* -------------------------------------------------------------------------- */
87: /*
88: PCSetFromOptions_GEO
90: Input Parameter:
91: . pc -
92: */
93: PetscErrorCode PCSetFromOptions_GEO(PetscOptionItems *PetscOptionsObject,PC pc)
94: {
98: PetscOptionsHead(PetscOptionsObject,"GAMG-GEO options");
99: {
100: /* -pc_gamg_sa_nsmooths */
101: /* pc_gamg_sa->smooths = 0; */
102: /* PetscOptionsInt("-pc_gamg_agg_nsmooths", */
103: /* "smoothing steps for smoothed aggregation, usually 1 (0)", */
104: /* "PCGAMGSetNSmooths_AGG", */
105: /* pc_gamg_sa->smooths, */
106: /* &pc_gamg_sa->smooths, */
107: /* &flag); */
108: /* */
109: }
110: PetscOptionsTail();
111: return(0);
112: }
114: /* -------------------------------------------------------------------------- */
115: /*
116: triangulateAndFormProl
118: Input Parameter:
119: . selected_2 - list of selected local ID, includes selected ghosts
120: . data_stride -
121: . coords[2*data_stride] - column vector of local coordinates w/ ghosts
122: . nselected_1 - selected IDs that go with base (1) graph includes selected ghosts
123: . clid_lid_1[nselected_1] - lids of selected (c) nodes ???????????
124: . agg_lists_1 - list of aggregates selected_1 vertices of aggregate unselected vertices
125: . crsGID[selected.size()] - global index for prolongation operator
126: . bs - block size
127: Output Parameter:
128: . a_Prol - prolongation operator
129: . a_worst_best - measure of worst missed fine vertex, 0 is no misses
130: */
131: static PetscErrorCode triangulateAndFormProl(IS selected_2,PetscInt data_stride,PetscReal coords[],PetscInt nselected_1,const PetscInt clid_lid_1[],const PetscCoarsenData *agg_lists_1,
132: const PetscInt crsGID[],PetscInt bs,Mat a_Prol,PetscReal *a_worst_best)
133: {
134: #if defined(PETSC_HAVE_TRIANGLE)
135: PetscErrorCode ierr;
136: PetscInt jj,tid,tt,idx,nselected_2;
137: struct triangulateio in,mid;
138: const PetscInt *selected_idx_2;
139: PetscMPIInt rank;
140: PetscInt Istart,Iend,nFineLoc,myFine0;
141: int kk,nPlotPts,sid;
142: MPI_Comm comm;
143: PetscReal tm;
146: PetscObjectGetComm((PetscObject)a_Prol,&comm);
147: MPI_Comm_rank(comm,&rank);
148: ISGetSize(selected_2, &nselected_2);
149: if (nselected_2 == 1 || nselected_2 == 2) { /* 0 happens on idle processors */
150: *a_worst_best = 100.0; /* this will cause a stop, but not globalized (should not happen) */
151: } else *a_worst_best = 0.0;
152: MPIU_Allreduce(a_worst_best, &tm, 1, MPIU_REAL, MPIU_MAX, comm);
153: if (tm > 0.0) {
154: *a_worst_best = 100.0;
155: return(0);
156: }
157: MatGetOwnershipRange(a_Prol, &Istart, &Iend);
158: nFineLoc = (Iend-Istart)/bs; myFine0 = Istart/bs;
159: nPlotPts = nFineLoc; /* locals */
160: /* traingle */
161: /* Define input points - in*/
162: in.numberofpoints = nselected_2;
163: in.numberofpointattributes = 0;
164: /* get nselected points */
165: PetscMalloc1(2*nselected_2, &in.pointlist);
166: ISGetIndices(selected_2, &selected_idx_2);
168: for (kk=0,sid=0; kk<nselected_2; kk++,sid += 2) {
169: PetscInt lid = selected_idx_2[kk];
170: in.pointlist[sid] = coords[lid];
171: in.pointlist[sid+1] = coords[data_stride + lid];
172: if (lid>=nFineLoc) nPlotPts++;
173: }
174: if (sid != 2*nselected_2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %D != 2*nselected_2 %D",sid,nselected_2);
176: in.numberofsegments = 0;
177: in.numberofedges = 0;
178: in.numberofholes = 0;
179: in.numberofregions = 0;
180: in.trianglelist = 0;
181: in.segmentmarkerlist = 0;
182: in.pointattributelist = 0;
183: in.pointmarkerlist = 0;
184: in.triangleattributelist = 0;
185: in.trianglearealist = 0;
186: in.segmentlist = 0;
187: in.holelist = 0;
188: in.regionlist = 0;
189: in.edgelist = 0;
190: in.edgemarkerlist = 0;
191: in.normlist = 0;
193: /* triangulate */
194: mid.pointlist = 0; /* Not needed if -N switch used. */
195: /* Not needed if -N switch used or number of point attributes is zero: */
196: mid.pointattributelist = 0;
197: mid.pointmarkerlist = 0; /* Not needed if -N or -B switch used. */
198: mid.trianglelist = 0; /* Not needed if -E switch used. */
199: /* Not needed if -E switch used or number of triangle attributes is zero: */
200: mid.triangleattributelist = 0;
201: mid.neighborlist = 0; /* Needed only if -n switch used. */
202: /* Needed only if segments are output (-p or -c) and -P not used: */
203: mid.segmentlist = 0;
204: /* Needed only if segments are output (-p or -c) and -P and -B not used: */
205: mid.segmentmarkerlist = 0;
206: mid.edgelist = 0; /* Needed only if -e switch used. */
207: mid.edgemarkerlist = 0; /* Needed if -e used and -B not used. */
208: mid.numberoftriangles = 0;
210: /* Triangulate the points. Switches are chosen to read and write a */
211: /* PSLG (p), preserve the convex hull (c), number everything from */
212: /* zero (z), assign a regional attribute to each element (A), and */
213: /* produce an edge list (e), a Voronoi diagram (v), and a triangle */
214: /* neighbor list (n). */
215: if (nselected_2 != 0) { /* inactive processor */
216: char args[] = "npczQ"; /* c is needed ? */
217: triangulate(args, &in, &mid, (struct triangulateio*) NULL);
218: /* output .poly files for 'showme' */
219: if (!PETSC_TRUE) {
220: static int level = 1;
221: FILE *file; char fname[32];
223: sprintf(fname,"C%d_%d.poly",level,rank); file = fopen(fname, "w");
224: /*First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>*/
225: fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0);
226: /*Following lines: <vertex #> <x> <y> */
227: for (kk=0,sid=0; kk<in.numberofpoints; kk++,sid += 2) {
228: fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]);
229: }
230: /*One line: <# of segments> <# of boundary markers (0 or 1)> */
231: fprintf(file, "%d %d\n",0,0);
232: /*Following lines: <segment #> <endpoint> <endpoint> [boundary marker] */
233: /* One line: <# of holes> */
234: fprintf(file, "%d\n",0);
235: /* Following lines: <hole #> <x> <y> */
236: /* Optional line: <# of regional attributes and/or area constraints> */
237: /* Optional following lines: <region #> <x> <y> <attribute> <maximum area> */
238: fclose(file);
240: /* elems */
241: sprintf(fname,"C%d_%d.ele",level,rank); file = fopen(fname, "w");
242: /* First line: <# of triangles> <nodes per triangle> <# of attributes> */
243: fprintf(file, "%d %d %d\n",mid.numberoftriangles,3,0);
244: /* Remaining lines: <triangle #> <node> <node> <node> ... [attributes] */
245: for (kk=0,sid=0; kk<mid.numberoftriangles; kk++,sid += 3) {
246: fprintf(file, "%d %d %d %d\n",kk,mid.trianglelist[sid],mid.trianglelist[sid+1],mid.trianglelist[sid+2]);
247: }
248: fclose(file);
250: sprintf(fname,"C%d_%d.node",level,rank); file = fopen(fname, "w");
251: /* First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> */
252: /* fprintf(file, "%d %d %d %d\n",in.numberofpoints,2,0,0); */
253: fprintf(file, "%d %d %d %d\n",nPlotPts,2,0,0);
254: /*Following lines: <vertex #> <x> <y> */
255: for (kk=0,sid=0; kk<in.numberofpoints; kk++,sid+=2) {
256: fprintf(file, "%d %e %e\n",kk,in.pointlist[sid],in.pointlist[sid+1]);
257: }
259: sid /= 2;
260: for (jj=0; jj<nFineLoc; jj++) {
261: PetscBool sel = PETSC_TRUE;
262: for (kk=0; kk<nselected_2 && sel; kk++) {
263: PetscInt lid = selected_idx_2[kk];
264: if (lid == jj) sel = PETSC_FALSE;
265: }
266: if (sel) fprintf(file, "%d %e %e\n",sid++,coords[jj],coords[data_stride + jj]);
267: }
268: fclose(file);
269: if (sid != nPlotPts) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"sid %D != nPlotPts %D",sid,nPlotPts);
270: level++;
271: }
272: }
273: #if defined PETSC_GAMG_USE_LOG
274: PetscLogEventBegin(petsc_gamg_setup_events[FIND_V],0,0,0,0);
275: #endif
276: { /* form P - setup some maps */
277: PetscInt clid,mm,*nTri,*node_tri;
279: PetscMalloc2(nselected_2, &node_tri,nselected_2, &nTri);
281: /* need list of triangles on node */
282: for (kk=0; kk<nselected_2; kk++) nTri[kk] = 0;
283: for (tid=0,kk=0; tid<mid.numberoftriangles; tid++) {
284: for (jj=0; jj<3; jj++) {
285: PetscInt cid = mid.trianglelist[kk++];
286: if (nTri[cid] == 0) node_tri[cid] = tid;
287: nTri[cid]++;
288: }
289: }
290: #define EPS 1.e-12
291: /* find points and set prolongation */
292: for (mm = clid = 0; mm < nFineLoc; mm++) {
293: PetscBool ise;
294: PetscCDEmptyAt(agg_lists_1,mm,&ise);
295: if (!ise) {
296: const PetscInt lid = mm;
297: PetscScalar AA[3][3];
298: PetscBLASInt N=3,NRHS=1,LDA=3,IPIV[3],LDB=3,INFO;
299: PetscCDIntNd *pos;
301: PetscCDGetHeadPos(agg_lists_1,lid,&pos);
302: while (pos) {
303: PetscInt flid;
304: PetscCDIntNdGetID(pos, &flid);
305: PetscCDGetNextPos(agg_lists_1,lid,&pos);
307: if (flid < nFineLoc) { /* could be a ghost */
308: PetscInt bestTID = -1; PetscReal best_alpha = 1.e10;
309: const PetscInt fgid = flid + myFine0;
310: /* compute shape function for gid */
311: const PetscReal fcoord[3] = {coords[flid],coords[data_stride+flid],1.0};
312: PetscBool haveit =PETSC_FALSE; PetscScalar alpha[3]; PetscInt clids[3];
314: /* look for it */
315: for (tid = node_tri[clid], jj=0;
316: jj < 5 && !haveit && tid != -1;
317: jj++) {
318: for (tt=0; tt<3; tt++) {
319: PetscInt cid2 = mid.trianglelist[3*tid + tt];
320: PetscInt lid2 = selected_idx_2[cid2];
321: AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
322: clids[tt] = cid2; /* store for interp */
323: }
325: for (tt=0; tt<3; tt++) alpha[tt] = (PetscScalar)fcoord[tt];
327: /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
328: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
329: {
330: PetscBool have=PETSC_TRUE; PetscReal lowest=1.e10;
331: for (tt = 0, idx = 0; tt < 3; tt++) {
332: if (PetscRealPart(alpha[tt]) > (1.0+EPS) || PetscRealPart(alpha[tt]) < -EPS) have = PETSC_FALSE;
333: if (PetscRealPart(alpha[tt]) < lowest) {
334: lowest = PetscRealPart(alpha[tt]);
335: idx = tt;
336: }
337: }
338: haveit = have;
339: }
340: tid = mid.neighborlist[3*tid + idx];
341: }
343: if (!haveit) {
344: /* brute force */
345: for (tid=0; tid<mid.numberoftriangles && !haveit; tid++) {
346: for (tt=0; tt<3; tt++) {
347: PetscInt cid2 = mid.trianglelist[3*tid + tt];
348: PetscInt lid2 = selected_idx_2[cid2];
349: AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
350: clids[tt] = cid2; /* store for interp */
351: }
352: for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt];
353: /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
354: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
355: {
356: PetscBool have=PETSC_TRUE; PetscReal worst=0.0, v;
357: for (tt=0; tt<3 && have; tt++) {
358: if (PetscRealPart(alpha[tt]) > 1.0+EPS || PetscRealPart(alpha[tt]) < -EPS) have=PETSC_FALSE;
359: if ((v=PetscAbs(PetscRealPart(alpha[tt])-0.5)) > worst) worst = v;
360: }
361: if (worst < best_alpha) {
362: best_alpha = worst; bestTID = tid;
363: }
364: haveit = have;
365: }
366: }
367: }
368: if (!haveit) {
369: if (best_alpha > *a_worst_best) *a_worst_best = best_alpha;
370: /* use best one */
371: for (tt=0; tt<3; tt++) {
372: PetscInt cid2 = mid.trianglelist[3*bestTID + tt];
373: PetscInt lid2 = selected_idx_2[cid2];
374: AA[tt][0] = coords[lid2]; AA[tt][1] = coords[data_stride + lid2]; AA[tt][2] = 1.0;
375: clids[tt] = cid2; /* store for interp */
376: }
377: for (tt=0; tt<3; tt++) alpha[tt] = fcoord[tt];
378: /* SUBROUTINE DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO) */
379: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&N, &NRHS, (PetscScalar*)AA, &LDA, IPIV, alpha, &LDB, &INFO));
380: }
382: /* put in row of P */
383: for (idx=0; idx<3; idx++) {
384: PetscScalar shp = alpha[idx];
385: if (PetscAbs(PetscRealPart(shp)) > 1.e-6) {
386: PetscInt cgid = crsGID[clids[idx]];
387: PetscInt jj = cgid*bs, ii = fgid*bs; /* need to gloalize */
388: for (tt=0; tt < bs; tt++, ii++, jj++) {
389: MatSetValues(a_Prol,1,&ii,1,&jj,&shp,INSERT_VALUES);
390: }
391: }
392: }
393: }
394: } /* aggregates iterations */
395: clid++;
396: } /* a coarse agg */
397: } /* for all fine nodes */
399: ISRestoreIndices(selected_2, &selected_idx_2);
400: MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
401: MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);
403: PetscFree2(node_tri,nTri);
404: }
405: #if defined PETSC_GAMG_USE_LOG
406: PetscLogEventEnd(petsc_gamg_setup_events[FIND_V],0,0,0,0);
407: #endif
408: free(mid.trianglelist);
409: free(mid.neighborlist);
410: PetscFree(in.pointlist);
411: return(0);
412: #else
413: SETERRQ(PetscObjectComm((PetscObject)a_Prol),PETSC_ERR_PLIB,"configure with TRIANGLE to use geometric MG");
414: #endif
415: }
416: /* -------------------------------------------------------------------------- */
417: /*
418: getGIDsOnSquareGraph - square graph, get
420: Input Parameter:
421: . nselected_1 - selected local indices (includes ghosts in input Gmat1)
422: . clid_lid_1 - [nselected_1] lids of selected nodes
423: . Gmat1 - graph that goes with 'selected_1'
424: Output Parameter:
425: . a_selected_2 - selected local indices (includes ghosts in output a_Gmat_2)
426: . a_Gmat_2 - graph that is squared of 'Gmat_1'
427: . a_crsGID[a_selected_2.size()] - map of global IDs of coarse grid nodes
428: */
429: static PetscErrorCode getGIDsOnSquareGraph(PetscInt nselected_1,const PetscInt clid_lid_1[],const Mat Gmat1,IS *a_selected_2,Mat *a_Gmat_2,PetscInt **a_crsGID)
430: {
432: PetscMPIInt size;
433: PetscInt *crsGID, kk,my0,Iend,nloc;
434: MPI_Comm comm;
437: PetscObjectGetComm((PetscObject)Gmat1,&comm);
438: MPI_Comm_size(comm,&size);
439: MatGetOwnershipRange(Gmat1,&my0,&Iend); /* AIJ */
440: nloc = Iend - my0; /* this does not change */
442: if (size == 1) { /* not much to do in serial */
443: PetscMalloc1(nselected_1, &crsGID);
444: for (kk=0; kk<nselected_1; kk++) crsGID[kk] = kk;
445: *a_Gmat_2 = 0;
446: ISCreateGeneral(PETSC_COMM_SELF,nselected_1,clid_lid_1,PETSC_COPY_VALUES,a_selected_2);
447: } else {
448: PetscInt idx,num_fine_ghosts,num_crs_ghost,myCrs0;
449: Mat_MPIAIJ *mpimat2;
450: Mat Gmat2;
451: Vec locState;
452: PetscScalar *cpcol_state;
454: /* scan my coarse zero gid, set 'lid_state' with coarse GID */
455: kk = nselected_1;
456: MPI_Scan(&kk, &myCrs0, 1, MPIU_INT, MPI_SUM, comm);
457: myCrs0 -= nselected_1;
459: if (a_Gmat_2) { /* output */
460: /* grow graph to get wider set of selected vertices to cover fine grid, invalidates 'llist' */
461: MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);
462: *a_Gmat_2 = Gmat2; /* output */
463: } else Gmat2 = Gmat1; /* use local to get crsGIDs at least */
464: /* get coarse grid GIDS for selected (locals and ghosts) */
465: mpimat2 = (Mat_MPIAIJ*)Gmat2->data;
466: MatCreateVecs(Gmat2, &locState, 0);
467: VecSet(locState, (PetscScalar)(PetscReal)(-1)); /* set with UNKNOWN state */
468: for (kk=0; kk<nselected_1; kk++) {
469: PetscInt fgid = clid_lid_1[kk] + my0;
470: PetscScalar v = (PetscScalar)(kk+myCrs0);
471: VecSetValues(locState, 1, &fgid, &v, INSERT_VALUES); /* set with PID */
472: }
473: VecAssemblyBegin(locState);
474: VecAssemblyEnd(locState);
475: VecScatterBegin(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);
476: VecScatterEnd(mpimat2->Mvctx,locState,mpimat2->lvec,INSERT_VALUES,SCATTER_FORWARD);
477: VecGetLocalSize(mpimat2->lvec, &num_fine_ghosts);
478: VecGetArray(mpimat2->lvec, &cpcol_state);
479: for (kk=0,num_crs_ghost=0; kk<num_fine_ghosts; kk++) {
480: if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) num_crs_ghost++;
481: }
482: PetscMalloc1(nselected_1+num_crs_ghost, &crsGID); /* output */
483: {
484: PetscInt *selected_set;
485: PetscMalloc1(nselected_1+num_crs_ghost, &selected_set);
486: /* do ghost of 'crsGID' */
487: for (kk=0,idx=nselected_1; kk<num_fine_ghosts; kk++) {
488: if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
489: PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
490: selected_set[idx] = nloc + kk;
491: crsGID[idx++] = cgid;
492: }
493: }
494: if (idx != (nselected_1+num_crs_ghost)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"idx %D != (nselected_1 %D + num_crs_ghost %D)",idx,nselected_1,num_crs_ghost);
495: VecRestoreArray(mpimat2->lvec, &cpcol_state);
496: /* do locals in 'crsGID' */
497: VecGetArray(locState, &cpcol_state);
498: for (kk=0,idx=0; kk<nloc; kk++) {
499: if ((PetscInt)PetscRealPart(cpcol_state[kk]) != -1) {
500: PetscInt cgid = (PetscInt)PetscRealPart(cpcol_state[kk]);
501: selected_set[idx] = kk;
502: crsGID[idx++] = cgid;
503: }
504: }
505: if (idx != nselected_1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"idx %D != nselected_1 %D",idx,nselected_1);
506: VecRestoreArray(locState, &cpcol_state);
508: if (a_selected_2 != 0) { /* output */
509: ISCreateGeneral(PETSC_COMM_SELF,(nselected_1+num_crs_ghost),selected_set,PETSC_OWN_POINTER,a_selected_2);
510: } else {
511: PetscFree(selected_set);
512: }
513: }
514: VecDestroy(&locState);
515: }
516: *a_crsGID = crsGID; /* output */
517: return(0);
518: }
520: /* -------------------------------------------------------------------------- */
521: /*
522: PCGAMGGraph_GEO
524: Input Parameter:
525: . pc - this
526: . Amat - matrix on this fine level
527: Output Parameter:
528: . a_Gmat
529: */
530: PetscErrorCode PCGAMGGraph_GEO(PC pc,Mat Amat,Mat *a_Gmat)
531: {
532: PetscErrorCode ierr;
533: PC_MG *mg = (PC_MG*)pc->data;
534: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
535: const PetscReal vfilter = pc_gamg->threshold[0];
536: MPI_Comm comm;
537: Mat Gmat;
538: PetscBool set,flg,symm;
541: PetscObjectGetComm((PetscObject)Amat,&comm);
542: PetscLogEventBegin(PC_GAMGGraph_GEO,0,0,0,0);
544: MatIsSymmetricKnown(Amat, &set, &flg);
545: symm = (PetscBool)!(set && flg);
547: PCGAMGCreateGraph(Amat, &Gmat);
548: PCGAMGFilterGraph(&Gmat, vfilter, symm);
550: *a_Gmat = Gmat;
551: PetscLogEventEnd(PC_GAMGGraph_GEO,0,0,0,0);
552: return(0);
553: }
555: /* -------------------------------------------------------------------------- */
556: /*
557: PCGAMGCoarsen_GEO
559: Input Parameter:
560: . a_pc - this
561: . a_Gmat - graph
562: Output Parameter:
563: . a_llist_parent - linked list from selected indices for data locality only
564: */
565: PetscErrorCode PCGAMGCoarsen_GEO(PC a_pc,Mat *a_Gmat,PetscCoarsenData **a_llist_parent)
566: {
568: PetscInt Istart,Iend,nloc,kk,Ii,ncols;
569: IS perm;
570: GAMGNode *gnodes;
571: PetscInt *permute;
572: Mat Gmat = *a_Gmat;
573: MPI_Comm comm;
574: MatCoarsen crs;
577: PetscObjectGetComm((PetscObject)a_pc,&comm);
578: PetscLogEventBegin(PC_GAMGCoarsen_GEO,0,0,0,0);
579: MatGetOwnershipRange(Gmat, &Istart, &Iend);
580: nloc = (Iend-Istart);
582: /* create random permutation with sort for geo-mg */
583: PetscMalloc1(nloc, &gnodes);
584: PetscMalloc1(nloc, &permute);
586: for (Ii=Istart; Ii<Iend; Ii++) { /* locals only? */
587: MatGetRow(Gmat,Ii,&ncols,0,0);
588: {
589: PetscInt lid = Ii - Istart;
590: gnodes[lid].lid = lid;
591: gnodes[lid].degree = ncols;
592: }
593: MatRestoreRow(Gmat,Ii,&ncols,0,0);
594: }
595: if (PETSC_TRUE) {
596: PetscRandom rand;
597: PetscBool *bIndexSet;
598: PetscReal rr;
599: PetscInt iSwapIndex;
601: PetscRandomCreate(comm,&rand);
602: PetscCalloc1(nloc, &bIndexSet);
603: for (Ii = 0; Ii < nloc; Ii++) {
604: PetscRandomGetValueReal(rand,&rr);
605: iSwapIndex = (PetscInt) (rr*nloc);
606: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
607: GAMGNode iTemp = gnodes[iSwapIndex];
608: gnodes[iSwapIndex] = gnodes[Ii];
609: gnodes[Ii] = iTemp;
610: bIndexSet[Ii] = PETSC_TRUE;
611: bIndexSet[iSwapIndex] = PETSC_TRUE;
612: }
613: }
614: PetscRandomDestroy(&rand);
615: PetscFree(bIndexSet);
616: }
617: /* only sort locals */
618: qsort(gnodes, nloc, sizeof(GAMGNode), petsc_geo_mg_compare);
619: /* create IS of permutation */
620: for (kk=0; kk<nloc; kk++) permute[kk] = gnodes[kk].lid; /* locals only */
621: PetscFree(gnodes);
622: ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_OWN_POINTER, &perm);
624: /* get MIS aggs */
626: MatCoarsenCreate(comm, &crs);
627: MatCoarsenSetType(crs, MATCOARSENMIS);
628: MatCoarsenSetGreedyOrdering(crs, perm);
629: MatCoarsenSetAdjacency(crs, Gmat);
630: MatCoarsenSetStrictAggs(crs, PETSC_FALSE);
631: MatCoarsenApply(crs);
632: MatCoarsenGetData(crs, a_llist_parent);
633: MatCoarsenDestroy(&crs);
635: ISDestroy(&perm);
636: PetscLogEventEnd(PC_GAMGCoarsen_GEO,0,0,0,0);
637: return(0);
638: }
640: /* -------------------------------------------------------------------------- */
641: /*
642: PCGAMGProlongator_GEO
644: Input Parameter:
645: . pc - this
646: . Amat - matrix on this fine level
647: . Graph - used to get ghost data for nodes in
648: . selected_1 - [nselected]
649: . agg_lists - [nselected]
650: Output Parameter:
651: . a_P_out - prolongation operator to the next level
652: */
653: PetscErrorCode PCGAMGProlongator_GEO(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
654: {
655: PC_MG *mg = (PC_MG*)pc->data;
656: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
657: const PetscInt dim = pc_gamg->data_cell_cols, data_cols = pc_gamg->data_cell_cols;
659: PetscInt Istart,Iend,nloc,my0,jj,kk,ncols,nLocalSelected,bs,*clid_flid;
660: Mat Prol;
661: PetscMPIInt rank, size;
662: MPI_Comm comm;
663: IS selected_2,selected_1;
664: const PetscInt *selected_idx;
665: MatType mtype;
668: PetscObjectGetComm((PetscObject)Amat,&comm);
669: PetscLogEventBegin(PC_GAMGProlongator_GEO,0,0,0,0);
670: MPI_Comm_rank(comm,&rank);
671: MPI_Comm_size(comm,&size);
672: MatGetOwnershipRange(Amat, &Istart, &Iend);
673: MatGetBlockSize(Amat, &bs);
674: nloc = (Iend-Istart)/bs; my0 = Istart/bs;
675: if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) % bs %D",Iend,Istart,bs);
677: /* get 'nLocalSelected' */
678: PetscCDGetMIS(agg_lists, &selected_1);
679: ISGetSize(selected_1, &jj);
680: PetscMalloc1(jj, &clid_flid);
681: ISGetIndices(selected_1, &selected_idx);
682: for (kk=0,nLocalSelected=0; kk<jj; kk++) {
683: PetscInt lid = selected_idx[kk];
684: if (lid<nloc) {
685: MatGetRow(Gmat,lid+my0,&ncols,0,0);
686: if (ncols>1) clid_flid[nLocalSelected++] = lid; /* fiter out singletons */
687: MatRestoreRow(Gmat,lid+my0,&ncols,0,0);
688: }
689: }
690: ISRestoreIndices(selected_1, &selected_idx);
691: ISDestroy(&selected_1); /* this is selected_1 in serial */
693: /* create prolongator matrix */
694: MatGetType(Amat,&mtype);
695: MatCreate(comm, &Prol);
696: MatSetSizes(Prol,nloc*bs,nLocalSelected*bs,PETSC_DETERMINE,PETSC_DETERMINE);
697: MatSetBlockSizes(Prol, bs, bs);
698: MatSetType(Prol, mtype);
699: MatSeqAIJSetPreallocation(Prol,3*data_cols,NULL);
700: MatMPIAIJSetPreallocation(Prol,3*data_cols,NULL,3*data_cols,NULL);
702: /* can get all points "removed" - but not on geomg */
703: MatGetSize(Prol, &kk, &jj);
704: if (!jj) {
705: PetscInfo(pc,"ERROE: no selected points on coarse grid\n");
706: PetscFree(clid_flid);
707: MatDestroy(&Prol);
708: *a_P_out = NULL; /* out */
709: return(0);
710: }
712: {
713: PetscReal *coords;
714: PetscInt data_stride;
715: PetscInt *crsGID = NULL;
716: Mat Gmat2;
718: if (dim != data_cols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"dim %D != data_cols %D",dim,data_cols);
719: /* grow ghost data for better coarse grid cover of fine grid */
720: #if defined PETSC_GAMG_USE_LOG
721: PetscLogEventBegin(petsc_gamg_setup_events[SET5],0,0,0,0);
722: #endif
723: /* messy method, squares graph and gets some data */
724: getGIDsOnSquareGraph(nLocalSelected, clid_flid, Gmat, &selected_2, &Gmat2, &crsGID);
725: /* llist is now not valid wrt squared graph, but will work as iterator in 'triangulateAndFormProl' */
726: #if defined PETSC_GAMG_USE_LOG
727: PetscLogEventEnd(petsc_gamg_setup_events[SET5],0,0,0,0);
728: #endif
729: /* create global vector of coorindates in 'coords' */
730: if (size > 1) {
731: PCGAMGGetDataWithGhosts(Gmat2, dim, pc_gamg->data, &data_stride, &coords);
732: } else {
733: coords = (PetscReal*)pc_gamg->data;
734: data_stride = pc_gamg->data_sz/pc_gamg->data_cell_cols;
735: }
736: MatDestroy(&Gmat2);
738: /* triangulate */
739: if (dim == 2) {
740: PetscReal metric,tm;
741: #if defined PETSC_GAMG_USE_LOG
742: PetscLogEventBegin(petsc_gamg_setup_events[SET6],0,0,0,0);
743: #endif
744: triangulateAndFormProl(selected_2, data_stride, coords,nLocalSelected, clid_flid, agg_lists, crsGID, bs, Prol, &metric);
745: #if defined PETSC_GAMG_USE_LOG
746: PetscLogEventEnd(petsc_gamg_setup_events[SET6],0,0,0,0);
747: #endif
748: PetscFree(crsGID);
750: /* clean up and create coordinates for coarse grid (output) */
751: if (size > 1) PetscFree(coords);
753: MPIU_Allreduce(&metric, &tm, 1, MPIU_REAL, MPIU_MAX, comm);
754: if (tm > 1.) { /* needs to be globalized - should not happen */
755: PetscInfo1(pc," failed metric for coarse grid %e\n",(double)tm);
756: MatDestroy(&Prol);
757: } else if (metric > .0) {
758: PetscInfo1(pc,"worst metric for coarse grid = %e\n",(double)metric);
759: }
760: } else SETERRQ(comm,PETSC_ERR_PLIB,"3D not implemented for 'geo' AMG");
761: { /* create next coords - output */
762: PetscReal *crs_crds;
763: PetscMalloc1(dim*nLocalSelected, &crs_crds);
764: for (kk=0; kk<nLocalSelected; kk++) { /* grab local select nodes to promote - output */
765: PetscInt lid = clid_flid[kk];
766: for (jj=0; jj<dim; jj++) crs_crds[jj*nLocalSelected + kk] = pc_gamg->data[jj*nloc + lid];
767: }
769: PetscFree(pc_gamg->data);
770: pc_gamg->data = crs_crds; /* out */
771: pc_gamg->data_sz = dim*nLocalSelected;
772: }
773: ISDestroy(&selected_2);
774: }
776: *a_P_out = Prol; /* out */
777: PetscFree(clid_flid);
778: PetscLogEventEnd(PC_GAMGProlongator_GEO,0,0,0,0);
779: return(0);
780: }
782: static PetscErrorCode PCDestroy_GAMG_GEO(PC pc)
783: {
787: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
788: return(0);
789: }
791: /* -------------------------------------------------------------------------- */
792: /*
793: PCCreateGAMG_GEO
795: Input Parameter:
796: . pc -
797: */
798: PetscErrorCode PCCreateGAMG_GEO(PC pc)
799: {
801: PC_MG *mg = (PC_MG*)pc->data;
802: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
805: pc_gamg->ops->setfromoptions = PCSetFromOptions_GEO;
806: pc_gamg->ops->destroy = PCDestroy_GAMG_GEO;
807: /* reset does not do anything; setup not virtual */
809: /* set internal function pointers */
810: pc_gamg->ops->graph = PCGAMGGraph_GEO;
811: pc_gamg->ops->coarsen = PCGAMGCoarsen_GEO;
812: pc_gamg->ops->prolongator = PCGAMGProlongator_GEO;
813: pc_gamg->ops->optprolongator = NULL;
814: pc_gamg->ops->createdefaultdata = PCSetData_GEO;
816: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_GEO);
817: return(0);
818: }