Actual source code: vinv.c
1: /*
2: Some useful vector utility functions.
3: */
4: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
6: /*@
7: VecStrideSet - Sets a subvector of a vector defined
8: by a starting point and a stride with a given value
10: Logically Collective
12: Input Parameters:
13: + v - the vector
14: . start - starting point of the subvector (defined by a stride)
15: - s - value to set for each entry in that subvector
17: Level: advanced
19: Notes:
20: One must call `VecSetBlockSize()` before this routine to set the stride
21: information, or use a vector created from a multicomponent `DMDA`.
23: This will only work if the desire subvector is a stride subvector
25: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideScale()`
26: @*/
27: PetscErrorCode VecStrideSet(Vec v, PetscInt start, PetscScalar s)
28: {
29: PetscInt i, n, bs;
30: PetscScalar *x;
32: PetscFunctionBegin;
34: PetscCall(VecGetLocalSize(v, &n));
35: PetscCall(VecGetBlockSize(v, &bs));
36: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
37: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
38: PetscCall(VecGetArray(v, &x));
39: for (i = start; i < n; i += bs) x[i] = s;
40: PetscCall(VecRestoreArray(v, &x));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: /*@
45: VecStrideScale - Scales a subvector of a vector defined
46: by a starting point and a stride.
48: Logically Collective
50: Input Parameters:
51: + v - the vector
52: . start - starting point of the subvector (defined by a stride)
53: - scale - value to multiply each subvector entry by
55: Level: advanced
57: Notes:
58: One must call `VecSetBlockSize()` before this routine to set the stride
59: information, or use a vector created from a multicomponent `DMDA`.
61: This will only work if the desire subvector is a stride subvector
63: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
64: @*/
65: PetscErrorCode VecStrideScale(Vec v, PetscInt start, PetscScalar scale)
66: {
67: PetscInt i, n, bs;
68: PetscScalar *x;
70: PetscFunctionBegin;
72: PetscCall(VecGetLocalSize(v, &n));
73: PetscCall(VecGetBlockSize(v, &bs));
74: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
75: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
76: PetscCall(VecGetArray(v, &x));
77: for (i = start; i < n; i += bs) x[i] *= scale;
78: PetscCall(VecRestoreArray(v, &x));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: /*@
83: VecStrideNorm - Computes the norm of subvector of a vector defined
84: by a starting point and a stride.
86: Collective
88: Input Parameters:
89: + v - the vector
90: . start - starting point of the subvector (defined by a stride)
91: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`
93: Output Parameter:
94: . nrm - the norm
96: Level: advanced
98: Notes:
99: One must call `VecSetBlockSize()` before this routine to set the stride
100: information, or use a vector created from a multicomponent `DMDA`.
102: If x is the array representing the vector x then this computes the norm
103: of the array (x[start],x[start+stride],x[start+2*stride], ....)
105: This is useful for computing, say the norm of the pressure variable when
106: the pressure is stored (interlaced) with other variables, say density etc.
108: This will only work if the desire subvector is a stride subvector
110: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
111: @*/
112: PetscErrorCode VecStrideNorm(Vec v, PetscInt start, NormType ntype, PetscReal *nrm)
113: {
114: PetscInt i, n, bs;
115: const PetscScalar *x;
116: PetscReal tnorm;
118: PetscFunctionBegin;
121: PetscAssertPointer(nrm, 4);
122: PetscCall(VecGetLocalSize(v, &n));
123: PetscCall(VecGetBlockSize(v, &bs));
124: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
125: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
126: PetscCall(VecGetArrayRead(v, &x));
127: if (ntype == NORM_2) {
128: PetscScalar sum = 0.0;
129: for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
130: tnorm = PetscRealPart(sum);
131: PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
132: *nrm = PetscSqrtReal(*nrm);
133: } else if (ntype == NORM_1) {
134: tnorm = 0.0;
135: for (i = start; i < n; i += bs) tnorm += PetscAbsScalar(x[i]);
136: PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
137: } else if (ntype == NORM_INFINITY) {
138: tnorm = 0.0;
139: for (i = start; i < n; i += bs) {
140: if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
141: }
142: PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
143: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
144: PetscCall(VecRestoreArrayRead(v, &x));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*@
149: VecStrideMax - Computes the maximum of subvector of a vector defined
150: by a starting point and a stride and optionally its location.
152: Collective
154: Input Parameters:
155: + v - the vector
156: - start - starting point of the subvector (defined by a stride)
158: Output Parameters:
159: + idex - the location where the maximum occurred (pass `NULL` if not required)
160: - nrm - the maximum value in the subvector
162: Level: advanced
164: Notes:
165: One must call `VecSetBlockSize()` before this routine to set the stride
166: information, or use a vector created from a multicomponent `DMDA`.
168: If xa is the array representing the vector x, then this computes the max
169: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
171: This is useful for computing, say the maximum of the pressure variable when
172: the pressure is stored (interlaced) with other variables, e.g., density, etc.
173: This will only work if the desire subvector is a stride subvector.
175: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
176: @*/
177: PetscErrorCode VecStrideMax(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
178: {
179: PetscInt i, n, bs, id = -1;
180: const PetscScalar *x;
181: PetscReal max = PETSC_MIN_REAL;
183: PetscFunctionBegin;
185: PetscAssertPointer(nrm, 4);
186: PetscCall(VecGetLocalSize(v, &n));
187: PetscCall(VecGetBlockSize(v, &bs));
188: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
189: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
190: PetscCall(VecGetArrayRead(v, &x));
191: for (i = start; i < n; i += bs) {
192: if (PetscRealPart(x[i]) > max) {
193: max = PetscRealPart(x[i]);
194: id = i;
195: }
196: }
197: PetscCall(VecRestoreArrayRead(v, &x));
198: #if defined(PETSC_HAVE_MPIUNI)
199: *nrm = max;
200: if (idex) *idex = id;
201: #else
202: if (!idex) {
203: PetscCall(MPIU_Allreduce(&max, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
204: } else {
205: struct {
206: PetscReal v;
207: PetscInt i;
208: } in, out;
209: PetscInt rstart;
211: PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
212: in.v = max;
213: in.i = rstart + id;
214: PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v)));
215: *nrm = out.v;
216: *idex = out.i;
217: }
218: #endif
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@
223: VecStrideMin - Computes the minimum of subvector of a vector defined
224: by a starting point and a stride and optionally its location.
226: Collective
228: Input Parameters:
229: + v - the vector
230: - start - starting point of the subvector (defined by a stride)
232: Output Parameters:
233: + idex - the location where the minimum occurred. (pass `NULL` if not required)
234: - nrm - the minimum value in the subvector
236: Level: advanced
238: Notes:
239: One must call `VecSetBlockSize()` before this routine to set the stride
240: information, or use a vector created from a multicomponent `DMDA`.
242: If xa is the array representing the vector x, then this computes the min
243: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
245: This is useful for computing, say the minimum of the pressure variable when
246: the pressure is stored (interlaced) with other variables, e.g., density, etc.
247: This will only work if the desire subvector is a stride subvector.
249: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
250: @*/
251: PetscErrorCode VecStrideMin(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
252: {
253: PetscInt i, n, bs, id = -1;
254: const PetscScalar *x;
255: PetscReal min = PETSC_MAX_REAL;
257: PetscFunctionBegin;
259: PetscAssertPointer(nrm, 4);
260: PetscCall(VecGetLocalSize(v, &n));
261: PetscCall(VecGetBlockSize(v, &bs));
262: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
263: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\nHave you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
264: PetscCall(VecGetArrayRead(v, &x));
265: for (i = start; i < n; i += bs) {
266: if (PetscRealPart(x[i]) < min) {
267: min = PetscRealPart(x[i]);
268: id = i;
269: }
270: }
271: PetscCall(VecRestoreArrayRead(v, &x));
272: #if defined(PETSC_HAVE_MPIUNI)
273: *nrm = min;
274: if (idex) *idex = id;
275: #else
276: if (!idex) {
277: PetscCall(MPIU_Allreduce(&min, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v)));
278: } else {
279: struct {
280: PetscReal v;
281: PetscInt i;
282: } in, out;
283: PetscInt rstart;
285: PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
286: in.v = min;
287: in.i = rstart + id;
288: PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v)));
289: *nrm = out.v;
290: *idex = out.i;
291: }
292: #endif
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@
297: VecStrideSum - Computes the sum of subvector of a vector defined
298: by a starting point and a stride.
300: Collective
302: Input Parameters:
303: + v - the vector
304: - start - starting point of the subvector (defined by a stride)
306: Output Parameter:
307: . sum - the sum
309: Level: advanced
311: Notes:
312: One must call `VecSetBlockSize()` before this routine to set the stride
313: information, or use a vector created from a multicomponent `DMDA`.
315: If x is the array representing the vector x then this computes the sum
316: of the array (x[start],x[start+stride],x[start+2*stride], ....)
318: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
319: @*/
320: PetscErrorCode VecStrideSum(Vec v, PetscInt start, PetscScalar *sum)
321: {
322: PetscInt i, n, bs;
323: const PetscScalar *x;
324: PetscScalar local_sum = 0.0;
326: PetscFunctionBegin;
328: PetscAssertPointer(sum, 3);
329: PetscCall(VecGetLocalSize(v, &n));
330: PetscCall(VecGetBlockSize(v, &bs));
331: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
332: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
333: PetscCall(VecGetArrayRead(v, &x));
334: for (i = start; i < n; i += bs) local_sum += x[i];
335: PetscCall(MPIU_Allreduce(&local_sum, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
336: PetscCall(VecRestoreArrayRead(v, &x));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@
341: VecStrideScaleAll - Scales the subvectors of a vector defined
342: by a starting point and a stride.
344: Logically Collective
346: Input Parameters:
347: + v - the vector
348: - scales - values to multiply each subvector entry by
350: Level: advanced
352: Notes:
353: One must call `VecSetBlockSize()` before this routine to set the stride
354: information, or use a vector created from a multicomponent `DMDA`.
356: The dimension of scales must be the same as the vector block size
358: .seealso: `Vec`, `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
359: @*/
360: PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
361: {
362: PetscInt i, j, n, bs;
363: PetscScalar *x;
365: PetscFunctionBegin;
367: PetscAssertPointer(scales, 2);
368: PetscCall(VecGetLocalSize(v, &n));
369: PetscCall(VecGetBlockSize(v, &bs));
370: PetscCall(VecGetArray(v, &x));
371: /* need to provide optimized code for each bs */
372: for (i = 0; i < n; i += bs) {
373: for (j = 0; j < bs; j++) x[i + j] *= scales[j];
374: }
375: PetscCall(VecRestoreArray(v, &x));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /*@
380: VecStrideNormAll - Computes the norms of subvectors of a vector defined
381: by a starting point and a stride.
383: Collective
385: Input Parameters:
386: + v - the vector
387: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`
389: Output Parameter:
390: . nrm - the norms
392: Level: advanced
394: Notes:
395: One must call `VecSetBlockSize()` before this routine to set the stride
396: information, or use a vector created from a multicomponent `DMDA`.
398: If x is the array representing the vector x then this computes the norm
399: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
401: The dimension of nrm must be the same as the vector block size
403: This will only work if the desire subvector is a stride subvector
405: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
406: @*/
407: PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
408: {
409: PetscInt i, j, n, bs;
410: const PetscScalar *x;
411: PetscReal tnorm[128];
412: MPI_Comm comm;
414: PetscFunctionBegin;
417: PetscAssertPointer(nrm, 3);
418: PetscCall(VecGetLocalSize(v, &n));
419: PetscCall(VecGetArrayRead(v, &x));
420: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
422: PetscCall(VecGetBlockSize(v, &bs));
423: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
425: if (ntype == NORM_2) {
426: PetscScalar sum[128];
427: for (j = 0; j < bs; j++) sum[j] = 0.0;
428: for (i = 0; i < n; i += bs) {
429: for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
430: }
431: for (j = 0; j < bs; j++) tnorm[j] = PetscRealPart(sum[j]);
433: PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
434: for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
435: } else if (ntype == NORM_1) {
436: for (j = 0; j < bs; j++) tnorm[j] = 0.0;
438: for (i = 0; i < n; i += bs) {
439: for (j = 0; j < bs; j++) tnorm[j] += PetscAbsScalar(x[i + j]);
440: }
442: PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
443: } else if (ntype == NORM_INFINITY) {
444: PetscReal tmp;
445: for (j = 0; j < bs; j++) tnorm[j] = 0.0;
447: for (i = 0; i < n; i += bs) {
448: for (j = 0; j < bs; j++) {
449: if ((tmp = PetscAbsScalar(x[i + j])) > tnorm[j]) tnorm[j] = tmp;
450: /* check special case of tmp == NaN */
451: if (tmp != tmp) {
452: tnorm[j] = tmp;
453: break;
454: }
455: }
456: }
457: PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_MAX, comm));
458: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
459: PetscCall(VecRestoreArrayRead(v, &x));
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*@
464: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
465: by a starting point and a stride and optionally its location.
467: Collective
469: Input Parameter:
470: . v - the vector
472: Output Parameters:
473: + idex - the location where the maximum occurred (not supported, pass `NULL`,
474: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
475: - nrm - the maximum values of each subvector
477: Level: advanced
479: Notes:
480: One must call `VecSetBlockSize()` before this routine to set the stride
481: information, or use a vector created from a multicomponent `DMDA`.
483: The dimension of nrm must be the same as the vector block size
485: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
486: @*/
487: PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
488: {
489: PetscInt i, j, n, bs;
490: const PetscScalar *x;
491: PetscReal max[128], tmp;
492: MPI_Comm comm;
494: PetscFunctionBegin;
496: PetscAssertPointer(nrm, 3);
497: PetscCheck(!idex, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
498: PetscCall(VecGetLocalSize(v, &n));
499: PetscCall(VecGetArrayRead(v, &x));
500: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
502: PetscCall(VecGetBlockSize(v, &bs));
503: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
505: if (!n) {
506: for (j = 0; j < bs; j++) max[j] = PETSC_MIN_REAL;
507: } else {
508: for (j = 0; j < bs; j++) max[j] = PetscRealPart(x[j]);
510: for (i = bs; i < n; i += bs) {
511: for (j = 0; j < bs; j++) {
512: if ((tmp = PetscRealPart(x[i + j])) > max[j]) max[j] = tmp;
513: }
514: }
515: }
516: PetscCall(MPIU_Allreduce(max, nrm, bs, MPIU_REAL, MPIU_MAX, comm));
518: PetscCall(VecRestoreArrayRead(v, &x));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*@
523: VecStrideMinAll - Computes the minimum of subvector of a vector defined
524: by a starting point and a stride and optionally its location.
526: Collective
528: Input Parameter:
529: . v - the vector
531: Output Parameters:
532: + idex - the location where the minimum occurred (not supported, pass `NULL`,
533: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
534: - nrm - the minimums of each subvector
536: Level: advanced
538: Notes:
539: One must call `VecSetBlockSize()` before this routine to set the stride
540: information, or use a vector created from a multicomponent `DMDA`.
542: The dimension of `nrm` must be the same as the vector block size
544: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
545: @*/
546: PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
547: {
548: PetscInt i, n, bs, j;
549: const PetscScalar *x;
550: PetscReal min[128], tmp;
551: MPI_Comm comm;
553: PetscFunctionBegin;
555: PetscAssertPointer(nrm, 3);
556: PetscCheck(!idex, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
557: PetscCall(VecGetLocalSize(v, &n));
558: PetscCall(VecGetArrayRead(v, &x));
559: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
561: PetscCall(VecGetBlockSize(v, &bs));
562: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
564: if (!n) {
565: for (j = 0; j < bs; j++) min[j] = PETSC_MAX_REAL;
566: } else {
567: for (j = 0; j < bs; j++) min[j] = PetscRealPart(x[j]);
569: for (i = bs; i < n; i += bs) {
570: for (j = 0; j < bs; j++) {
571: if ((tmp = PetscRealPart(x[i + j])) < min[j]) min[j] = tmp;
572: }
573: }
574: }
575: PetscCall(MPIU_Allreduce(min, nrm, bs, MPIU_REAL, MPIU_MIN, comm));
577: PetscCall(VecRestoreArrayRead(v, &x));
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*@
582: VecStrideSumAll - Computes the sums of subvectors of a vector defined by a stride.
584: Collective
586: Input Parameter:
587: . v - the vector
589: Output Parameter:
590: . sums - the sums
592: Level: advanced
594: Notes:
595: One must call `VecSetBlockSize()` before this routine to set the stride
596: information, or use a vector created from a multicomponent `DMDA`.
598: If x is the array representing the vector x then this computes the sum
599: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
601: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
602: @*/
603: PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
604: {
605: PetscInt i, j, n, bs;
606: const PetscScalar *x;
607: PetscScalar local_sums[128];
608: MPI_Comm comm;
610: PetscFunctionBegin;
612: PetscAssertPointer(sums, 2);
613: PetscCall(VecGetLocalSize(v, &n));
614: PetscCall(VecGetArrayRead(v, &x));
615: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
617: PetscCall(VecGetBlockSize(v, &bs));
618: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
620: for (j = 0; j < bs; j++) local_sums[j] = 0.0;
621: for (i = 0; i < n; i += bs) {
622: for (j = 0; j < bs; j++) local_sums[j] += x[i + j];
623: }
624: PetscCall(MPIU_Allreduce(local_sums, sums, bs, MPIU_SCALAR, MPIU_SUM, comm));
626: PetscCall(VecRestoreArrayRead(v, &x));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: /*----------------------------------------------------------------------------------------------*/
631: /*@
632: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
633: separate vectors.
635: Collective
637: Input Parameters:
638: + v - the vector
639: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
641: Output Parameter:
642: . s - the location where the subvectors are stored
644: Level: advanced
646: Notes:
647: One must call `VecSetBlockSize()` before this routine to set the stride
648: information, or use a vector created from a multicomponent `DMDA`.
650: If x is the array representing the vector x then this gathers
651: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
652: for start=0,1,2,...bs-1
654: The parallel layout of the vector and the subvector must be the same;
655: i.e., nlocal of v = stride*(nlocal of s)
657: Not optimized; could be easily
659: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
660: `VecStrideScatterAll()`
661: @*/
662: PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
663: {
664: PetscInt i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
665: PetscScalar **y;
666: const PetscScalar *x;
668: PetscFunctionBegin;
670: PetscAssertPointer(s, 2);
672: PetscCall(VecGetLocalSize(v, &n));
673: PetscCall(VecGetLocalSize(s[0], &n2));
674: PetscCall(VecGetArrayRead(v, &x));
675: PetscCall(VecGetBlockSize(v, &bs));
676: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");
678: PetscCall(PetscMalloc2(bs, &y, bs, &bss));
679: nv = 0;
680: nvc = 0;
681: for (i = 0; i < bs; i++) {
682: PetscCall(VecGetBlockSize(s[i], &bss[i]));
683: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
684: PetscCall(VecGetArray(s[i], &y[i]));
685: nvc += bss[i];
686: nv++;
687: PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
688: if (nvc == bs) break;
689: }
691: n = n / bs;
693: jj = 0;
694: if (addv == INSERT_VALUES) {
695: for (j = 0; j < nv; j++) {
696: for (k = 0; k < bss[j]; k++) {
697: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
698: }
699: jj += bss[j];
700: }
701: } else if (addv == ADD_VALUES) {
702: for (j = 0; j < nv; j++) {
703: for (k = 0; k < bss[j]; k++) {
704: for (i = 0; i < n; i++) y[j][i * bss[j] + k] += x[bs * i + jj + k];
705: }
706: jj += bss[j];
707: }
708: #if !defined(PETSC_USE_COMPLEX)
709: } else if (addv == MAX_VALUES) {
710: for (j = 0; j < nv; j++) {
711: for (k = 0; k < bss[j]; k++) {
712: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
713: }
714: jj += bss[j];
715: }
716: #endif
717: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
719: PetscCall(VecRestoreArrayRead(v, &x));
720: for (i = 0; i < nv; i++) PetscCall(VecRestoreArray(s[i], &y[i]));
722: PetscCall(PetscFree2(y, bss));
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*@
727: VecStrideScatterAll - Scatters all the single components from separate vectors into
728: a multi-component vector.
730: Collective
732: Input Parameters:
733: + s - the location where the subvectors are stored
734: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
736: Output Parameter:
737: . v - the multicomponent vector
739: Level: advanced
741: Notes:
742: One must call `VecSetBlockSize()` before this routine to set the stride
743: information, or use a vector created from a multicomponent `DMDA`.
745: The parallel layout of the vector and the subvector must be the same;
746: i.e., nlocal of v = stride*(nlocal of s)
748: Not optimized; could be easily
750: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
752: @*/
753: PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
754: {
755: PetscInt i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
756: PetscScalar *x;
757: PetscScalar const **y;
759: PetscFunctionBegin;
761: PetscAssertPointer(s, 1);
763: PetscCall(VecGetLocalSize(v, &n));
764: PetscCall(VecGetLocalSize(s[0], &n2));
765: PetscCall(VecGetArray(v, &x));
766: PetscCall(VecGetBlockSize(v, &bs));
767: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");
769: PetscCall(PetscMalloc2(bs, (PetscScalar ***)&y, bs, &bss));
770: nv = 0;
771: nvc = 0;
772: for (i = 0; i < bs; i++) {
773: PetscCall(VecGetBlockSize(s[i], &bss[i]));
774: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
775: PetscCall(VecGetArrayRead(s[i], &y[i]));
776: nvc += bss[i];
777: nv++;
778: PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
779: if (nvc == bs) break;
780: }
782: n = n / bs;
784: jj = 0;
785: if (addv == INSERT_VALUES) {
786: for (j = 0; j < nv; j++) {
787: for (k = 0; k < bss[j]; k++) {
788: for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
789: }
790: jj += bss[j];
791: }
792: } else if (addv == ADD_VALUES) {
793: for (j = 0; j < nv; j++) {
794: for (k = 0; k < bss[j]; k++) {
795: for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
796: }
797: jj += bss[j];
798: }
799: #if !defined(PETSC_USE_COMPLEX)
800: } else if (addv == MAX_VALUES) {
801: for (j = 0; j < nv; j++) {
802: for (k = 0; k < bss[j]; k++) {
803: for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
804: }
805: jj += bss[j];
806: }
807: #endif
808: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
810: PetscCall(VecRestoreArray(v, &x));
811: for (i = 0; i < nv; i++) PetscCall(VecRestoreArrayRead(s[i], &y[i]));
812: PetscCall(PetscFree2(*(PetscScalar ***)&y, bss));
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
816: /*@
817: VecStrideGather - Gathers a single component from a multi-component vector into
818: another vector.
820: Collective
822: Input Parameters:
823: + v - the vector
824: . start - starting point of the subvector (defined by a stride)
825: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
827: Output Parameter:
828: . s - the location where the subvector is stored
830: Level: advanced
832: Notes:
833: One must call `VecSetBlockSize()` before this routine to set the stride
834: information, or use a vector created from a multicomponent `DMDA`.
836: If x is the array representing the vector x then this gathers
837: the array (x[start],x[start+stride],x[start+2*stride], ....)
839: The parallel layout of the vector and the subvector must be the same;
840: i.e., nlocal of v = stride*(nlocal of s)
842: Not optimized; could be easily
844: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
845: `VecStrideScatterAll()`
846: @*/
847: PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
848: {
849: PetscFunctionBegin;
852: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
853: PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
854: v->map->bs);
855: PetscUseTypeMethod(v, stridegather, start, s, addv);
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: /*@
860: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
862: Collective
864: Input Parameters:
865: + s - the single-component vector
866: . start - starting point of the subvector (defined by a stride)
867: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
869: Output Parameter:
870: . v - the location where the subvector is scattered (the multi-component vector)
872: Level: advanced
874: Notes:
875: One must call `VecSetBlockSize()` on the multi-component vector before this
876: routine to set the stride information, or use a vector created from a multicomponent `DMDA`.
878: The parallel layout of the vector and the subvector must be the same;
879: i.e., nlocal of v = stride*(nlocal of s)
881: Not optimized; could be easily
883: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
884: `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
885: @*/
886: PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
887: {
888: PetscFunctionBegin;
891: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
892: PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
893: v->map->bs);
894: PetscCall((*v->ops->stridescatter)(s, start, v, addv));
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: /*@
899: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
900: another vector.
902: Collective
904: Input Parameters:
905: + v - the vector
906: . nidx - the number of indices
907: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
908: . idxs - the indices of the components 0 <= idxs[0] ...idxs[nidx-1] < bs(s), they need not be sorted, may be null if nidx == bs(s) or is `PETSC_DETERMINE`
909: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
911: Output Parameter:
912: . s - the location where the subvector is stored
914: Level: advanced
916: Notes:
917: One must call `VecSetBlockSize()` on both vectors before this routine to set the stride
918: information, or use a vector created from a multicomponent `DMDA`.
920: The parallel layout of the vector and the subvector must be the same;
922: Not optimized; could be easily
924: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
925: `VecStrideScatterAll()`
926: @*/
927: PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
928: {
929: PetscFunctionBegin;
932: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
933: PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: /*@
938: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
940: Collective
942: Input Parameters:
943: + s - the smaller-component vector
944: . nidx - the number of indices in idx
945: . idxs - the indices of the components in the smaller-component vector, 0 <= idxs[0] ...idxs[nidx-1] < bs(s) they need not be sorted, may be null if nidx == bs(s) or is `PETSC_DETERMINE`
946: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
947: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
949: Output Parameter:
950: . v - the location where the subvector is into scattered (the multi-component vector)
952: Level: advanced
954: Notes:
955: One must call `VecSetBlockSize()` on the vectors before this
956: routine to set the stride information, or use a vector created from a multicomponent `DMDA`.
958: The parallel layout of the vector and the subvector must be the same;
960: Not optimized; could be easily
962: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
963: `VecStrideScatterAll()`
964: @*/
965: PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
966: {
967: PetscFunctionBegin;
970: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
971: PetscCall((*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv));
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
976: {
977: PetscInt i, n, bs, ns;
978: const PetscScalar *x;
979: PetscScalar *y;
981: PetscFunctionBegin;
982: PetscCall(VecGetLocalSize(v, &n));
983: PetscCall(VecGetLocalSize(s, &ns));
984: PetscCall(VecGetArrayRead(v, &x));
985: PetscCall(VecGetArray(s, &y));
987: bs = v->map->bs;
988: PetscCheck(n == ns * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Subvector length * blocksize %" PetscInt_FMT " not correct for gather from original vector %" PetscInt_FMT, ns * bs, n);
989: x += start;
990: n = n / bs;
992: if (addv == INSERT_VALUES) {
993: for (i = 0; i < n; i++) y[i] = x[bs * i];
994: } else if (addv == ADD_VALUES) {
995: for (i = 0; i < n; i++) y[i] += x[bs * i];
996: #if !defined(PETSC_USE_COMPLEX)
997: } else if (addv == MAX_VALUES) {
998: for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
999: #endif
1000: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1002: PetscCall(VecRestoreArrayRead(v, &x));
1003: PetscCall(VecRestoreArray(s, &y));
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
1008: {
1009: PetscInt i, n, bs, ns;
1010: PetscScalar *x;
1011: const PetscScalar *y;
1013: PetscFunctionBegin;
1014: PetscCall(VecGetLocalSize(v, &n));
1015: PetscCall(VecGetLocalSize(s, &ns));
1016: PetscCall(VecGetArray(v, &x));
1017: PetscCall(VecGetArrayRead(s, &y));
1019: bs = v->map->bs;
1020: PetscCheck(n == ns * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Subvector length * blocksize %" PetscInt_FMT " not correct for scatter to multicomponent vector %" PetscInt_FMT, ns * bs, n);
1021: x += start;
1022: n = n / bs;
1024: if (addv == INSERT_VALUES) {
1025: for (i = 0; i < n; i++) x[bs * i] = y[i];
1026: } else if (addv == ADD_VALUES) {
1027: for (i = 0; i < n; i++) x[bs * i] += y[i];
1028: #if !defined(PETSC_USE_COMPLEX)
1029: } else if (addv == MAX_VALUES) {
1030: for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
1031: #endif
1032: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1034: PetscCall(VecRestoreArray(v, &x));
1035: PetscCall(VecRestoreArrayRead(s, &y));
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
1040: {
1041: PetscInt i, j, n, bs, bss, ns;
1042: const PetscScalar *x;
1043: PetscScalar *y;
1045: PetscFunctionBegin;
1046: PetscCall(VecGetLocalSize(v, &n));
1047: PetscCall(VecGetLocalSize(s, &ns));
1048: PetscCall(VecGetArrayRead(v, &x));
1049: PetscCall(VecGetArray(s, &y));
1051: bs = v->map->bs;
1052: bss = s->map->bs;
1053: n = n / bs;
1055: if (PetscDefined(USE_DEBUG)) {
1056: PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1057: for (j = 0; j < nidx; j++) {
1058: PetscCheck(idxv[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxv[j]);
1059: PetscCheck(idxv[j] < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is greater than or equal to vector blocksize %" PetscInt_FMT, j, idxv[j], bs);
1060: }
1061: PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not gathering into all locations");
1062: }
1064: if (addv == INSERT_VALUES) {
1065: if (!idxs) {
1066: for (i = 0; i < n; i++) {
1067: for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
1068: }
1069: } else {
1070: for (i = 0; i < n; i++) {
1071: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
1072: }
1073: }
1074: } else if (addv == ADD_VALUES) {
1075: if (!idxs) {
1076: for (i = 0; i < n; i++) {
1077: for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
1078: }
1079: } else {
1080: for (i = 0; i < n; i++) {
1081: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
1082: }
1083: }
1084: #if !defined(PETSC_USE_COMPLEX)
1085: } else if (addv == MAX_VALUES) {
1086: if (!idxs) {
1087: for (i = 0; i < n; i++) {
1088: for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1089: }
1090: } else {
1091: for (i = 0; i < n; i++) {
1092: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1093: }
1094: }
1095: #endif
1096: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1098: PetscCall(VecRestoreArrayRead(v, &x));
1099: PetscCall(VecRestoreArray(s, &y));
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
1104: {
1105: PetscInt j, i, n, bs, ns, bss;
1106: PetscScalar *x;
1107: const PetscScalar *y;
1109: PetscFunctionBegin;
1110: PetscCall(VecGetLocalSize(v, &n));
1111: PetscCall(VecGetLocalSize(s, &ns));
1112: PetscCall(VecGetArray(v, &x));
1113: PetscCall(VecGetArrayRead(s, &y));
1115: bs = v->map->bs;
1116: bss = s->map->bs;
1117: n = n / bs;
1119: if (PetscDefined(USE_DEBUG)) {
1120: PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1121: for (j = 0; j < bss; j++) {
1122: if (idxs) {
1123: PetscCheck(idxs[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxs[j]);
1124: PetscCheck(idxs[j] < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is greater than or equal to vector blocksize %" PetscInt_FMT, j, idxs[j], bs);
1125: }
1126: }
1127: PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not scattering from all locations");
1128: }
1130: if (addv == INSERT_VALUES) {
1131: if (!idxs) {
1132: for (i = 0; i < n; i++) {
1133: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
1134: }
1135: } else {
1136: for (i = 0; i < n; i++) {
1137: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
1138: }
1139: }
1140: } else if (addv == ADD_VALUES) {
1141: if (!idxs) {
1142: for (i = 0; i < n; i++) {
1143: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
1144: }
1145: } else {
1146: for (i = 0; i < n; i++) {
1147: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
1148: }
1149: }
1150: #if !defined(PETSC_USE_COMPLEX)
1151: } else if (addv == MAX_VALUES) {
1152: if (!idxs) {
1153: for (i = 0; i < n; i++) {
1154: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1155: }
1156: } else {
1157: for (i = 0; i < n; i++) {
1158: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1159: }
1160: }
1161: #endif
1162: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1164: PetscCall(VecRestoreArray(v, &x));
1165: PetscCall(VecRestoreArrayRead(s, &y));
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: static PetscErrorCode VecApplyUnary_Private(Vec v, PetscDeviceContext dctx, const char async_name[], PetscErrorCode (*unary_op)(Vec), PetscScalar (*UnaryFunc)(PetscScalar))
1170: {
1171: PetscFunctionBegin;
1173: PetscCall(VecSetErrorIfLocked(v, 1));
1174: if (dctx) {
1175: PetscErrorCode (*unary_op_async)(Vec, PetscDeviceContext);
1177: PetscCall(PetscObjectQueryFunction((PetscObject)v, async_name, &unary_op_async));
1178: if (unary_op_async) {
1179: PetscCall((*unary_op_async)(v, dctx));
1180: PetscFunctionReturn(PETSC_SUCCESS);
1181: }
1182: }
1183: if (unary_op) {
1185: PetscCall((*unary_op)(v));
1186: } else {
1187: PetscInt n;
1188: PetscScalar *x;
1191: PetscCall(VecGetLocalSize(v, &n));
1192: PetscCall(VecGetArray(v, &x));
1193: for (PetscInt i = 0; i < n; ++i) x[i] = UnaryFunc(x[i]);
1194: PetscCall(VecRestoreArray(v, &x));
1195: }
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: static PetscScalar ScalarReciprocal_Fn(PetscScalar x)
1200: {
1201: const PetscScalar zero = 0.0;
1203: return x == zero ? zero : ((PetscScalar)1.0) / x;
1204: }
1206: PetscErrorCode VecReciprocalAsync_Private(Vec v, PetscDeviceContext dctx)
1207: {
1208: PetscFunctionBegin;
1209: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Reciprocal), v->ops->reciprocal, ScalarReciprocal_Fn));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: PetscErrorCode VecReciprocal_Default(Vec v)
1214: {
1215: PetscFunctionBegin;
1216: PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarReciprocal_Fn));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: static PetscScalar ScalarExp_Fn(PetscScalar x)
1221: {
1222: return PetscExpScalar(x);
1223: }
1225: PetscErrorCode VecExpAsync_Private(Vec v, PetscDeviceContext dctx)
1226: {
1227: PetscFunctionBegin;
1229: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Exp), v->ops->exp, ScalarExp_Fn));
1230: PetscFunctionReturn(PETSC_SUCCESS);
1231: }
1233: /*@
1234: VecExp - Replaces each component of a vector by e^x_i
1236: Not Collective
1238: Input Parameter:
1239: . v - The vector
1241: Output Parameter:
1242: . v - The vector of exponents
1244: Level: beginner
1246: .seealso: `Vec`, `VecLog()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1248: @*/
1249: PetscErrorCode VecExp(Vec v)
1250: {
1251: PetscFunctionBegin;
1252: PetscCall(VecExpAsync_Private(v, NULL));
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: static PetscScalar ScalarLog_Fn(PetscScalar x)
1257: {
1258: return PetscLogScalar(x);
1259: }
1261: PetscErrorCode VecLogAsync_Private(Vec v, PetscDeviceContext dctx)
1262: {
1263: PetscFunctionBegin;
1265: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Log), v->ops->log, ScalarLog_Fn));
1266: PetscFunctionReturn(PETSC_SUCCESS);
1267: }
1269: /*@
1270: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1272: Not Collective
1274: Input Parameter:
1275: . v - The vector
1277: Output Parameter:
1278: . v - The vector of logs
1280: Level: beginner
1282: .seealso: `Vec`, `VecExp()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1284: @*/
1285: PetscErrorCode VecLog(Vec v)
1286: {
1287: PetscFunctionBegin;
1288: PetscCall(VecLogAsync_Private(v, NULL));
1289: PetscFunctionReturn(PETSC_SUCCESS);
1290: }
1292: static PetscScalar ScalarAbs_Fn(PetscScalar x)
1293: {
1294: return PetscAbsScalar(x);
1295: }
1297: PetscErrorCode VecAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1298: {
1299: PetscFunctionBegin;
1301: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Abs), v->ops->abs, ScalarAbs_Fn));
1302: PetscFunctionReturn(PETSC_SUCCESS);
1303: }
1305: /*@
1306: VecAbs - Replaces every element in a vector with its absolute value.
1308: Logically Collective
1310: Input Parameter:
1311: . v - the vector
1313: Level: intermediate
1315: .seealso: `Vec`, `VecExp()`, `VecSqrtAbs()`, `VecReciprocal()`, `VecLog()`
1316: @*/
1317: PetscErrorCode VecAbs(Vec v)
1318: {
1319: PetscFunctionBegin;
1320: PetscCall(VecAbsAsync_Private(v, NULL));
1321: PetscFunctionReturn(PETSC_SUCCESS);
1322: }
1324: static PetscScalar ScalarConjugate_Fn(PetscScalar x)
1325: {
1326: return PetscConj(x);
1327: }
1329: PetscErrorCode VecConjugateAsync_Private(Vec v, PetscDeviceContext dctx)
1330: {
1331: PetscFunctionBegin;
1333: if (PetscDefined(USE_COMPLEX)) PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Conjugate), v->ops->conjugate, ScalarConjugate_Fn));
1334: PetscFunctionReturn(PETSC_SUCCESS);
1335: }
1337: /*@
1338: VecConjugate - Conjugates a vector. That is, replace every entry in a vector with its complex conjugate
1340: Logically Collective
1342: Input Parameter:
1343: . x - the vector
1345: Level: intermediate
1347: .seealso: [](ch_vectors), `Vec`, `VecSet()`
1348: @*/
1349: PetscErrorCode VecConjugate(Vec x)
1350: {
1351: PetscFunctionBegin;
1352: PetscCall(VecConjugateAsync_Private(x, NULL));
1353: PetscFunctionReturn(PETSC_SUCCESS);
1354: }
1356: static PetscScalar ScalarSqrtAbs_Fn(PetscScalar x)
1357: {
1358: return PetscSqrtScalar(ScalarAbs_Fn(x));
1359: }
1361: PetscErrorCode VecSqrtAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1362: {
1363: PetscFunctionBegin;
1365: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(SqrtAbs), v->ops->sqrt, ScalarSqrtAbs_Fn));
1366: PetscFunctionReturn(PETSC_SUCCESS);
1367: }
1369: /*@
1370: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1372: Not Collective
1374: Input Parameter:
1375: . v - The vector
1377: Level: beginner
1379: Note:
1380: The actual function is sqrt(|x_i|)
1382: .seealso: `Vec`, `VecLog()`, `VecExp()`, `VecReciprocal()`, `VecAbs()`
1384: @*/
1385: PetscErrorCode VecSqrtAbs(Vec v)
1386: {
1387: PetscFunctionBegin;
1388: PetscCall(VecSqrtAbsAsync_Private(v, NULL));
1389: PetscFunctionReturn(PETSC_SUCCESS);
1390: }
1392: static PetscScalar ScalarImaginaryPart_Fn(PetscScalar x)
1393: {
1394: const PetscReal imag = PetscImaginaryPart(x);
1396: #if PetscDefined(USE_COMPLEX)
1397: return PetscCMPLX(imag, 0.0);
1398: #else
1399: return imag;
1400: #endif
1401: }
1403: /*@
1404: VecImaginaryPart - Replaces a complex vector with its imginary part
1406: Collective
1408: Input Parameter:
1409: . v - the vector
1411: Level: beginner
1413: .seealso: `Vec`, `VecNorm()`, `VecRealPart()`
1414: @*/
1415: PetscErrorCode VecImaginaryPart(Vec v)
1416: {
1417: PetscFunctionBegin;
1419: PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarImaginaryPart_Fn));
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1423: static PetscScalar ScalarRealPart_Fn(PetscScalar x)
1424: {
1425: const PetscReal real = PetscRealPart(x);
1427: #if PetscDefined(USE_COMPLEX)
1428: return PetscCMPLX(real, 0.0);
1429: #else
1430: return real;
1431: #endif
1432: }
1434: /*@
1435: VecRealPart - Replaces a complex vector with its real part
1437: Collective
1439: Input Parameter:
1440: . v - the vector
1442: Level: beginner
1444: .seealso: `Vec`, `VecNorm()`, `VecImaginaryPart()`
1445: @*/
1446: PetscErrorCode VecRealPart(Vec v)
1447: {
1448: PetscFunctionBegin;
1450: PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarRealPart_Fn));
1451: PetscFunctionReturn(PETSC_SUCCESS);
1452: }
1454: /*@
1455: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1457: Collective
1459: Input Parameters:
1460: + s - first vector
1461: - t - second vector
1463: Output Parameters:
1464: + dp - s'conj(t)
1465: - nm - t'conj(t)
1467: Level: advanced
1469: Note:
1470: conj(x) is the complex conjugate of x when x is complex
1472: .seealso: `Vec`, `VecDot()`, `VecNorm()`, `VecDotBegin()`, `VecNormBegin()`, `VecDotEnd()`, `VecNormEnd()`
1474: @*/
1475: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1476: {
1477: PetscScalar work[] = {0.0, 0.0};
1479: PetscFunctionBegin;
1482: PetscAssertPointer(dp, 3);
1483: PetscAssertPointer(nm, 4);
1486: PetscCheckSameTypeAndComm(s, 1, t, 2);
1487: PetscCheck(s->map->N == t->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector global lengths");
1488: PetscCheck(s->map->n == t->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector local lengths");
1490: PetscCall(PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0));
1491: if (s->ops->dotnorm2) {
1492: PetscUseTypeMethod(s, dotnorm2, t, work, work + 1);
1493: } else {
1494: const PetscScalar *sx, *tx;
1495: PetscInt n;
1497: PetscCall(VecGetLocalSize(s, &n));
1498: PetscCall(VecGetArrayRead(s, &sx));
1499: PetscCall(VecGetArrayRead(t, &tx));
1500: for (PetscInt i = 0; i < n; ++i) {
1501: const PetscScalar txconj = PetscConj(tx[i]);
1503: work[0] += sx[i] * txconj;
1504: work[1] += tx[i] * txconj;
1505: }
1506: PetscCall(VecRestoreArrayRead(t, &tx));
1507: PetscCall(VecRestoreArrayRead(s, &sx));
1508: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, work, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
1509: PetscCall(PetscLogFlops(4.0 * n));
1510: }
1511: PetscCall(PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0));
1512: *dp = work[0];
1513: *nm = PetscRealPart(work[1]);
1514: PetscFunctionReturn(PETSC_SUCCESS);
1515: }
1517: /*@
1518: VecSum - Computes the sum of all the components of a vector.
1520: Collective
1522: Input Parameter:
1523: . v - the vector
1525: Output Parameter:
1526: . sum - the result
1528: Level: beginner
1530: .seealso: `Vec`, `VecMean()`, `VecNorm()`
1531: @*/
1532: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1533: {
1534: PetscScalar tmp = 0.0;
1536: PetscFunctionBegin;
1538: PetscAssertPointer(sum, 2);
1539: if (v->ops->sum) {
1540: PetscUseTypeMethod(v, sum, &tmp);
1541: } else {
1542: const PetscScalar *x;
1543: PetscInt n;
1545: PetscCall(VecGetLocalSize(v, &n));
1546: PetscCall(VecGetArrayRead(v, &x));
1547: for (PetscInt i = 0; i < n; ++i) tmp += x[i];
1548: PetscCall(VecRestoreArrayRead(v, &x));
1549: }
1550: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &tmp, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
1551: *sum = tmp;
1552: PetscFunctionReturn(PETSC_SUCCESS);
1553: }
1555: /*@
1556: VecMean - Computes the arithmetic mean of all the components of a vector.
1558: Collective
1560: Input Parameter:
1561: . v - the vector
1563: Output Parameter:
1564: . mean - the result
1566: Level: beginner
1568: .seealso: `Vec`, `VecSum()`, `VecNorm()`
1569: @*/
1570: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1571: {
1572: PetscInt n;
1574: PetscFunctionBegin;
1576: PetscAssertPointer(mean, 2);
1577: PetscCall(VecGetSize(v, &n));
1578: PetscCall(VecSum(v, mean));
1579: *mean /= n;
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: PetscErrorCode VecShiftAsync_Private(Vec v, PetscScalar shift, PetscDeviceContext dctx)
1584: {
1585: PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext) = NULL;
1587: PetscFunctionBegin;
1590: PetscCall(VecSetErrorIfLocked(v, 1));
1591: if (shift == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
1593: if (dctx) {
1594: PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext);
1596: PetscCall(PetscObjectQueryFunction((PetscObject)v, VecAsyncFnName(Shift), &shift_async));
1597: }
1598: if (shift_async) {
1599: PetscCall((*shift_async)(v, shift, dctx));
1600: } else if (v->ops->shift) {
1601: PetscUseTypeMethod(v, shift, shift);
1602: } else {
1603: PetscInt n;
1604: PetscScalar *x;
1606: PetscCall(VecGetLocalSize(v, &n));
1607: PetscCall(VecGetArray(v, &x));
1608: for (PetscInt i = 0; i < n; ++i) x[i] += shift;
1609: PetscCall(VecRestoreArray(v, &x));
1610: PetscCall(PetscLogFlops(n));
1611: }
1612: PetscFunctionReturn(PETSC_SUCCESS);
1613: }
1615: /*@
1616: VecShift - Shifts all of the components of a vector by computing
1617: `x[i] = x[i] + shift`.
1619: Logically Collective
1621: Input Parameters:
1622: + v - the vector
1623: - shift - the shift
1625: Level: intermediate
1627: .seealso: `Vec`
1628: @*/
1629: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1630: {
1631: PetscFunctionBegin;
1632: PetscCall(VecShiftAsync_Private(v, shift, NULL));
1633: PetscFunctionReturn(PETSC_SUCCESS);
1634: }
1636: /*@
1637: VecPermute - Permutes a vector in place using the given ordering.
1639: Input Parameters:
1640: + x - The vector
1641: . row - The ordering
1642: - inv - The flag for inverting the permutation
1644: Level: beginner
1646: Note:
1647: This function does not yet support parallel Index Sets with non-local permutations
1649: .seealso: `Vec`, `MatPermute()`
1650: @*/
1651: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1652: {
1653: const PetscScalar *array;
1654: PetscScalar *newArray;
1655: const PetscInt *idx;
1656: PetscInt i, rstart, rend;
1658: PetscFunctionBegin;
1661: PetscCall(VecSetErrorIfLocked(x, 1));
1662: PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1663: PetscCall(ISGetIndices(row, &idx));
1664: PetscCall(VecGetArrayRead(x, &array));
1665: PetscCall(PetscMalloc1(x->map->n, &newArray));
1666: if (PetscDefined(USE_DEBUG)) {
1667: for (i = 0; i < x->map->n; i++) PetscCheck(!(idx[i] < rstart) && !(idx[i] >= rend), PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Permutation index %" PetscInt_FMT " is out of bounds: %" PetscInt_FMT, i, idx[i]);
1668: }
1669: if (!inv) {
1670: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i] - rstart];
1671: } else {
1672: for (i = 0; i < x->map->n; i++) newArray[idx[i] - rstart] = array[i];
1673: }
1674: PetscCall(VecRestoreArrayRead(x, &array));
1675: PetscCall(ISRestoreIndices(row, &idx));
1676: PetscCall(VecReplaceArray(x, newArray));
1677: PetscFunctionReturn(PETSC_SUCCESS);
1678: }
1680: /*@
1681: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1682: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1683: Does NOT take round-off errors into account.
1685: Collective
1687: Input Parameters:
1688: + vec1 - the first vector
1689: - vec2 - the second vector
1691: Output Parameter:
1692: . flg - `PETSC_TRUE` if the vectors are equal; `PETSC_FALSE` otherwise.
1694: Level: intermediate
1696: .seealso: `Vec`
1697: @*/
1698: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1699: {
1700: const PetscScalar *v1, *v2;
1701: PetscInt n1, n2, N1, N2;
1702: PetscBool flg1;
1704: PetscFunctionBegin;
1707: PetscAssertPointer(flg, 3);
1708: if (vec1 == vec2) *flg = PETSC_TRUE;
1709: else {
1710: PetscCall(VecGetSize(vec1, &N1));
1711: PetscCall(VecGetSize(vec2, &N2));
1712: if (N1 != N2) flg1 = PETSC_FALSE;
1713: else {
1714: PetscCall(VecGetLocalSize(vec1, &n1));
1715: PetscCall(VecGetLocalSize(vec2, &n2));
1716: if (n1 != n2) flg1 = PETSC_FALSE;
1717: else {
1718: PetscCall(VecGetArrayRead(vec1, &v1));
1719: PetscCall(VecGetArrayRead(vec2, &v2));
1720: PetscCall(PetscArraycmp(v1, v2, n1, &flg1));
1721: PetscCall(VecRestoreArrayRead(vec1, &v1));
1722: PetscCall(VecRestoreArrayRead(vec2, &v2));
1723: }
1724: }
1725: /* combine results from all processors */
1726: PetscCall(MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1)));
1727: }
1728: PetscFunctionReturn(PETSC_SUCCESS);
1729: }
1731: /*@
1732: VecUniqueEntries - Compute the number of unique entries, and those entries
1734: Collective
1736: Input Parameter:
1737: . vec - the vector
1739: Output Parameters:
1740: + n - The number of unique entries
1741: - e - The entries
1743: Level: intermediate
1745: .seealso: `Vec`
1746: @*/
1747: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1748: {
1749: const PetscScalar *v;
1750: PetscScalar *tmp, *vals;
1751: PetscMPIInt *N, *displs, l;
1752: PetscInt ng, m, i, j, p;
1753: PetscMPIInt size;
1755: PetscFunctionBegin;
1757: PetscAssertPointer(n, 2);
1758: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1759: PetscCall(VecGetLocalSize(vec, &m));
1760: PetscCall(VecGetArrayRead(vec, &v));
1761: PetscCall(PetscMalloc2(m, &tmp, size, &N));
1762: for (i = 0, j = 0, l = 0; i < m; ++i) {
1763: /* Can speed this up with sorting */
1764: for (j = 0; j < l; ++j) {
1765: if (v[i] == tmp[j]) break;
1766: }
1767: if (j == l) {
1768: tmp[j] = v[i];
1769: ++l;
1770: }
1771: }
1772: PetscCall(VecRestoreArrayRead(vec, &v));
1773: /* Gather serial results */
1774: PetscCallMPI(MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec)));
1775: for (p = 0, ng = 0; p < size; ++p) ng += N[p];
1776: PetscCall(PetscMalloc2(ng, &vals, size + 1, &displs));
1777: for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
1778: PetscCallMPI(MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec)));
1779: /* Find unique entries */
1780: #ifdef PETSC_USE_COMPLEX
1781: SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1782: #else
1783: *n = displs[size];
1784: PetscCall(PetscSortRemoveDupsReal(n, (PetscReal *)vals));
1785: if (e) {
1786: PetscAssertPointer(e, 3);
1787: PetscCall(PetscMalloc1(*n, e));
1788: for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
1789: }
1790: PetscCall(PetscFree2(vals, displs));
1791: PetscCall(PetscFree2(tmp, N));
1792: PetscFunctionReturn(PETSC_SUCCESS);
1793: #endif
1794: }