1 #include <petsc/private/vecimpl.h> /*I "petscvec.h" I*/
2
3 /*@
4 VecWhichEqual - Creates an index set containing the indices
5 where the vectors `Vec1` and `Vec2` have identical elements.
6
7 Collective
8
9 Input Parameters:
10 + Vec1 - the first vector to compare
11 - Vec2 - the second two vector to compare
12
13 Output Parameter:
14 . S - The index set containing the indices i where vec1[i] == vec2[i]
15
16 Level: advanced
17
18 Note:
19 The two vectors must have the same parallel layout
20
21 .seealso: `Vec`
22 @*/
VecWhichEqual(Vec Vec1,Vec Vec2,IS * S)23 PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS *S)
24 {
25 PetscInt i, n_same = 0;
26 PetscInt n, low, high;
27 PetscInt *same = NULL;
28 const PetscScalar *v1, *v2;
29
30 PetscFunctionBegin;
31 PetscValidHeaderSpecific(Vec1, VEC_CLASSID, 1);
32 PetscValidHeaderSpecific(Vec2, VEC_CLASSID, 2);
33 PetscCheckSameComm(Vec1, 1, Vec2, 2);
34 VecCheckSameSize(Vec1, 1, Vec2, 2);
35
36 PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
37 PetscCall(VecGetLocalSize(Vec1, &n));
38 if (n > 0) {
39 if (Vec1 == Vec2) {
40 PetscCall(VecGetArrayRead(Vec1, &v1));
41 v2 = v1;
42 } else {
43 PetscCall(VecGetArrayRead(Vec1, &v1));
44 PetscCall(VecGetArrayRead(Vec2, &v2));
45 }
46
47 PetscCall(PetscMalloc1(n, &same));
48
49 for (i = 0; i < n; ++i) {
50 if (v1[i] == v2[i]) {
51 same[n_same] = low + i;
52 ++n_same;
53 }
54 }
55
56 if (Vec1 == Vec2) {
57 PetscCall(VecRestoreArrayRead(Vec1, &v1));
58 } else {
59 PetscCall(VecRestoreArrayRead(Vec1, &v1));
60 PetscCall(VecRestoreArrayRead(Vec2, &v2));
61 }
62 }
63 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_same, same, PETSC_OWN_POINTER, S));
64 PetscFunctionReturn(PETSC_SUCCESS);
65 }
66
67 /*@
68 VecWhichLessThan - Creates an index set containing the indices
69 where the vectors `Vec1` < `Vec2`
70
71 Collective
72
73 Input Parameters:
74 + Vec1 - the first vector to compare
75 - Vec2 - the second vector to compare
76
77 Output Parameter:
78 . S - The index set containing the indices i where vec1[i] < vec2[i]
79
80 Level: advanced
81
82 Notes:
83 The two vectors must have the same parallel layout
84
85 For complex numbers this only compares the real part
86
87 .seealso: `Vec`
88 @*/
VecWhichLessThan(Vec Vec1,Vec Vec2,IS * S)89 PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
90 {
91 PetscInt i, n_lt = 0;
92 PetscInt n, low, high;
93 PetscInt *lt = NULL;
94 const PetscScalar *v1, *v2;
95
96 PetscFunctionBegin;
97 PetscValidHeaderSpecific(Vec1, VEC_CLASSID, 1);
98 PetscValidHeaderSpecific(Vec2, VEC_CLASSID, 2);
99 PetscCheckSameComm(Vec1, 1, Vec2, 2);
100 VecCheckSameSize(Vec1, 1, Vec2, 2);
101
102 PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
103 PetscCall(VecGetLocalSize(Vec1, &n));
104 if (n > 0) {
105 if (Vec1 == Vec2) {
106 PetscCall(VecGetArrayRead(Vec1, &v1));
107 v2 = v1;
108 } else {
109 PetscCall(VecGetArrayRead(Vec1, &v1));
110 PetscCall(VecGetArrayRead(Vec2, &v2));
111 }
112
113 PetscCall(PetscMalloc1(n, <));
114
115 for (i = 0; i < n; ++i) {
116 if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {
117 lt[n_lt] = low + i;
118 ++n_lt;
119 }
120 }
121
122 if (Vec1 == Vec2) {
123 PetscCall(VecRestoreArrayRead(Vec1, &v1));
124 } else {
125 PetscCall(VecRestoreArrayRead(Vec1, &v1));
126 PetscCall(VecRestoreArrayRead(Vec2, &v2));
127 }
128 }
129 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_lt, lt, PETSC_OWN_POINTER, S));
130 PetscFunctionReturn(PETSC_SUCCESS);
131 }
132
133 /*@
134 VecWhichGreaterThan - Creates an index set containing the indices
135 where the vectors `Vec1` > `Vec2`
136
137 Collective
138
139 Input Parameters:
140 + Vec1 - the first vector to compare
141 - Vec2 - the second vector to compare
142
143 Output Parameter:
144 . S - The index set containing the indices i where vec1[i] > vec2[i]
145
146 Level: advanced
147
148 Notes:
149 The two vectors must have the same parallel layout
150
151 For complex numbers this only compares the real part
152
153 .seealso: `Vec`
154 @*/
VecWhichGreaterThan(Vec Vec1,Vec Vec2,IS * S)155 PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
156 {
157 PetscInt i, n_gt = 0;
158 PetscInt n, low, high;
159 PetscInt *gt = NULL;
160 const PetscScalar *v1, *v2;
161
162 PetscFunctionBegin;
163 PetscValidHeaderSpecific(Vec1, VEC_CLASSID, 1);
164 PetscValidHeaderSpecific(Vec2, VEC_CLASSID, 2);
165 PetscCheckSameComm(Vec1, 1, Vec2, 2);
166 VecCheckSameSize(Vec1, 1, Vec2, 2);
167
168 PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
169 PetscCall(VecGetLocalSize(Vec1, &n));
170 if (n > 0) {
171 if (Vec1 == Vec2) {
172 PetscCall(VecGetArrayRead(Vec1, &v1));
173 v2 = v1;
174 } else {
175 PetscCall(VecGetArrayRead(Vec1, &v1));
176 PetscCall(VecGetArrayRead(Vec2, &v2));
177 }
178
179 PetscCall(PetscMalloc1(n, >));
180
181 for (i = 0; i < n; ++i) {
182 if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {
183 gt[n_gt] = low + i;
184 ++n_gt;
185 }
186 }
187
188 if (Vec1 == Vec2) {
189 PetscCall(VecRestoreArrayRead(Vec1, &v1));
190 } else {
191 PetscCall(VecRestoreArrayRead(Vec1, &v1));
192 PetscCall(VecRestoreArrayRead(Vec2, &v2));
193 }
194 }
195 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_gt, gt, PETSC_OWN_POINTER, S));
196 PetscFunctionReturn(PETSC_SUCCESS);
197 }
198
199 /*@
200 VecWhichBetween - Creates an index set containing the indices
201 where `VecLow` < `V` < `VecHigh`
202
203 Collective
204
205 Input Parameters:
206 + VecLow - lower bound
207 . V - Vector to compare
208 - VecHigh - higher bound
209
210 Output Parameter:
211 . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]
212
213 Level: advanced
214
215 Notes:
216 The vectors must have the same parallel layout
217
218 For complex numbers this only compares the real part
219
220 .seealso: `Vec`
221 @*/
VecWhichBetween(Vec VecLow,Vec V,Vec VecHigh,IS * S)222 PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
223 {
224 PetscInt i, n_vm = 0;
225 PetscInt n, low, high;
226 PetscInt *vm = NULL;
227 const PetscScalar *v1, *v2, *vmiddle;
228
229 PetscFunctionBegin;
230 PetscValidHeaderSpecific(VecLow, VEC_CLASSID, 1);
231 PetscValidHeaderSpecific(V, VEC_CLASSID, 2);
232 PetscValidHeaderSpecific(VecHigh, VEC_CLASSID, 3);
233 PetscCheckSameComm(V, 2, VecLow, 1);
234 PetscCheckSameComm(V, 2, VecHigh, 3);
235 VecCheckSameSize(V, 2, VecLow, 1);
236 VecCheckSameSize(V, 2, VecHigh, 3);
237
238 PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
239 PetscCall(VecGetLocalSize(VecLow, &n));
240 if (n > 0) {
241 PetscCall(VecGetArrayRead(VecLow, &v1));
242 if (VecLow != VecHigh) {
243 PetscCall(VecGetArrayRead(VecHigh, &v2));
244 } else {
245 v2 = v1;
246 }
247 if (V != VecLow && V != VecHigh) {
248 PetscCall(VecGetArrayRead(V, &vmiddle));
249 } else if (V == VecLow) {
250 vmiddle = v1;
251 } else {
252 vmiddle = v2;
253 }
254
255 PetscCall(PetscMalloc1(n, &vm));
256
257 for (i = 0; i < n; ++i) {
258 if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {
259 vm[n_vm] = low + i;
260 ++n_vm;
261 }
262 }
263
264 PetscCall(VecRestoreArrayRead(VecLow, &v1));
265 if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
266 if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &vmiddle));
267 }
268 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
269 PetscFunctionReturn(PETSC_SUCCESS);
270 }
271
272 /*@
273 VecWhichBetweenOrEqual - Creates an index set containing the indices
274 where `VecLow` <= `V` <= `VecHigh`
275
276 Collective
277
278 Input Parameters:
279 + VecLow - lower bound
280 . V - Vector to compare
281 - VecHigh - higher bound
282
283 Output Parameter:
284 . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
285
286 Level: advanced
287
288 .seealso: `Vec`
289 @*/
VecWhichBetweenOrEqual(Vec VecLow,Vec V,Vec VecHigh,IS * S)290 PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS *S)
291 {
292 PetscInt i, n_vm = 0;
293 PetscInt n, low, high;
294 PetscInt *vm = NULL;
295 const PetscScalar *v1, *v2, *vmiddle;
296
297 PetscFunctionBegin;
298 PetscValidHeaderSpecific(VecLow, VEC_CLASSID, 1);
299 PetscValidHeaderSpecific(V, VEC_CLASSID, 2);
300 PetscValidHeaderSpecific(VecHigh, VEC_CLASSID, 3);
301 PetscCheckSameComm(V, 2, VecLow, 1);
302 PetscCheckSameComm(V, 2, VecHigh, 3);
303 VecCheckSameSize(V, 2, VecLow, 1);
304 VecCheckSameSize(V, 2, VecHigh, 3);
305
306 PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
307 PetscCall(VecGetLocalSize(VecLow, &n));
308 if (n > 0) {
309 PetscCall(VecGetArrayRead(VecLow, &v1));
310 if (VecLow != VecHigh) {
311 PetscCall(VecGetArrayRead(VecHigh, &v2));
312 } else {
313 v2 = v1;
314 }
315 if (V != VecLow && V != VecHigh) {
316 PetscCall(VecGetArrayRead(V, &vmiddle));
317 } else if (V == VecLow) {
318 vmiddle = v1;
319 } else {
320 vmiddle = v2;
321 }
322
323 PetscCall(PetscMalloc1(n, &vm));
324
325 for (i = 0; i < n; ++i) {
326 if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {
327 vm[n_vm] = low + i;
328 ++n_vm;
329 }
330 }
331
332 PetscCall(VecRestoreArrayRead(VecLow, &v1));
333 if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
334 if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &vmiddle));
335 }
336 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
337 PetscFunctionReturn(PETSC_SUCCESS);
338 }
339
340 /*@
341 VecWhichInactive - Creates an `IS` based on a set of vectors
342
343 Collective
344
345 Input Parameters:
346 + VecLow - lower bound
347 . V - Vector to compare
348 . D - Direction to compare
349 . VecHigh - higher bound
350 - Strong - indicator for applying strongly inactive test
351
352 Output Parameter:
353 . S - The index set containing the indices i where the bound is inactive
354
355 Level: advanced
356
357 Notes:
358 Creates an index set containing the indices where one of the following holds\:
359 .vb
360 - VecLow(i) < V(i) < VecHigh(i)
361 - VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
362 - VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
363 .ve
364
365 .seealso: `Vec`
366 @*/
VecWhichInactive(Vec VecLow,Vec V,Vec D,Vec VecHigh,PetscBool Strong,IS * S)367 PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS *S)
368 {
369 PetscInt i, n_vm = 0;
370 PetscInt n, low, high;
371 PetscInt *vm = NULL;
372 const PetscScalar *v1, *v2, *v, *d;
373
374 PetscFunctionBegin;
375 if (!VecLow && !VecHigh) {
376 PetscCall(VecGetOwnershipRange(V, &low, &high));
377 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)V), high - low, low, 1, S));
378 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380 PetscValidHeaderSpecific(VecLow, VEC_CLASSID, 1);
381 PetscValidHeaderSpecific(V, VEC_CLASSID, 2);
382 PetscValidHeaderSpecific(D, VEC_CLASSID, 3);
383 PetscValidHeaderSpecific(VecHigh, VEC_CLASSID, 4);
384 PetscCheckSameComm(V, 2, VecLow, 1);
385 PetscCheckSameComm(V, 2, D, 3);
386 PetscCheckSameComm(V, 2, VecHigh, 4);
387 VecCheckSameSize(V, 2, VecLow, 1);
388 VecCheckSameSize(V, 2, D, 3);
389 VecCheckSameSize(V, 2, VecHigh, 4);
390
391 PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
392 PetscCall(VecGetLocalSize(VecLow, &n));
393 if (n > 0) {
394 PetscCall(VecGetArrayRead(VecLow, &v1));
395 if (VecLow != VecHigh) {
396 PetscCall(VecGetArrayRead(VecHigh, &v2));
397 } else {
398 v2 = v1;
399 }
400 if (V != VecLow && V != VecHigh) {
401 PetscCall(VecGetArrayRead(V, &v));
402 } else if (V == VecLow) {
403 v = v1;
404 } else {
405 v = v2;
406 }
407 if (D != VecLow && D != VecHigh && D != V) {
408 PetscCall(VecGetArrayRead(D, &d));
409 } else if (D == VecLow) {
410 d = v1;
411 } else if (D == VecHigh) {
412 d = v2;
413 } else {
414 d = v;
415 }
416
417 PetscCall(PetscMalloc1(n, &vm));
418
419 if (Strong) {
420 for (i = 0; i < n; ++i) {
421 if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
422 vm[n_vm] = low + i;
423 ++n_vm;
424 } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
425 vm[n_vm] = low + i;
426 ++n_vm;
427 } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
428 vm[n_vm] = low + i;
429 ++n_vm;
430 }
431 }
432 } else {
433 for (i = 0; i < n; ++i) {
434 if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
435 vm[n_vm] = low + i;
436 ++n_vm;
437 } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
438 vm[n_vm] = low + i;
439 ++n_vm;
440 } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
441 vm[n_vm] = low + i;
442 ++n_vm;
443 }
444 }
445 }
446
447 PetscCall(VecRestoreArrayRead(VecLow, &v1));
448 if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
449 if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &v));
450 if (D != VecLow && D != VecHigh && D != V) PetscCall(VecRestoreArrayRead(D, &d));
451 }
452 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
453 PetscFunctionReturn(PETSC_SUCCESS);
454 }
455
456 /*@
457 VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
458 vfull[is[i]] += alpha*vreduced[i]
459
460 Logically Collective
461
462 Input Parameters:
463 + vfull - the full-space vector
464 . is - the index set for the reduced space
465 . alpha - the scalar coefficient
466 - vreduced - the reduced-space vector
467
468 Output Parameter:
469 . vfull - the sum of the full-space vector and reduced-space vector
470
471 Level: advanced
472
473 Notes:
474 The index set identifies entries in the global vector.
475 Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
476
477 .seealso: `VecISCopy()`, `VecISSet()`, `VecAXPY()`
478 @*/
VecISAXPY(Vec vfull,IS is,PetscScalar alpha,Vec vreduced)479 PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
480 {
481 PetscInt nfull, nreduced;
482 PetscBool sorted = PETSC_FALSE;
483
484 PetscFunctionBegin;
485 PetscValidHeaderSpecific(vfull, VEC_CLASSID, 1);
486 PetscValidHeaderSpecific(is, IS_CLASSID, 2);
487 PetscCheckSameComm(vfull, 1, is, 2);
488 PetscValidLogicalCollectiveScalar(vfull, alpha, 3);
489 PetscValidHeaderSpecific(vreduced, VEC_CLASSID, 4);
490 PetscCall(ISGetSize(is, &nfull));
491 if (!nfull) PetscFunctionReturn(PETSC_SUCCESS);
492 PetscCall(VecGetSize(vfull, &nfull));
493 PetscCall(VecGetSize(vreduced, &nreduced));
494 if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
495 if (sorted) PetscCall(VecAXPY(vfull, alpha, vreduced));
496 else {
497 PetscScalar *y;
498 const PetscScalar *x;
499 PetscInt i, n, m, rstart, rend;
500 const PetscInt *id;
501
502 PetscCall(VecGetArray(vfull, &y));
503 PetscCall(VecGetArrayRead(vreduced, &x));
504 PetscCall(ISGetIndices(is, &id));
505 PetscCall(ISGetLocalSize(is, &n));
506 PetscCall(VecGetLocalSize(vreduced, &m));
507 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
508 PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
509 y = PetscSafePointerPlusOffset(y, -rstart);
510 for (i = 0; i < n; ++i) {
511 if (id[i] < 0) continue;
512 PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
513 y[id[i]] += alpha * x[i];
514 }
515 y = PetscSafePointerPlusOffset(y, rstart);
516 PetscCall(ISRestoreIndices(is, &id));
517 PetscCall(VecRestoreArray(vfull, &y));
518 PetscCall(VecRestoreArrayRead(vreduced, &x));
519 }
520 PetscFunctionReturn(PETSC_SUCCESS);
521 }
522
523 /*@
524 VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
525
526 Logically Collective
527
528 Input Parameters:
529 + vfull - the full-space vector
530 . is - the index set for the reduced space
531 . mode - the direction of copying, `SCATTER_FORWARD` or `SCATTER_REVERSE`
532 - vreduced - the reduced-space vector
533
534 Output Parameter:
535 . vfull - the sum of the full-space vector and reduced-space vector
536
537 Level: advanced
538
539 Notes:
540 The index set identifies entries in the global vector.
541 Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
542 .vb
543 mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
544 mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
545 .ve
546
547 .seealso: `VecISSet()`, `VecISAXPY()`, `VecCopy()`
548 @*/
VecISCopy(Vec vfull,IS is,ScatterMode mode,Vec vreduced)549 PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
550 {
551 PetscInt nfull, nreduced;
552 PetscBool sorted = PETSC_FALSE;
553
554 PetscFunctionBegin;
555 PetscValidHeaderSpecific(vfull, VEC_CLASSID, 1);
556 PetscValidHeaderSpecific(is, IS_CLASSID, 2);
557 PetscCheckSameComm(vfull, 1, is, 2);
558 PetscValidLogicalCollectiveEnum(vfull, mode, 3);
559 PetscValidHeaderSpecific(vreduced, VEC_CLASSID, 4);
560 PetscCall(ISGetSize(is, &nfull));
561 if (!nfull) PetscFunctionReturn(PETSC_SUCCESS);
562 PetscCall(VecGetSize(vfull, &nfull));
563 PetscCall(VecGetSize(vreduced, &nreduced));
564 if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
565 if (sorted) {
566 if (mode == SCATTER_FORWARD) {
567 PetscCall(VecCopy(vreduced, vfull));
568 } else {
569 PetscCall(VecCopy(vfull, vreduced));
570 }
571 } else {
572 const PetscInt *id;
573 PetscInt i, n, m, rstart, rend;
574
575 PetscCall(ISGetIndices(is, &id));
576 PetscCall(ISGetLocalSize(is, &n));
577 PetscCall(VecGetLocalSize(vreduced, &m));
578 PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
579 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
580 if (mode == SCATTER_FORWARD) {
581 PetscScalar *y;
582 const PetscScalar *x;
583
584 PetscCall(VecGetArray(vfull, &y));
585 PetscCall(VecGetArrayRead(vreduced, &x));
586 y = PetscSafePointerPlusOffset(y, -rstart);
587 for (i = 0; i < n; ++i) {
588 if (id[i] < 0) continue;
589 PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
590 y[id[i]] = x[i];
591 }
592 y = PetscSafePointerPlusOffset(y, rstart);
593 PetscCall(VecRestoreArrayRead(vreduced, &x));
594 PetscCall(VecRestoreArray(vfull, &y));
595 } else if (mode == SCATTER_REVERSE) {
596 PetscScalar *x;
597 const PetscScalar *y;
598
599 PetscCall(VecGetArrayRead(vfull, &y));
600 PetscCall(VecGetArray(vreduced, &x));
601 for (i = 0; i < n; ++i) {
602 if (id[i] < 0) continue;
603 PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
604 x[i] = y[id[i] - rstart];
605 }
606 PetscCall(VecRestoreArray(vreduced, &x));
607 PetscCall(VecRestoreArrayRead(vfull, &y));
608 } else SETERRQ(PetscObjectComm((PetscObject)vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
609 PetscCall(ISRestoreIndices(is, &id));
610 }
611 PetscFunctionReturn(PETSC_SUCCESS);
612 }
613
614 /*@
615 ISComplementVec - Creates the complement of the index set relative to a layout defined by a `Vec`
616
617 Collective
618
619 Input Parameters:
620 + S - a PETSc `IS`
621 - V - the reference vector space
622
623 Output Parameter:
624 . T - the complement of S
625
626 Level: advanced
627
628 .seealso: `IS`, `Vec`, `ISCreateGeneral()`
629 @*/
ISComplementVec(IS S,Vec V,IS * T)630 PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
631 {
632 PetscInt start, end;
633
634 PetscFunctionBegin;
635 PetscCall(VecGetOwnershipRange(V, &start, &end));
636 PetscCall(ISComplement(S, start, end, T));
637 PetscFunctionReturn(PETSC_SUCCESS);
638 }
639
640 /*@
641 VecISSet - Sets the elements of a vector, specified by an index set, to a constant
642
643 Logically Collective
644
645 Input Parameters:
646 + V - the vector
647 . S - index set for the locations in the vector
648 - c - the constant
649
650 Level: advanced
651
652 Notes:
653 The index set identifies entries in the global vector.
654 Negative indices are skipped; indices outside the ownership range of V will raise an error.
655
656 .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISShift()`, `VecSet()`
657 @*/
VecISSet(Vec V,IS S,PetscScalar c)658 PetscErrorCode VecISSet(Vec V, IS S, PetscScalar c)
659 {
660 PetscInt nloc, low, high, i;
661 const PetscInt *s;
662 PetscScalar *v;
663
664 PetscFunctionBegin;
665 PetscValidHeaderSpecific(V, VEC_CLASSID, 1);
666 PetscValidHeaderSpecific(S, IS_CLASSID, 2);
667 PetscCheckSameComm(V, 1, S, 2);
668 PetscCall(ISGetSize(S, &nloc));
669 if (nloc) {
670 PetscCall(VecGetOwnershipRange(V, &low, &high));
671 PetscCall(ISGetLocalSize(S, &nloc));
672 PetscCall(ISGetIndices(S, &s));
673 PetscCall(VecGetArray(V, &v));
674 for (i = 0; i < nloc; ++i) {
675 if (s[i] < 0) continue;
676 PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
677 v[s[i] - low] = c;
678 }
679 PetscCall(ISRestoreIndices(S, &s));
680 PetscCall(VecRestoreArray(V, &v));
681 }
682 PetscFunctionReturn(PETSC_SUCCESS);
683 }
684
685 /*@
686 VecISShift - Shifts the elements of a vector, specified by an index set, by a constant
687
688 Logically Collective
689
690 Input Parameters:
691 + V - the vector
692 . S - index set for the locations in the vector
693 - c - the constant
694
695 Level: advanced
696
697 Notes:
698 The index set identifies entries in the global vector.
699 Negative indices are skipped; indices outside the ownership range of V will raise an error.
700
701 .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISSet()`, `VecShift()`
702 @*/
VecISShift(Vec V,IS S,PetscScalar c)703 PetscErrorCode VecISShift(Vec V, IS S, PetscScalar c)
704 {
705 PetscInt nloc, low, high, i;
706 const PetscInt *s;
707 PetscScalar *v;
708
709 PetscFunctionBegin;
710 PetscValidHeaderSpecific(V, VEC_CLASSID, 1);
711 PetscValidHeaderSpecific(S, IS_CLASSID, 2);
712 PetscCheckSameComm(V, 1, S, 2);
713 PetscCall(ISGetSize(S, &nloc));
714 if (nloc) {
715 PetscCall(VecGetOwnershipRange(V, &low, &high));
716 PetscCall(ISGetLocalSize(S, &nloc));
717 PetscCall(ISGetIndices(S, &s));
718 PetscCall(VecGetArray(V, &v));
719 for (i = 0; i < nloc; ++i) {
720 if (s[i] < 0) continue;
721 PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
722 v[s[i] - low] += c;
723 }
724 PetscCall(ISRestoreIndices(S, &s));
725 PetscCall(VecRestoreArray(V, &v));
726 }
727 PetscFunctionReturn(PETSC_SUCCESS);
728 }
729
730 /*@
731 VecBoundGradientProjection - Projects vector according to this definition.
732 If XL[i] < X[i] < XU[i], then GP[i] = G[i];
733 If X[i] <= XL[i], then GP[i] = min(G[i],0);
734 If X[i] >= XU[i], then GP[i] = max(G[i],0);
735
736 Input Parameters:
737 + G - current gradient vector
738 . X - current solution vector with XL[i] <= X[i] <= XU[i]
739 . XL - lower bounds
740 - XU - upper bounds
741
742 Output Parameter:
743 . GP - gradient projection vector
744
745 Level: advanced
746
747 Note:
748 `GP` may be the same vector as `G`
749
750 For complex numbers only the real part is used in the bounds.
751
752 .seealso: `Vec`
753 @*/
VecBoundGradientProjection(Vec G,Vec X,Vec XL,Vec XU,Vec GP)754 PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
755 {
756 PetscInt n, i;
757 const PetscScalar *xptr, *xlptr, *xuptr;
758 PetscScalar *gptr, *gpptr;
759 PetscScalar xval, gpval;
760
761 /* Project variables at the lower and upper bound */
762 PetscFunctionBegin;
763 PetscValidHeaderSpecific(G, VEC_CLASSID, 1);
764 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
765 if (XL) PetscValidHeaderSpecific(XL, VEC_CLASSID, 3);
766 if (XU) PetscValidHeaderSpecific(XU, VEC_CLASSID, 4);
767 PetscValidHeaderSpecific(GP, VEC_CLASSID, 5);
768 if (!XL && !XU) {
769 PetscCall(VecCopy(G, GP));
770 PetscFunctionReturn(PETSC_SUCCESS);
771 }
772
773 PetscCall(VecGetLocalSize(X, &n));
774
775 PetscCall(VecGetArrayRead(X, &xptr));
776 PetscCall(VecGetArrayRead(XL, &xlptr));
777 PetscCall(VecGetArrayRead(XU, &xuptr));
778 PetscCall(VecGetArrayPair(G, GP, &gptr, &gpptr));
779
780 for (i = 0; i < n; ++i) {
781 gpval = gptr[i];
782 xval = xptr[i];
783 if (PetscRealPart(gpval) > 0.0 && PetscRealPart(xval) <= PetscRealPart(xlptr[i])) {
784 gpval = 0.0;
785 } else if (PetscRealPart(gpval) < 0.0 && PetscRealPart(xval) >= PetscRealPart(xuptr[i])) {
786 gpval = 0.0;
787 }
788 gpptr[i] = gpval;
789 }
790
791 PetscCall(VecRestoreArrayRead(X, &xptr));
792 PetscCall(VecRestoreArrayRead(XL, &xlptr));
793 PetscCall(VecRestoreArrayRead(XU, &xuptr));
794 PetscCall(VecRestoreArrayPair(G, GP, &gptr, &gpptr));
795 PetscFunctionReturn(PETSC_SUCCESS);
796 }
797
798 /*@
799 VecStepMaxBounded - See below
800
801 Collective
802
803 Input Parameters:
804 + X - vector with no negative entries
805 . XL - lower bounds
806 . XU - upper bounds
807 - DX - step direction, can have negative, positive or zero entries
808
809 Output Parameter:
810 . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
811
812 Level: intermediate
813
814 .seealso: `Vec`
815 @*/
VecStepMaxBounded(Vec X,Vec DX,Vec XL,Vec XU,PetscReal * stepmax)816 PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
817 {
818 PetscInt i, nn;
819 const PetscScalar *xx, *dx, *xl, *xu;
820 PetscReal localmax = 0;
821
822 PetscFunctionBegin;
823 PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
824 PetscValidHeaderSpecific(DX, VEC_CLASSID, 2);
825 PetscValidHeaderSpecific(XL, VEC_CLASSID, 3);
826 PetscValidHeaderSpecific(XU, VEC_CLASSID, 4);
827
828 PetscCall(VecGetArrayRead(X, &xx));
829 PetscCall(VecGetArrayRead(XL, &xl));
830 PetscCall(VecGetArrayRead(XU, &xu));
831 PetscCall(VecGetArrayRead(DX, &dx));
832 PetscCall(VecGetLocalSize(X, &nn));
833 for (i = 0; i < nn; i++) {
834 if (PetscRealPart(dx[i]) > 0) {
835 localmax = PetscMax(localmax, PetscRealPart((xu[i] - xx[i]) / dx[i]));
836 } else if (PetscRealPart(dx[i]) < 0) {
837 localmax = PetscMax(localmax, PetscRealPart((xl[i] - xx[i]) / dx[i]));
838 }
839 }
840 PetscCall(VecRestoreArrayRead(X, &xx));
841 PetscCall(VecRestoreArrayRead(XL, &xl));
842 PetscCall(VecRestoreArrayRead(XU, &xu));
843 PetscCall(VecRestoreArrayRead(DX, &dx));
844 PetscCallMPI(MPIU_Allreduce(&localmax, stepmax, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)X)));
845 PetscFunctionReturn(PETSC_SUCCESS);
846 }
847
848 /*@
849 VecStepBoundInfo - See below
850
851 Collective
852
853 Input Parameters:
854 + X - vector with no negative entries
855 . XL - lower bounds
856 . XU - upper bounds
857 - DX - step direction, can have negative, positive or zero entries
858
859 Output Parameters:
860 + boundmin - (may be `NULL` this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
861 . wolfemin - (may be `NULL` this it is not computed)
862 - boundmax - (may be `NULL` this it is not computed) minimum value so that X[i] + boundmax*DX[i] <= XL[i] or XU[i] <= X[i] + boundmax*DX[i]
863
864 Level: advanced
865
866 Note:
867 For complex numbers only compares the real part
868
869 .seealso: `Vec`
870 @*/
VecStepBoundInfo(Vec X,Vec DX,Vec XL,Vec XU,PetscReal * boundmin,PetscReal * wolfemin,PetscReal * boundmax)871 PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
872 {
873 PetscInt n, i;
874 const PetscScalar *x, *xl, *xu, *dx;
875 PetscReal t;
876 PetscReal localmin = PETSC_INFINITY, localwolfemin = PETSC_INFINITY, localmax = -1;
877 MPI_Comm comm;
878
879 PetscFunctionBegin;
880 PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
881 PetscValidHeaderSpecific(XL, VEC_CLASSID, 3);
882 PetscValidHeaderSpecific(XU, VEC_CLASSID, 4);
883 PetscValidHeaderSpecific(DX, VEC_CLASSID, 2);
884
885 PetscCall(VecGetArrayRead(X, &x));
886 PetscCall(VecGetArrayRead(XL, &xl));
887 PetscCall(VecGetArrayRead(XU, &xu));
888 PetscCall(VecGetArrayRead(DX, &dx));
889 PetscCall(VecGetLocalSize(X, &n));
890 for (i = 0; i < n; ++i) {
891 if (PetscRealPart(dx[i]) > 0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
892 t = PetscRealPart((xu[i] - x[i]) / dx[i]);
893 localmin = PetscMin(t, localmin);
894 if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
895 localmax = PetscMax(t, localmax);
896 } else if (PetscRealPart(dx[i]) < 0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
897 t = PetscRealPart((xl[i] - x[i]) / dx[i]);
898 localmin = PetscMin(t, localmin);
899 if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
900 localmax = PetscMax(t, localmax);
901 }
902 }
903
904 PetscCall(VecRestoreArrayRead(X, &x));
905 PetscCall(VecRestoreArrayRead(XL, &xl));
906 PetscCall(VecRestoreArrayRead(XU, &xu));
907 PetscCall(VecRestoreArrayRead(DX, &dx));
908 PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
909
910 if (boundmin) {
911 PetscCallMPI(MPIU_Allreduce(&localmin, boundmin, 1, MPIU_REAL, MPIU_MIN, comm));
912 PetscCall(PetscInfo(X, "Step Bound Info: Closest Bound: %20.19e\n", (double)*boundmin));
913 }
914 if (wolfemin) {
915 PetscCallMPI(MPIU_Allreduce(&localwolfemin, wolfemin, 1, MPIU_REAL, MPIU_MIN, comm));
916 PetscCall(PetscInfo(X, "Step Bound Info: Wolfe: %20.19e\n", (double)*wolfemin));
917 }
918 if (boundmax) {
919 PetscCallMPI(MPIU_Allreduce(&localmax, boundmax, 1, MPIU_REAL, MPIU_MAX, comm));
920 if (*boundmax < 0) *boundmax = PETSC_INFINITY;
921 PetscCall(PetscInfo(X, "Step Bound Info: Max: %20.19e\n", (double)*boundmax));
922 }
923 PetscFunctionReturn(PETSC_SUCCESS);
924 }
925
926 /*@
927 VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
928
929 Collective
930
931 Input Parameters:
932 + X - vector with no negative entries
933 - DX - a step direction, can have negative, positive or zero entries
934
935 Output Parameter:
936 . step - largest value such that x[i] + step*DX[i] >= 0 for all i
937
938 Level: advanced
939
940 Note:
941 For complex numbers only compares the real part
942
943 .seealso: `Vec`
944 @*/
VecStepMax(Vec X,Vec DX,PetscReal * step)945 PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
946 {
947 PetscInt i, nn;
948 PetscReal stepmax = PETSC_INFINITY;
949 const PetscScalar *xx, *dx;
950
951 PetscFunctionBegin;
952 PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
953 PetscValidHeaderSpecific(DX, VEC_CLASSID, 2);
954
955 PetscCall(VecGetLocalSize(X, &nn));
956 PetscCall(VecGetArrayRead(X, &xx));
957 PetscCall(VecGetArrayRead(DX, &dx));
958 for (i = 0; i < nn; ++i) {
959 PetscCheck(PetscRealPart(xx[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector must be positive");
960 if (PetscRealPart(dx[i]) < 0) stepmax = PetscMin(stepmax, PetscRealPart(-xx[i] / dx[i]));
961 }
962 PetscCall(VecRestoreArrayRead(X, &xx));
963 PetscCall(VecRestoreArrayRead(DX, &dx));
964 PetscCallMPI(MPIU_Allreduce(&stepmax, step, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)X)));
965 PetscFunctionReturn(PETSC_SUCCESS);
966 }
967
968 /*@
969 VecPow - Replaces each component of a vector by $ x_i^p $
970
971 Logically Collective
972
973 Input Parameters:
974 + v - the vector
975 - p - the exponent to use on each element
976
977 Level: intermediate
978
979 Note:
980 This handles negative values, in infinity, and NaN in the expected IEEE floating pointing manner. For example, the square root of a negative real number is NaN
981 and 1/0.0 is infinity.
982
983 .seealso: `Vec`
984 @*/
VecPow(Vec v,PetscScalar p)985 PetscErrorCode VecPow(Vec v, PetscScalar p)
986 {
987 PetscInt n, i;
988 PetscScalar *v1;
989
990 PetscFunctionBegin;
991 PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
992
993 PetscCall(VecGetArray(v, &v1));
994 PetscCall(VecGetLocalSize(v, &n));
995
996 if (1.0 == p) {
997 } else if (-1.0 == p) {
998 for (i = 0; i < n; ++i) v1[i] = 1.0 / v1[i];
999 } else if (0.0 == p) {
1000 for (i = 0; i < n; ++i) {
1001 /* Not-a-number left alone
1002 Infinity set to one */
1003 if (!PetscIsNanScalar(v1[i])) v1[i] = 1.0;
1004 }
1005 } else if (0.5 == p) {
1006 for (i = 0; i < n; ++i) v1[i] = PetscSqrtScalar(v1[i]);
1007 } else if (-0.5 == p) {
1008 for (i = 0; i < n; ++i) v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
1009 } else if (2.0 == p) {
1010 for (i = 0; i < n; ++i) v1[i] *= v1[i];
1011 } else if (-2.0 == p) {
1012 for (i = 0; i < n; ++i) v1[i] = 1.0 / (v1[i] * v1[i]);
1013 } else {
1014 for (i = 0; i < n; ++i) v1[i] = PetscPowScalar(v1[i], p);
1015 }
1016 PetscCall(VecRestoreArray(v, &v1));
1017 PetscFunctionReturn(PETSC_SUCCESS);
1018 }
1019
1020 /*@
1021 VecMedian - Computes the componentwise median of three vectors
1022 and stores the result in this vector. Used primarily for projecting
1023 a vector within upper and lower bounds.
1024
1025 Logically Collective
1026
1027 Input Parameters:
1028 + Vec1 - The first vector
1029 . Vec2 - The second vector
1030 - Vec3 - The third vector
1031
1032 Output Parameter:
1033 . VMedian - The median vector (this can be any one of the input vectors)
1034
1035 Level: advanced
1036
1037 .seealso: `Vec`
1038 @*/
VecMedian(Vec Vec1,Vec Vec2,Vec Vec3,Vec VMedian)1039 PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
1040 {
1041 PetscInt i, n, low1, high1;
1042 const PetscScalar *v1, *v2, *v3;
1043 PetscScalar *vmed;
1044
1045 PetscFunctionBegin;
1046 if (Vec1) PetscValidHeaderSpecific(Vec1, VEC_CLASSID, 1);
1047 PetscValidHeaderSpecific(Vec2, VEC_CLASSID, 2);
1048 if (Vec3) PetscValidHeaderSpecific(Vec3, VEC_CLASSID, 3);
1049 PetscValidHeaderSpecific(VMedian, VEC_CLASSID, 4);
1050
1051 if (!Vec1 && !Vec3) {
1052 PetscCall(VecCopy(Vec2, VMedian));
1053 PetscFunctionReturn(PETSC_SUCCESS);
1054 }
1055 if (Vec1 == Vec2 || Vec1 == Vec3) {
1056 PetscCall(VecCopy(Vec1, VMedian));
1057 PetscFunctionReturn(PETSC_SUCCESS);
1058 }
1059 if (Vec2 == Vec3) {
1060 PetscCall(VecCopy(Vec2, VMedian));
1061 PetscFunctionReturn(PETSC_SUCCESS);
1062 }
1063
1064 /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1065 PetscValidType(Vec1, 1);
1066 PetscValidType(Vec2, 2);
1067 PetscValidType(Vec3, 3);
1068 PetscValidType(VMedian, 4);
1069 PetscCheckSameType(Vec1, 1, Vec2, 2);
1070 PetscCheckSameType(Vec1, 1, Vec3, 3);
1071 PetscCheckSameType(Vec1, 1, VMedian, 4);
1072 PetscCheckSameComm(Vec1, 1, Vec2, 2);
1073 PetscCheckSameComm(Vec1, 1, Vec3, 3);
1074 PetscCheckSameComm(Vec1, 1, VMedian, 4);
1075 VecCheckSameSize(Vec1, 1, Vec2, 2);
1076 VecCheckSameSize(Vec1, 1, Vec3, 3);
1077 VecCheckSameSize(Vec1, 1, VMedian, 4);
1078
1079 PetscCall(VecGetOwnershipRange(Vec1, &low1, &high1));
1080 PetscCall(VecGetLocalSize(Vec1, &n));
1081 if (n > 0) {
1082 PetscCall(VecGetArray(VMedian, &vmed));
1083 if (Vec1 != VMedian) {
1084 PetscCall(VecGetArrayRead(Vec1, &v1));
1085 } else {
1086 v1 = vmed;
1087 }
1088 if (Vec2 != VMedian) {
1089 PetscCall(VecGetArrayRead(Vec2, &v2));
1090 } else {
1091 v2 = vmed;
1092 }
1093 if (Vec3 != VMedian) {
1094 PetscCall(VecGetArrayRead(Vec3, &v3));
1095 } else {
1096 v3 = vmed;
1097 }
1098
1099 for (i = 0; i < n; ++i) vmed[i] = PetscMax(PetscMax(PetscMin(PetscRealPart(v1[i]), PetscRealPart(v2[i])), PetscMin(PetscRealPart(v1[i]), PetscRealPart(v3[i]))), PetscMin(PetscRealPart(v2[i]), PetscRealPart(v3[i])));
1100
1101 PetscCall(VecRestoreArray(VMedian, &vmed));
1102 if (VMedian != Vec1) PetscCall(VecRestoreArrayRead(Vec1, &v1));
1103 if (VMedian != Vec2) PetscCall(VecRestoreArrayRead(Vec2, &v2));
1104 if (VMedian != Vec3) PetscCall(VecRestoreArrayRead(Vec3, &v3));
1105 }
1106 PetscFunctionReturn(PETSC_SUCCESS);
1107 }
1108