1 static char help[] = "VecTagger interface routines.\n\n";
2
3 #include <petscvec.h>
4
ISGetBlockGlobalIS(IS is,Vec vec,PetscInt bs,IS * isBlockGlobal)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;
10
11 PetscFunctionBegin;
12 PetscCall(VecGetLayout(vec, &map));
13 rstart = map->rstart / bs;
14 PetscCall(ISGetLocalSize(is, &n));
15 PetscCall(PetscMalloc1(n, &idxout));
16 PetscCall(ISGetIndices(is, &idxin));
17 for (i = 0; i < n; i++) idxout[i] = rstart + idxin[i];
18 PetscCall(ISRestoreIndices(is, &idxin));
19 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)vec), bs, n, idxout, PETSC_OWN_POINTER, isBlockGlobal));
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
main(int argc,char ** argv)23 int main(int argc, char **argv)
24 {
25 Vec vec, tagged, untagged;
26 VecScatter taggedScatter, untaggedScatter;
27 PetscInt bs;
28 PetscInt n, nloc, nint, i, j, k, localStart, localEnd, ntagged, nuntagged;
29 MPI_Comm comm;
30 VecTagger tagger;
31 PetscScalar *array;
32 PetscRandom rand;
33 VecTaggerBox *defaultBox;
34 VecTaggerBox *boxes;
35 IS is, isBlockGlobal, isComp;
36 PetscBool listed;
37
38 PetscFunctionBeginUser;
39 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
40 n = 10.;
41 bs = 1;
42 comm = PETSC_COMM_WORLD;
43 PetscOptionsBegin(comm, "", "VecTagger Test Options", "Vec");
44 PetscCall(PetscOptionsInt("-bs", "The block size of the vector", "ex1.c", bs, &bs, NULL));
45 PetscCall(PetscOptionsInt("-n", "The size of the vector (in blocks)", "ex1.c", n, &n, NULL));
46 PetscOptionsEnd();
47
48 PetscCall(PetscRandomCreate(comm, &rand));
49 PetscCall(PetscRandomSetFromOptions(rand));
50
51 PetscCall(VecCreate(comm, &vec));
52 PetscCall(PetscObjectSetName((PetscObject)vec, "Vec to Tag"));
53 PetscCall(VecSetBlockSize(vec, bs));
54 PetscCall(VecSetSizes(vec, PETSC_DECIDE, n));
55 PetscCall(VecSetUp(vec));
56 PetscCall(VecGetLocalSize(vec, &nloc));
57 PetscCall(VecGetArray(vec, &array));
58 for (i = 0; i < nloc; i++) {
59 PetscScalar val;
60
61 PetscCall(PetscRandomGetValue(rand, &val));
62 array[i] = val;
63 }
64 PetscCall(VecRestoreArray(vec, &array));
65 PetscCall(PetscRandomDestroy(&rand));
66 PetscCall(VecViewFromOptions(vec, NULL, "-vec_view"));
67
68 PetscCall(VecTaggerCreate(comm, &tagger));
69 PetscCall(VecTaggerSetBlockSize(tagger, bs));
70 PetscCall(VecTaggerSetType(tagger, VECTAGGERABSOLUTE));
71 PetscCall(PetscMalloc1(bs, &defaultBox));
72 for (i = 0; i < bs; i++) {
73 #if !defined(PETSC_USE_COMPLEX)
74 defaultBox[i].min = 0.1;
75 defaultBox[i].max = 1.5;
76 #else
77 defaultBox[i].min = PetscCMPLX(0.1, 0.1);
78 defaultBox[i].max = PetscCMPLX(1.5, 1.5);
79 #endif
80 }
81 PetscCall(VecTaggerAbsoluteSetBox(tagger, defaultBox));
82 PetscCall(PetscFree(defaultBox));
83 PetscCall(VecTaggerSetFromOptions(tagger));
84 PetscCall(VecTaggerSetUp(tagger));
85 PetscCall(PetscObjectViewFromOptions((PetscObject)tagger, NULL, "-vec_tagger_view"));
86 PetscCall(VecTaggerGetBlockSize(tagger, &bs));
87
88 PetscCall(VecTaggerComputeBoxes(tagger, vec, &nint, &boxes, &listed));
89 if (listed) {
90 PetscViewer viewer = NULL;
91
92 PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-vec_tagger_boxes_view", &viewer, NULL, NULL));
93 if (viewer) {
94 PetscBool isascii;
95
96 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
97 if (isascii) {
98 PetscCall(PetscViewerASCIIPrintf(viewer, "Num boxes: %" PetscInt_FMT "\n", nint));
99 PetscCall(PetscViewerASCIIPushTab(viewer));
100 for (i = 0, k = 0; i < nint; i++) {
101 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", i));
102 for (j = 0; j < bs; j++, k++) {
103 if (j) PetscCall(PetscViewerASCIIPrintf(viewer, " x "));
104 #if !defined(PETSC_USE_COMPLEX)
105 PetscCall(PetscViewerASCIIPrintf(viewer, "[%g,%g]", (double)boxes[k].min, (double)boxes[k].max));
106 #else
107 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)));
108 #endif
109 }
110 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
111 }
112 PetscCall(PetscViewerASCIIPopTab(viewer));
113 }
114 }
115 PetscCall(PetscViewerDestroy(&viewer));
116 PetscCall(PetscFree(boxes));
117 }
118
119 PetscCall(VecTaggerComputeIS(tagger, vec, &is, &listed));
120 PetscCall(ISGetBlockGlobalIS(is, vec, bs, &isBlockGlobal));
121 PetscCall(PetscObjectSetName((PetscObject)isBlockGlobal, "Tagged IS (block global)"));
122 PetscCall(ISViewFromOptions(isBlockGlobal, NULL, "-tagged_is_view"));
123
124 PetscCall(VecGetOwnershipRange(vec, &localStart, &localEnd));
125 PetscCall(ISComplement(isBlockGlobal, localStart, localEnd, &isComp));
126 PetscCall(PetscObjectSetName((PetscObject)isComp, "Untagged IS (global)"));
127 PetscCall(ISViewFromOptions(isComp, NULL, "-untagged_is_view"));
128
129 PetscCall(ISGetLocalSize(isBlockGlobal, &ntagged));
130 PetscCall(ISGetLocalSize(isComp, &nuntagged));
131
132 PetscCall(VecCreate(comm, &tagged));
133 PetscCall(PetscObjectSetName((PetscObject)tagged, "Tagged selection"));
134 PetscCall(VecSetSizes(tagged, ntagged, PETSC_DETERMINE));
135 PetscCall(VecSetUp(tagged));
136
137 PetscCall(VecCreate(comm, &untagged));
138 PetscCall(PetscObjectSetName((PetscObject)untagged, "Untagged selection"));
139 PetscCall(VecSetSizes(untagged, nuntagged, PETSC_DETERMINE));
140 PetscCall(VecSetUp(untagged));
141
142 PetscCall(VecScatterCreate(vec, isBlockGlobal, tagged, NULL, &taggedScatter));
143 PetscCall(VecScatterCreate(vec, isComp, untagged, NULL, &untaggedScatter));
144
145 PetscCall(VecScatterBegin(taggedScatter, vec, tagged, INSERT_VALUES, SCATTER_FORWARD));
146 PetscCall(VecScatterEnd(taggedScatter, vec, tagged, INSERT_VALUES, SCATTER_FORWARD));
147 PetscCall(VecScatterBegin(untaggedScatter, vec, untagged, INSERT_VALUES, SCATTER_FORWARD));
148 PetscCall(VecScatterEnd(untaggedScatter, vec, untagged, INSERT_VALUES, SCATTER_FORWARD));
149
150 PetscCall(VecViewFromOptions(tagged, NULL, "-tagged_vec_view"));
151 PetscCall(VecViewFromOptions(untagged, NULL, "-untagged_vec_view"));
152
153 PetscCall(VecScatterDestroy(&untaggedScatter));
154 PetscCall(VecScatterDestroy(&taggedScatter));
155
156 PetscCall(VecDestroy(&untagged));
157 PetscCall(VecDestroy(&tagged));
158 PetscCall(ISDestroy(&isComp));
159 PetscCall(ISDestroy(&isBlockGlobal));
160 PetscCall(ISDestroy(&is));
161
162 PetscCall(VecTaggerDestroy(&tagger));
163 PetscCall(VecDestroy(&vec));
164 PetscCall(PetscFinalize());
165 return 0;
166 }
167
168 /*TEST
169
170 test:
171 suffix: 0
172 requires: !complex
173 args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view
174
175 test:
176 suffix: 1
177 requires: !complex
178 nsize: 3
179 args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view
180
181 test:
182 suffix: 2
183 requires: !complex
184 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
185
186 test:
187 suffix: 3
188 requires: !complex
189 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
190
191 test:
192 suffix: 4
193 requires: !complex
194 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
195
196 test:
197 suffix: 5
198 requires: !complex
199 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
200
201 test:
202 suffix: 6
203 requires: !complex
204 nsize: 3
205 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
206
207 test:
208 suffix: 7
209 requires: !complex
210 nsize: 3
211 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
212
213 test:
214 suffix: 8
215 requires: !complex
216 nsize: 3
217 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
218 filter: sed -e s~Inf~inf~g
219
220 test:
221 suffix: 9
222 requires: !complex
223 nsize: 3
224 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
225 filter: sed -e s~Inf~inf~g
226
227 test:
228 suffix: 10
229 requires: complex
230 args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view
231
232 test:
233 suffix: 11
234 requires: complex
235 nsize: 3
236 args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view
237
238 test:
239 suffix: 12
240 requires: complex
241 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
242
243 test:
244 suffix: 13
245 requires: complex
246 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
247
248 test:
249 suffix: 14
250 requires: complex
251 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
252
253 test:
254 suffix: 15
255 requires: complex
256 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
257
258 test:
259 suffix: 16
260 requires: complex
261 nsize: 3
262 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
263
264 test:
265 suffix: 17
266 requires: complex
267 nsize: 3
268 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
269
270 test:
271 suffix: 18
272 requires: complex
273 nsize: 3
274 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
275 filter: sed -e s~Inf~inf~g
276
277 test:
278 suffix: 19
279 requires: complex
280 nsize: 3
281 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
282 filter: sed -e s~Inf~inf~g
283
284 TEST*/
285