Actual source code: ex18.c

  1: static char help[] = "Test PetscSFConcatenate()\n\n";

  3: #include <petscsf.h>

  5: typedef struct {
  6:   MPI_Comm                   comm;
  7:   PetscMPIInt                rank, size;
  8:   PetscInt                   leaveStep, nsfs, n;
  9:   PetscBool                  sparseLeaves;
 10:   PetscBool                  compare;
 11:   PetscBool                  irregular;
 12:   PetscSFConcatenateRootMode rootMode;
 13:   PetscViewer                viewer;
 14: } AppCtx;

 16: static PetscErrorCode GetOptions(MPI_Comm comm, AppCtx *ctx)
 17: {
 18:   PetscViewerFormat format;

 20:   PetscFunctionBegin;
 21:   ctx->comm         = comm;
 22:   ctx->nsfs         = 3;
 23:   ctx->n            = 1;
 24:   ctx->leaveStep    = 1;
 25:   ctx->sparseLeaves = PETSC_FALSE;
 26:   ctx->compare      = PETSC_FALSE;
 27:   ctx->irregular    = PETSC_FALSE;
 28:   ctx->rootMode     = PETSCSF_CONCATENATE_ROOTMODE_LOCAL;
 29:   ctx->viewer       = NULL;
 30:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nsfs", &ctx->nsfs, NULL));
 31:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &ctx->n, NULL));
 32:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-leave_step", &ctx->leaveStep, NULL));
 33:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-irregular", &ctx->irregular, NULL));
 34:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-compare_to_reference", &ctx->compare, NULL));
 35:   PetscCall(PetscOptionsGetEnum(NULL, NULL, "-root_mode", PetscSFConcatenateRootModes, (PetscEnum *)&ctx->rootMode, NULL));
 36:   PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-sf_view", &ctx->viewer, &format, NULL));
 37:   if (ctx->viewer) PetscCall(PetscViewerPushFormat(ctx->viewer, format));
 38:   ctx->sparseLeaves = (PetscBool)(ctx->leaveStep != 1);
 39:   PetscCallMPI(MPI_Comm_size(comm, &ctx->size));
 40:   PetscCallMPI(MPI_Comm_rank(comm, &ctx->rank));
 41:   PetscFunctionReturn(PETSC_SUCCESS);
 42: }

 44: static PetscErrorCode PetscSFCheckEqual_Private(PetscSF sf0, PetscSF sf1)
 45: {
 46:   PetscInt  nRoot, nLeave;
 47:   Vec       vecRoot0, vecLeave0, vecRoot1, vecLeave1;
 48:   MPI_Comm  comm;
 49:   PetscBool flg;

 51:   PetscFunctionBegin;
 52:   PetscCall(PetscObjectGetComm((PetscObject)sf0, &comm));
 53:   PetscCall(PetscSFGetGraph(sf0, &nRoot, NULL, NULL, NULL));
 54:   PetscCall(PetscSFGetLeafRange(sf0, NULL, &nLeave));
 55:   nLeave++;
 56:   PetscCall(VecCreateFromOptions(comm, NULL, 1, nRoot, PETSC_DECIDE, &vecRoot0));
 57:   PetscCall(VecCreateFromOptions(comm, NULL, 1, nLeave, PETSC_DECIDE, &vecLeave0));
 58:   PetscCall(VecDuplicate(vecRoot0, &vecRoot1));
 59:   PetscCall(VecDuplicate(vecLeave0, &vecLeave1));
 60:   {
 61:     PetscRandom rand;

 63:     PetscCall(PetscRandomCreate(comm, &rand));
 64:     PetscCall(PetscRandomSetFromOptions(rand));
 65:     PetscCall(VecSetRandom(vecRoot0, rand));
 66:     PetscCall(VecSetRandom(vecLeave0, rand));
 67:     PetscCall(VecCopy(vecRoot0, vecRoot1));
 68:     PetscCall(VecCopy(vecLeave0, vecLeave1));
 69:     PetscCall(PetscRandomDestroy(&rand));
 70:   }

 72:   PetscCall(VecScatterBegin(sf0, vecRoot0, vecLeave0, ADD_VALUES, SCATTER_FORWARD));
 73:   PetscCall(VecScatterEnd(sf0, vecRoot0, vecLeave0, ADD_VALUES, SCATTER_FORWARD));
 74:   PetscCall(VecScatterBegin(sf1, vecRoot1, vecLeave1, ADD_VALUES, SCATTER_FORWARD));
 75:   PetscCall(VecScatterEnd(sf1, vecRoot1, vecLeave1, ADD_VALUES, SCATTER_FORWARD));
 76:   PetscCall(VecEqual(vecLeave0, vecLeave1, &flg));
 77:   PetscCheck(flg, comm, PETSC_ERR_PLIB, "leave vectors differ");

 79:   PetscCall(VecScatterBegin(sf0, vecLeave0, vecRoot0, ADD_VALUES, SCATTER_REVERSE));
 80:   PetscCall(VecScatterEnd(sf0, vecLeave0, vecRoot0, ADD_VALUES, SCATTER_REVERSE));
 81:   PetscCall(VecScatterBegin(sf1, vecLeave1, vecRoot1, ADD_VALUES, SCATTER_REVERSE));
 82:   PetscCall(VecScatterEnd(sf1, vecLeave1, vecRoot1, ADD_VALUES, SCATTER_REVERSE));
 83:   PetscCall(VecEqual(vecRoot0, vecRoot1, &flg));
 84:   PetscCheck(flg, comm, PETSC_ERR_PLIB, "root vectors differ");

 86:   PetscCall(VecDestroy(&vecRoot0));
 87:   PetscCall(VecDestroy(&vecRoot1));
 88:   PetscCall(VecDestroy(&vecLeave0));
 89:   PetscCall(VecDestroy(&vecLeave1));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: static PetscErrorCode PetscSFViewCustom(PetscSF sf, PetscViewer viewer)
 94: {
 95:   PetscMPIInt        rank, nranks;
 96:   PetscInt           i, nroots, nleaves;
 97:   const PetscInt    *ilocal;
 98:   const PetscSFNode *iremote;
 99:   PetscLayout        rootLayout;
100:   PetscInt          *gremote;

102:   PetscFunctionBegin;
103:   PetscCall(PetscSFSetUp(sf));
104:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, &iremote));
105:   PetscCall(PetscSFGetRootRanks(sf, &nranks, NULL, NULL, NULL, NULL));
106:   PetscCall(PetscSFGetGraphLayout(sf, &rootLayout, NULL, NULL, &gremote));
107:   PetscCheck(nroots == rootLayout->n, PetscObjectComm((PetscObject)sf), PETSC_ERR_PLIB, "Assertion failed: nroots == rootLayout->n");
108:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sf, viewer));
109:   PetscCall(PetscViewerASCIIPushTab(viewer));
110:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sf), &rank));
111:   PetscCall(PetscViewerASCIIPushSynchronized(viewer));
112:   if (rank == 0) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "rank #leaves #roots\n"));
113:   PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%2d] %7" PetscInt_FMT " %6" PetscInt_FMT "\n", rank, nleaves, nroots));
114:   PetscCall(PetscViewerFlush(viewer));
115:   if (rank == 0) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "leaves      roots       roots in global numbering\n"));
116:   for (i = 0; i < nleaves; i++)
117:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "(%2d, %2d) <- (%2d, %2" PetscInt_FMT ")  = %3" PetscInt_FMT "\n", rank, (PetscMPIInt)(ilocal ? ilocal[i] : i), (PetscMPIInt)iremote[i].rank, iremote[i].index, gremote[i]));
118:   PetscCall(PetscViewerFlush(viewer));
119:   PetscCall(PetscViewerASCIIPopSynchronized(viewer));
120:   PetscCall(PetscViewerASCIIPopTab(viewer));
121:   PetscCall(PetscLayoutDestroy(&rootLayout));
122:   PetscCall(PetscFree(gremote));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: PetscErrorCode CreateReferenceSF_Regular(AppCtx *ctx, PetscSF *refSF)
127: {
128:   PetscInt  j;
129:   PetscInt *ilocal  = NULL;
130:   PetscInt  nLeaves = ctx->nsfs * ctx->n * ctx->size;
131:   PetscInt  nroots  = ctx->n * ctx->nsfs;
132:   PetscSF   sf;

134:   PetscFunctionBegin;
135:   ilocal = NULL;
136:   if (ctx->sparseLeaves) PetscCall(PetscCalloc1(nLeaves + 1, &ilocal));
137:   PetscCall(PetscSFCreate(ctx->comm, &sf));
138:   for (j = 0; j < nLeaves; j++) {
139:     if (ctx->sparseLeaves) ilocal[j + 1] = ilocal[j] + ctx->leaveStep;
140:   }
141:   switch (ctx->rootMode) {
142:   case PETSCSF_CONCATENATE_ROOTMODE_SHARED:
143:   case PETSCSF_CONCATENATE_ROOTMODE_LOCAL: {
144:     PetscInt     i, k;
145:     PetscMPIInt  r;
146:     PetscSFNode *iremote;

148:     PetscCall(PetscCalloc1(nLeaves, &iremote));
149:     for (i = 0, j = 0; i < ctx->nsfs; i++) {
150:       for (r = 0; r < ctx->size; r++) {
151:         for (k = 0; k < ctx->n; k++, j++) {
152:           iremote[j].rank  = r;
153:           iremote[j].index = k + i * ctx->n;
154:         }
155:       }
156:     }
157:     PetscCall(PetscSFSetGraph(sf, nroots, nLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
158:   } break;
159:   case PETSCSF_CONCATENATE_ROOTMODE_GLOBAL: {
160:     PetscLayout map = NULL;
161:     PetscInt   *gremote;

163:     PetscCall(PetscLayoutCreateFromSizes(ctx->comm, nroots, PETSC_DECIDE, 1, &map));
164:     PetscCall(PetscMalloc1(nLeaves, &gremote));
165:     for (j = 0; j < nLeaves; j++) gremote[j] = j;
166:     PetscCall(PetscSFSetGraphLayout(sf, map, nLeaves, ilocal, PETSC_OWN_POINTER, gremote));
167:     PetscCall(PetscFree(gremote));
168:     PetscCall(PetscLayoutDestroy(&map));
169:   } break;
170:   default:
171:     SETERRQ(ctx->comm, PETSC_ERR_SUP, "unsupported rootmode %d", ctx->rootMode);
172:   }
173:   PetscCall(PetscObjectSetName((PetscObject)sf, "reference_sf"));
174:   if (ctx->viewer) PetscCall(PetscSFViewCustom(sf, ctx->viewer));
175:   *refSF = sf;
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: PetscErrorCode CreateSFs_Irregular(AppCtx *ctx, PetscSF *newSFs[], PetscInt *leafOffsets[])
180: {
181:   PetscInt *lOffsets = NULL;
182:   PetscSF  *sfs;
183:   PetscInt  nLeaves = ctx->n * ctx->size + (ctx->size - 1) * ctx->size / 2;
184:   PetscInt  nroots  = ctx->n + ctx->rank + ctx->nsfs - 1 + ctx->size - 1;

186:   PetscFunctionBegin;
187:   if (ctx->sparseLeaves) PetscCall(PetscCalloc1(ctx->nsfs + 1, &lOffsets));
188:   PetscCall(PetscMalloc1(ctx->nsfs, &sfs));
189:   for (PetscInt i = 0; i < ctx->nsfs; i++) {
190:     PetscSF      sf;
191:     PetscInt     j, k;
192:     PetscMPIInt  r;
193:     PetscInt    *ilocal = NULL;
194:     PetscSFNode *iremote;
195:     char         name[32];

197:     if (ctx->sparseLeaves) PetscCall(PetscCalloc1(nLeaves + 1, &ilocal));
198:     PetscCall(PetscMalloc1(nLeaves, &iremote));
199:     for (r = ctx->size - 1, j = 0; r >= 0; r--) {
200:       for (k = 0; k < ctx->n + r; k++, j++) {
201:         if (ctx->sparseLeaves) ilocal[j + 1] = ilocal[j] + ctx->leaveStep;
202:         iremote[j].rank  = r;
203:         iremote[j].index = k + i + ctx->rank;
204:       }
205:     }
206:     if (ctx->sparseLeaves) lOffsets[i + 1] = lOffsets[i] + ilocal[j];

208:     PetscCall(PetscSFCreate(ctx->comm, &sf));
209:     PetscCall(PetscSFSetGraph(sf, nroots, nLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
210:     PetscCall(PetscSNPrintf(name, sizeof(name), "sf_%" PetscInt_FMT, i));
211:     PetscCall(PetscObjectSetName((PetscObject)sf, name));
212:     if (ctx->viewer) PetscCall(PetscSFViewCustom(sf, ctx->viewer));
213:     sfs[i] = sf;
214:   }
215:   *newSFs      = sfs;
216:   *leafOffsets = lOffsets;
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: PetscErrorCode CreateSFs_Regular(AppCtx *ctx, PetscSF *newSFs[], PetscInt *leafOffsets[])
221: {
222:   PetscInt                  *lOffsets = NULL;
223:   PetscInt                   nLeaves  = ctx->n * ctx->size;
224:   PetscSF                   *sfs;
225:   PetscSFConcatenateRootMode mode = ctx->compare ? ctx->rootMode : PETSCSF_CONCATENATE_ROOTMODE_LOCAL;

227:   PetscFunctionBegin;
228:   if (ctx->sparseLeaves) PetscCall(PetscCalloc1(ctx->nsfs + 1, &lOffsets));
229:   PetscCall(PetscCalloc1(ctx->nsfs, &sfs));
230:   for (PetscInt i = 0; i < ctx->nsfs; i++) {
231:     PetscSF   sf;
232:     PetscInt  j;
233:     PetscInt *ilocal = NULL;
234:     char      name[32];

236:     PetscCall(PetscSFCreate(ctx->comm, &sf));
237:     if (ctx->sparseLeaves) {
238:       PetscCall(PetscCalloc1(nLeaves + 1, &ilocal));
239:       for (j = 0; j < nLeaves; j++) ilocal[j + 1] = ilocal[j] + ctx->leaveStep;
240:       lOffsets[i + 1] = lOffsets[i] + ilocal[nLeaves];
241:     }
242:     switch (mode) {
243:     case PETSCSF_CONCATENATE_ROOTMODE_LOCAL: {
244:       PetscInt     k, nroots = ctx->n;
245:       PetscMPIInt  r;
246:       PetscSFNode *iremote;

248:       PetscCall(PetscMalloc1(nLeaves, &iremote));
249:       for (r = 0, j = 0; r < ctx->size; r++) {
250:         for (k = 0; k < ctx->n; k++, j++) {
251:           iremote[j].rank  = r;
252:           iremote[j].index = k;
253:         }
254:       }
255:       PetscCall(PetscSFSetGraph(sf, nroots, nLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
256:     } break;
257:     case PETSCSF_CONCATENATE_ROOTMODE_SHARED: {
258:       PetscInt     k, nroots = ctx->n * ctx->nsfs;
259:       PetscMPIInt  r;
260:       PetscSFNode *iremote;

262:       PetscCall(PetscMalloc1(nLeaves, &iremote));
263:       for (r = 0, j = 0; r < ctx->size; r++) {
264:         for (k = 0; k < ctx->n; k++, j++) {
265:           iremote[j].rank  = r;
266:           iremote[j].index = k + i * ctx->n;
267:         }
268:       }
269:       PetscCall(PetscSFSetGraph(sf, nroots, nLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
270:     } break;
271:     case PETSCSF_CONCATENATE_ROOTMODE_GLOBAL: {
272:       PetscInt    nroots = ctx->n;
273:       PetscLayout map    = NULL;
274:       PetscInt   *gremote;

276:       PetscCall(PetscLayoutCreateFromSizes(ctx->comm, nroots, PETSC_DECIDE, 1, &map));
277:       PetscCall(PetscMalloc1(nLeaves, &gremote));
278:       for (j = 0; j < nLeaves; j++) gremote[j] = j;
279:       PetscCall(PetscSFSetGraphLayout(sf, map, nLeaves, ilocal, PETSC_OWN_POINTER, gremote));
280:       PetscCall(PetscFree(gremote));
281:       PetscCall(PetscLayoutDestroy(&map));
282:     } break;
283:     default:
284:       SETERRQ(ctx->comm, PETSC_ERR_SUP, "unsupported rootmode %d", ctx->rootMode);
285:     }
286:     PetscCall(PetscSNPrintf(name, sizeof(name), "sf_%" PetscInt_FMT, i));
287:     PetscCall(PetscObjectSetName((PetscObject)sf, name));
288:     if (ctx->viewer) PetscCall(PetscSFViewCustom(sf, ctx->viewer));
289:     sfs[i] = sf;
290:   }
291:   *newSFs      = sfs;
292:   *leafOffsets = lOffsets;
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: PetscErrorCode DestroySFs(AppCtx *ctx, PetscSF *sfs[])
297: {
298:   PetscFunctionBegin;
299:   for (PetscInt i = 0; i < ctx->nsfs; i++) PetscCall(PetscSFDestroy(&(*sfs)[i]));
300:   PetscCall(PetscFree(*sfs));
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: int main(int argc, char **argv)
305: {
306:   AppCtx    ctx_;
307:   AppCtx   *ctx = &ctx_;
308:   PetscSF   sf;
309:   PetscSF  *sfs         = NULL;
310:   PetscInt *leafOffsets = NULL;
311:   MPI_Comm  comm;

313:   PetscFunctionBeginUser;
314:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
315:   comm = PETSC_COMM_WORLD;
316:   PetscCall(GetOptions(comm, ctx));

318:   if (ctx->irregular) {
319:     PetscCall(CreateSFs_Irregular(ctx, &sfs, &leafOffsets));
320:   } else {
321:     PetscCall(CreateSFs_Regular(ctx, &sfs, &leafOffsets));
322:   }
323:   PetscCall(PetscSFConcatenate(comm, ctx->nsfs, sfs, ctx->rootMode, leafOffsets, &sf));
324:   PetscCall(PetscObjectSetName((PetscObject)sf, "result_sf"));
325:   if (ctx->viewer) {
326:     PetscCall(PetscPrintf(comm, "rootMode = %s:\n", PetscSFConcatenateRootModes[ctx->rootMode]));
327:     PetscCall(PetscSFViewCustom(sf, ctx->viewer));
328:   }
329:   if (ctx->compare) {
330:     PetscSF sfRef;

332:     PetscAssert(!ctx->irregular, comm, PETSC_ERR_SUP, "Combination  -compare_to_reference true -irregular true  not implemented");
333:     PetscCall(CreateReferenceSF_Regular(ctx, &sfRef));
334:     PetscCall(PetscSFCheckEqual_Private(sf, sfRef));
335:     PetscCall(PetscSFDestroy(&sfRef));
336:   }
337:   PetscCall(DestroySFs(ctx, &sfs));
338:   PetscCall(PetscFree(leafOffsets));
339:   PetscCall(PetscSFDestroy(&sf));
340:   if (ctx->viewer) {
341:     PetscCall(PetscViewerPopFormat(ctx->viewer));
342:     PetscCall(PetscViewerDestroy(&ctx->viewer));
343:   }
344:   PetscCall(PetscFinalize());
345:   return 0;
346: }

348: /*TEST
349:   test:
350:     nsize: {{1 3}}
351:     args: -compare_to_reference -nsfs {{1 3}} -n {{0 1 5}} -leave_step {{1 3}} -root_mode {{local shared global}}
352:     output_file: output/empty.out

354:   test:
355:     suffix: 2
356:     nsize: 2
357:     args: -irregular {{false true}separate output} -sf_view -nsfs 3 -n 1 -leave_step {{1 3}separate output} -root_mode {{local shared global}separate output}
358: TEST*/