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