Actual source code: ex1.c

  1: static char help[] = "VecTagger interface routines.\n\n";

  3: #include <petscvec.h>

  5: static PetscErrorCode ISGetBlockGlobalIS(IS is, Vec vec, PetscInt bs, IS *isBlockGlobal)
  6: {
  7:   const PetscInt *idxin;
  8:   PetscInt       *idxout, i, n, rstart;
  9:   PetscLayout     map;

 11:   PetscFunctionBegin;

 13:   PetscCall(VecGetLayout(vec, &map));
 14:   rstart = map->rstart / bs;
 15:   PetscCall(ISGetLocalSize(is, &n));
 16:   PetscCall(PetscMalloc1(n, &idxout));
 17:   PetscCall(ISGetIndices(is, &idxin));
 18:   for (i = 0; i < n; i++) idxout[i] = rstart + idxin[i];
 19:   PetscCall(ISRestoreIndices(is, &idxin));
 20:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)vec), bs, n, idxout, PETSC_OWN_POINTER, isBlockGlobal));
 21:   PetscFunctionReturn(PETSC_SUCCESS);
 22: }

 24: int main(int argc, char **argv)
 25: {
 26:   Vec           vec, tagged, untagged;
 27:   VecScatter    taggedScatter, untaggedScatter;
 28:   PetscInt      bs;
 29:   PetscInt      n, nloc, nint, i, j, k, localStart, localEnd, ntagged, nuntagged;
 30:   MPI_Comm      comm;
 31:   VecTagger     tagger;
 32:   PetscScalar  *array;
 33:   PetscRandom   rand;
 34:   VecTaggerBox *defaultBox;
 35:   VecTaggerBox *boxes;
 36:   IS            is, isBlockGlobal, isComp;
 37:   PetscBool     listed;

 39:   PetscFunctionBeginUser;
 40:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 41:   n    = 10.;
 42:   bs   = 1;
 43:   comm = PETSC_COMM_WORLD;
 44:   PetscOptionsBegin(comm, "", "VecTagger Test Options", "Vec");
 45:   PetscCall(PetscOptionsInt("-bs", "The block size of the vector", "ex1.c", bs, &bs, NULL));
 46:   PetscCall(PetscOptionsInt("-n", "The size of the vector (in blocks)", "ex1.c", n, &n, NULL));
 47:   PetscOptionsEnd();

 49:   PetscCall(PetscRandomCreate(comm, &rand));
 50:   PetscCall(PetscRandomSetFromOptions(rand));

 52:   PetscCall(VecCreate(comm, &vec));
 53:   PetscCall(PetscObjectSetName((PetscObject)vec, "Vec to Tag"));
 54:   PetscCall(VecSetBlockSize(vec, bs));
 55:   PetscCall(VecSetSizes(vec, PETSC_DECIDE, n));
 56:   PetscCall(VecSetUp(vec));
 57:   PetscCall(VecGetLocalSize(vec, &nloc));
 58:   PetscCall(VecGetArray(vec, &array));
 59:   for (i = 0; i < nloc; i++) {
 60:     PetscScalar val;

 62:     PetscCall(PetscRandomGetValue(rand, &val));
 63:     array[i] = val;
 64:   }
 65:   PetscCall(VecRestoreArray(vec, &array));
 66:   PetscCall(PetscRandomDestroy(&rand));
 67:   PetscCall(VecViewFromOptions(vec, NULL, "-vec_view"));

 69:   PetscCall(VecTaggerCreate(comm, &tagger));
 70:   PetscCall(VecTaggerSetBlockSize(tagger, bs));
 71:   PetscCall(VecTaggerSetType(tagger, VECTAGGERABSOLUTE));
 72:   PetscCall(PetscMalloc1(bs, &defaultBox));
 73:   for (i = 0; i < bs; i++) {
 74: #if !defined(PETSC_USE_COMPLEX)
 75:     defaultBox[i].min = 0.1;
 76:     defaultBox[i].max = 1.5;
 77: #else
 78:     defaultBox[i].min = PetscCMPLX(0.1, 0.1);
 79:     defaultBox[i].max = PetscCMPLX(1.5, 1.5);
 80: #endif
 81:   }
 82:   PetscCall(VecTaggerAbsoluteSetBox(tagger, defaultBox));
 83:   PetscCall(PetscFree(defaultBox));
 84:   PetscCall(VecTaggerSetFromOptions(tagger));
 85:   PetscCall(VecTaggerSetUp(tagger));
 86:   PetscCall(PetscObjectViewFromOptions((PetscObject)tagger, NULL, "-vec_tagger_view"));
 87:   PetscCall(VecTaggerGetBlockSize(tagger, &bs));

 89:   PetscCall(VecTaggerComputeBoxes(tagger, vec, &nint, &boxes, &listed));
 90:   if (listed) {
 91:     PetscViewer viewer = NULL;

 93:     PetscCall(PetscOptionsGetViewer(comm, NULL, NULL, "-vec_tagger_boxes_view", &viewer, NULL, NULL));
 94:     if (viewer) {
 95:       PetscBool iascii;

 97:       PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 98:       if (iascii) {
 99:         PetscCall(PetscViewerASCIIPrintf(viewer, "Num boxes: %" PetscInt_FMT "\n", nint));
100:         PetscCall(PetscViewerASCIIPushTab(viewer));
101:         for (i = 0, k = 0; i < nint; i++) {
102:           PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", i));
103:           for (j = 0; j < bs; j++, k++) {
104:             if (j) PetscCall(PetscViewerASCIIPrintf(viewer, " x "));
105: #if !defined(PETSC_USE_COMPLEX)
106:             PetscCall(PetscViewerASCIIPrintf(viewer, "[%g,%g]", (double)boxes[k].min, (double)boxes[k].max));
107: #else
108:             PetscCall(PetscViewerASCIIPrintf(viewer, "[%g+%gi,%g+%gi]", (double)PetscRealPart(boxes[k].min), (double)PetscImaginaryPart(boxes[k].min), (double)PetscRealPart(boxes[k].max), (double)PetscImaginaryPart(boxes[k].max)));
109: #endif
110:           }
111:           PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
112:         }
113:         PetscCall(PetscViewerASCIIPopTab(viewer));
114:       }
115:     }
116:     PetscCall(PetscViewerDestroy(&viewer));
117:     PetscCall(PetscFree(boxes));
118:   }

120:   PetscCall(VecTaggerComputeIS(tagger, vec, &is, &listed));
121:   PetscCall(ISGetBlockGlobalIS(is, vec, bs, &isBlockGlobal));
122:   PetscCall(PetscObjectSetName((PetscObject)isBlockGlobal, "Tagged IS (block global)"));
123:   PetscCall(ISViewFromOptions(isBlockGlobal, NULL, "-tagged_is_view"));

125:   PetscCall(VecGetOwnershipRange(vec, &localStart, &localEnd));
126:   PetscCall(ISComplement(isBlockGlobal, localStart, localEnd, &isComp));
127:   PetscCall(PetscObjectSetName((PetscObject)isComp, "Untagged IS (global)"));
128:   PetscCall(ISViewFromOptions(isComp, NULL, "-untagged_is_view"));

130:   PetscCall(ISGetLocalSize(isBlockGlobal, &ntagged));
131:   PetscCall(ISGetLocalSize(isComp, &nuntagged));

133:   PetscCall(VecCreate(comm, &tagged));
134:   PetscCall(PetscObjectSetName((PetscObject)tagged, "Tagged selection"));
135:   PetscCall(VecSetSizes(tagged, ntagged, PETSC_DETERMINE));
136:   PetscCall(VecSetUp(tagged));

138:   PetscCall(VecCreate(comm, &untagged));
139:   PetscCall(PetscObjectSetName((PetscObject)untagged, "Untagged selection"));
140:   PetscCall(VecSetSizes(untagged, nuntagged, PETSC_DETERMINE));
141:   PetscCall(VecSetUp(untagged));

143:   PetscCall(VecScatterCreate(vec, isBlockGlobal, tagged, NULL, &taggedScatter));
144:   PetscCall(VecScatterCreate(vec, isComp, untagged, NULL, &untaggedScatter));

146:   PetscCall(VecScatterBegin(taggedScatter, vec, tagged, INSERT_VALUES, SCATTER_FORWARD));
147:   PetscCall(VecScatterEnd(taggedScatter, vec, tagged, INSERT_VALUES, SCATTER_FORWARD));
148:   PetscCall(VecScatterBegin(untaggedScatter, vec, untagged, INSERT_VALUES, SCATTER_FORWARD));
149:   PetscCall(VecScatterEnd(untaggedScatter, vec, untagged, INSERT_VALUES, SCATTER_FORWARD));

151:   PetscCall(VecViewFromOptions(tagged, NULL, "-tagged_vec_view"));
152:   PetscCall(VecViewFromOptions(untagged, NULL, "-untagged_vec_view"));

154:   PetscCall(VecScatterDestroy(&untaggedScatter));
155:   PetscCall(VecScatterDestroy(&taggedScatter));

157:   PetscCall(VecDestroy(&untagged));
158:   PetscCall(VecDestroy(&tagged));
159:   PetscCall(ISDestroy(&isComp));
160:   PetscCall(ISDestroy(&isBlockGlobal));
161:   PetscCall(ISDestroy(&is));

163:   PetscCall(VecTaggerDestroy(&tagger));
164:   PetscCall(VecDestroy(&vec));
165:   PetscCall(PetscFinalize());
166:   return 0;
167: }

169: /*TEST

171:   test:
172:     suffix: 0
173:     requires: !complex
174:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

176:   test:
177:     suffix: 1
178:     requires: !complex
179:     nsize: 3
180:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

182:   test:
183:     suffix: 2
184:     requires: !complex
185:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -bs 2

187:   test:
188:     suffix: 3
189:     requires: !complex
190:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_block_size 2 -vec_tagger_box 0.1,1.5,0.1,1.5

192:   test:
193:     suffix: 4
194:     requires: !complex
195:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_invert

197:   test:
198:     suffix: 5
199:     requires: !complex
200:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type relative -vec_tagger_box 0.25,0.75

202:   test:
203:     suffix: 6
204:     requires: !complex
205:     nsize: 3
206:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25,0.75

208:   test:
209:     suffix: 7
210:     requires: !complex
211:     nsize: 3
212:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25,0.75 -vec_tagger_cdf_method iterative -vec_tagger_cdf_max_it 10

214:   test:
215:     suffix: 8
216:     requires: !complex
217:     nsize: 3
218:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type or -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box 0.0,0.25 -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.75,inf
219:     filter: sed -e s~Inf~inf~g

221:   test:
222:     suffix: 9
223:     requires: !complex
224:     nsize: 3
225:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type and -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box -inf,0.5 -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.25,0.75
226:     filter: sed -e s~Inf~inf~g

228:   test:
229:     suffix: 10
230:     requires: complex
231:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

233:   test:
234:     suffix: 11
235:     requires: complex
236:     nsize: 3
237:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

239:   test:
240:     suffix: 12
241:     requires: complex
242:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -bs 2

244:   test:
245:     suffix: 13
246:     requires: complex
247:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_block_size 2 -vec_tagger_box 0.1+0.1i,1.5+1.5i,0.1+0.1i,1.5+1.5i

249:   test:
250:     suffix: 14
251:     requires: complex
252:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_invert

254:   test:
255:     suffix: 15
256:     requires: complex
257:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type relative -vec_tagger_box 0.25+0.25i,0.75+0.75i

259:   test:
260:     suffix: 16
261:     requires: complex
262:     nsize: 3
263:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25+0.25i,0.75+0.75i

265:   test:
266:     suffix: 17
267:     requires: complex
268:     nsize: 3
269:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25+0.25i,0.75+0.75i -vec_tagger_cdf_method iterative -vec_tagger_cdf_max_it 10

271:   test:
272:     suffix: 18
273:     requires: complex
274:     nsize: 3
275:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type or -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box 0.0+0.0i,0.25+0.25i -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.75+0.75i,inf+infi
276:     filter: sed -e s~Inf~inf~g

278:   test:
279:     suffix: 19
280:     requires: complex
281:     nsize: 3
282:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type and -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box -inf-infi,0.75+0.75i -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.25+0.25i,1.+1.i
283:     filter: sed -e s~Inf~inf~g

285: TEST*/