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;
35: PetscCall(VecGetLocalSize(v, &n));
36: PetscCall(VecGetBlockSize(v, &bs));
37: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
38: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
39: PetscCall(VecGetArray(v, &x));
40: for (i = start; i < n; i += bs) x[i] = s;
41: PetscCall(VecRestoreArray(v, &x));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: /*@
46: VecStrideScale - Scales a subvector of a vector defined
47: by a starting point and a stride.
49: Logically Collective
51: Input Parameters:
52: + v - the vector
53: . start - starting point of the subvector (defined by a stride)
54: - scale - value to multiply each subvector entry by
56: Level: advanced
58: Notes:
59: One must call `VecSetBlockSize()` before this routine to set the stride
60: information, or use a vector created from a multicomponent `DMDA`.
62: This will only work if the desire subvector is a stride subvector
64: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
65: @*/
66: PetscErrorCode VecStrideScale(Vec v, PetscInt start, PetscScalar scale)
67: {
68: PetscInt i, n, bs;
69: PetscScalar *x;
71: PetscFunctionBegin;
75: PetscCall(VecGetLocalSize(v, &n));
76: PetscCall(VecGetBlockSize(v, &bs));
77: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
78: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
79: PetscCall(VecGetArray(v, &x));
80: for (i = start; i < n; i += bs) x[i] *= scale;
81: PetscCall(VecRestoreArray(v, &x));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: VecStrideNorm - Computes the norm of subvector of a vector defined
87: by a starting point and a stride.
89: Collective
91: Input Parameters:
92: + v - the vector
93: . start - starting point of the subvector (defined by a stride)
94: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`
96: Output Parameter:
97: . nrm - the norm
99: Level: advanced
101: Notes:
102: One must call `VecSetBlockSize()` before this routine to set the stride
103: information, or use a vector created from a multicomponent `DMDA`.
105: If x is the array representing the vector x then this computes the norm
106: of the array (x[start],x[start+stride],x[start+2*stride], ....)
108: This is useful for computing, say the norm of the pressure variable when
109: the pressure is stored (interlaced) with other variables, say density etc.
111: This will only work if the desire subvector is a stride subvector
113: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
114: @*/
115: PetscErrorCode VecStrideNorm(Vec v, PetscInt start, NormType ntype, PetscReal *nrm)
116: {
117: PetscInt i, n, bs;
118: const PetscScalar *x;
120: PetscFunctionBegin;
124: PetscAssertPointer(nrm, 4);
125: PetscCall(VecGetLocalSize(v, &n));
126: PetscCall(VecGetBlockSize(v, &bs));
127: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
128: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
129: PetscCall(VecGetArrayRead(v, &x));
130: if (ntype == NORM_2) {
131: PetscScalar sum = 0.0;
132: for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
133: *nrm = PetscRealPart(sum);
134: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
135: *nrm = PetscSqrtReal(*nrm);
136: } else if (ntype == NORM_1) {
137: *nrm = 0.0;
138: for (i = start; i < n; i += bs) *nrm += PetscAbsScalar(x[i]);
139: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
140: } else if (ntype == NORM_INFINITY) {
141: *nrm = 0.0;
142: for (i = start; i < n; i += bs) {
143: if (PetscAbsScalar(x[i]) > *nrm) *nrm = PetscAbsScalar(x[i]);
144: }
145: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
146: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
147: PetscCall(VecRestoreArrayRead(v, &x));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: VecStrideMax - Computes the maximum of subvector of a vector defined
153: by a starting point and a stride and optionally its location.
155: Collective
157: Input Parameters:
158: + v - the vector
159: - start - starting point of the subvector (defined by a stride)
161: Output Parameters:
162: + idex - the location where the maximum occurred (pass `NULL` if not required)
163: - nrm - the maximum value in the subvector
165: Level: advanced
167: Notes:
168: One must call `VecSetBlockSize()` before this routine to set the stride
169: information, or use a vector created from a multicomponent `DMDA`.
171: If xa is the array representing the vector x, then this computes the max
172: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
174: This is useful for computing, say the maximum of the pressure variable when
175: the pressure is stored (interlaced) with other variables, e.g., density, etc.
176: This will only work if the desire subvector is a stride subvector.
178: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
179: @*/
180: PetscErrorCode VecStrideMax(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
181: {
182: PetscInt i, n, bs, id = -1;
183: const PetscScalar *x;
185: PetscFunctionBegin;
188: PetscAssertPointer(nrm, 4);
189: *nrm = PETSC_MIN_REAL;
190: PetscCall(VecGetLocalSize(v, &n));
191: PetscCall(VecGetBlockSize(v, &bs));
192: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
193: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
194: PetscCall(VecGetArrayRead(v, &x));
195: for (i = start; i < n; i += bs) {
196: if (PetscRealPart(x[i]) > *nrm) {
197: *nrm = PetscRealPart(x[i]);
198: id = i;
199: }
200: }
201: PetscCall(VecRestoreArrayRead(v, &x));
202: #if PetscDefined(HAVE_MPIUNI)
203: if (idex) *idex = id;
204: #else
205: if (!idex) {
206: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
207: } else {
208: struct {
209: PetscReal v;
210: PetscInt i;
211: } out;
212: PetscInt rstart;
214: PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
215: out.v = *nrm;
216: out.i = rstart + id;
217: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v)));
218: *nrm = out.v;
219: *idex = out.i;
220: }
221: #endif
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /*@
226: VecStrideMin - Computes the minimum of subvector of a vector defined
227: by a starting point and a stride and optionally its location.
229: Collective
231: Input Parameters:
232: + v - the vector
233: - start - starting point of the subvector (defined by a stride)
235: Output Parameters:
236: + idex - the location where the minimum occurred. (pass `NULL` if not required)
237: - nrm - the minimum value in the subvector
239: Level: advanced
241: Notes:
242: One must call `VecSetBlockSize()` before this routine to set the stride
243: information, or use a vector created from a multicomponent `DMDA`.
245: If xa is the array representing the vector x, then this computes the min
246: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
248: This is useful for computing, say the minimum of the pressure variable when
249: the pressure is stored (interlaced) with other variables, e.g., density, etc.
250: This will only work if the desire subvector is a stride subvector.
252: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
253: @*/
254: PetscErrorCode VecStrideMin(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
255: {
256: PetscInt i, n, bs, id = -1;
257: const PetscScalar *x;
259: PetscFunctionBegin;
262: PetscAssertPointer(nrm, 4);
263: *nrm = PETSC_MAX_REAL;
264: PetscCall(VecGetLocalSize(v, &n));
265: PetscCall(VecGetBlockSize(v, &bs));
266: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
267: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
268: PetscCall(VecGetArrayRead(v, &x));
269: for (i = start; i < n; i += bs) {
270: if (PetscRealPart(x[i]) < *nrm) {
271: *nrm = PetscRealPart(x[i]);
272: id = i;
273: }
274: }
275: PetscCall(VecRestoreArrayRead(v, &x));
276: #if PetscDefined(HAVE_MPIUNI)
277: if (idex) *idex = id;
278: #else
279: if (!idex) {
280: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v)));
281: } else {
282: struct {
283: PetscReal v;
284: PetscInt i;
285: } out;
286: PetscInt rstart;
288: PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
289: out.v = *nrm;
290: out.i = rstart + id;
291: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v)));
292: *nrm = out.v;
293: *idex = out.i;
294: }
295: #endif
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: /*@
300: VecStrideSum - Computes the sum of subvector of a vector defined
301: by a starting point and a stride.
303: Collective
305: Input Parameters:
306: + v - the vector
307: - start - starting point of the subvector (defined by a stride)
309: Output Parameter:
310: . sum - the sum
312: Level: advanced
314: Notes:
315: One must call `VecSetBlockSize()` before this routine to set the stride
316: information, or use a vector created from a multicomponent `DMDA`.
318: If x is the array representing the vector x then this computes the sum
319: of the array (x[start],x[start+stride],x[start+2*stride], ....)
321: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
322: @*/
323: PetscErrorCode VecStrideSum(Vec v, PetscInt start, PetscScalar *sum)
324: {
325: PetscInt i, n, bs;
326: const PetscScalar *x;
328: PetscFunctionBegin;
331: PetscAssertPointer(sum, 3);
332: PetscCall(VecGetLocalSize(v, &n));
333: PetscCall(VecGetBlockSize(v, &bs));
334: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
335: PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
336: *sum = 0.0;
337: PetscCall(VecGetArrayRead(v, &x));
338: for (i = start; i < n; i += bs) *sum += x[i];
339: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
340: PetscCall(VecRestoreArrayRead(v, &x));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: /*@
345: VecStrideScaleAll - Scales the subvectors of a vector defined
346: by a starting point and a stride.
348: Logically Collective
350: Input Parameters:
351: + v - the vector
352: - scales - values to multiply each subvector entry by
354: Level: advanced
356: Notes:
357: One must call `VecSetBlockSize()` before this routine to set the stride
358: information, or use a vector created from a multicomponent `DMDA`.
360: The dimension of scales must be the same as the vector block size
362: .seealso: `Vec`, `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
363: @*/
364: PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
365: {
366: PetscInt i, j, n, bs;
367: PetscScalar *x;
369: PetscFunctionBegin;
371: PetscAssertPointer(scales, 2);
372: PetscCall(VecGetLocalSize(v, &n));
373: PetscCall(VecGetBlockSize(v, &bs));
374: PetscCall(VecGetArray(v, &x));
375: /* need to provide optimized code for each bs */
376: for (i = 0; i < n; i += bs) {
377: for (j = 0; j < bs; j++) x[i + j] *= scales[j];
378: }
379: PetscCall(VecRestoreArray(v, &x));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: /*@
384: VecStrideNormAll - Computes the norms of subvectors of a vector defined
385: by a starting point and a stride.
387: Collective
389: Input Parameters:
390: + v - the vector
391: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`
393: Output Parameter:
394: . nrm - the norms
396: Level: advanced
398: Notes:
399: One must call `VecSetBlockSize()` before this routine to set the stride
400: information, or use a vector created from a multicomponent `DMDA`.
402: If x is the array representing the vector x then this computes the norm
403: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
405: The dimension of nrm must be the same as the vector block size
407: This will only work if the desire subvector is a stride subvector
409: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
410: @*/
411: PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
412: {
413: PetscInt i, j, n, bs;
414: const PetscScalar *x;
415: MPI_Comm comm;
416: PetscMPIInt ibs;
418: PetscFunctionBegin;
421: PetscAssertPointer(nrm, 3);
422: PetscCall(VecGetLocalSize(v, &n));
423: PetscCall(VecGetArrayRead(v, &x));
424: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
426: PetscCall(VecGetBlockSize(v, &bs));
427: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
428: PetscCall(PetscMPIIntCast(bs, &ibs));
429: if (ntype == NORM_2) {
430: PetscScalar sum[128];
431: for (j = 0; j < bs; j++) sum[j] = 0.0;
432: for (i = 0; i < n; i += bs) {
433: for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
434: }
435: for (j = 0; j < bs; j++) nrm[j] = PetscRealPart(sum[j]);
437: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, ibs, MPIU_REAL, MPIU_SUM, comm));
438: for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
439: } else if (ntype == NORM_1) {
440: for (j = 0; j < bs; j++) nrm[j] = 0.0;
442: for (i = 0; i < n; i += bs) {
443: for (j = 0; j < bs; j++) nrm[j] += PetscAbsScalar(x[i + j]);
444: }
446: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, ibs, MPIU_REAL, MPIU_SUM, comm));
447: } else if (ntype == NORM_INFINITY) {
448: PetscReal tmp;
449: for (j = 0; j < bs; j++) nrm[j] = 0.0;
451: for (i = 0; i < n; i += bs) {
452: for (j = 0; j < bs; j++) {
453: if ((tmp = PetscAbsScalar(x[i + j])) > nrm[j]) nrm[j] = tmp;
454: /* check special case of tmp == NaN */
455: if (tmp != tmp) {
456: nrm[j] = tmp;
457: break;
458: }
459: }
460: }
461: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, ibs, MPIU_REAL, MPIU_MAX, comm));
462: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
463: PetscCall(VecRestoreArrayRead(v, &x));
464: PetscFunctionReturn(PETSC_SUCCESS);
465: }
467: /*@
468: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
469: by a starting point and a stride and optionally its location.
471: Collective
473: Input Parameter:
474: . v - the vector
476: Output Parameters:
477: + idex - the location where the maximum occurred (not supported, pass `NULL`,
478: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
479: - nrm - the maximum values of each subvector
481: Level: advanced
483: Notes:
484: One must call `VecSetBlockSize()` before this routine to set the stride
485: information, or use a vector created from a multicomponent `DMDA`.
487: The dimension of nrm must be the same as the vector block size
489: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
490: @*/
491: PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
492: {
493: PetscInt i, j, n, bs;
494: const PetscScalar *x;
495: PetscReal tmp;
496: MPI_Comm comm;
497: PetscMPIInt ibs;
499: PetscFunctionBegin;
501: PetscAssertPointer(nrm, 3);
502: 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");
503: PetscCall(VecGetLocalSize(v, &n));
504: PetscCall(VecGetArrayRead(v, &x));
505: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
507: PetscCall(VecGetBlockSize(v, &bs));
508: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
509: PetscCall(PetscMPIIntCast(bs, &ibs));
511: if (!n) {
512: for (j = 0; j < bs; j++) nrm[j] = PETSC_MIN_REAL;
513: } else {
514: for (j = 0; j < bs; j++) nrm[j] = PetscRealPart(x[j]);
516: for (i = bs; i < n; i += bs) {
517: for (j = 0; j < bs; j++) {
518: if ((tmp = PetscRealPart(x[i + j])) > nrm[j]) nrm[j] = tmp;
519: }
520: }
521: }
522: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, ibs, MPIU_REAL, MPIU_MAX, comm));
524: PetscCall(VecRestoreArrayRead(v, &x));
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: /*@
529: VecStrideMinAll - Computes the minimum of subvector of a vector defined
530: by a starting point and a stride and optionally its location.
532: Collective
534: Input Parameter:
535: . v - the vector
537: Output Parameters:
538: + idex - the location where the minimum occurred (not supported, pass `NULL`,
539: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
540: - nrm - the minimums of each subvector
542: Level: advanced
544: Notes:
545: One must call `VecSetBlockSize()` before this routine to set the stride
546: information, or use a vector created from a multicomponent `DMDA`.
548: The dimension of `nrm` must be the same as the vector block size
550: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
551: @*/
552: PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
553: {
554: PetscInt i, n, bs, j;
555: const PetscScalar *x;
556: PetscReal tmp;
557: MPI_Comm comm;
558: PetscMPIInt ibs;
560: PetscFunctionBegin;
562: PetscAssertPointer(nrm, 3);
563: 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");
564: PetscCall(VecGetLocalSize(v, &n));
565: PetscCall(VecGetArrayRead(v, &x));
566: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
568: PetscCall(VecGetBlockSize(v, &bs));
569: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
570: PetscCall(PetscMPIIntCast(bs, &ibs));
572: if (!n) {
573: for (j = 0; j < bs; j++) nrm[j] = PETSC_MAX_REAL;
574: } else {
575: for (j = 0; j < bs; j++) nrm[j] = PetscRealPart(x[j]);
577: for (i = bs; i < n; i += bs) {
578: for (j = 0; j < bs; j++) {
579: if ((tmp = PetscRealPart(x[i + j])) < nrm[j]) nrm[j] = tmp;
580: }
581: }
582: }
583: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nrm, ibs, MPIU_REAL, MPIU_MIN, comm));
585: PetscCall(VecRestoreArrayRead(v, &x));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*@
590: VecStrideSumAll - Computes the sums of subvectors of a vector defined by a stride.
592: Collective
594: Input Parameter:
595: . v - the vector
597: Output Parameter:
598: . sums - the sums
600: Level: advanced
602: Notes:
603: One must call `VecSetBlockSize()` before this routine to set the stride
604: information, or use a vector created from a multicomponent `DMDA`.
606: If x is the array representing the vector x then this computes the sum
607: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
609: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
610: @*/
611: PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
612: {
613: PetscInt i, j, n, bs;
614: const PetscScalar *x;
615: MPI_Comm comm;
616: PetscMPIInt ibs;
618: PetscFunctionBegin;
620: PetscAssertPointer(sums, 2);
621: PetscCall(VecGetLocalSize(v, &n));
622: PetscCall(VecGetArrayRead(v, &x));
623: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
625: PetscCall(VecGetBlockSize(v, &bs));
626: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
627: PetscCall(PetscMPIIntCast(bs, &ibs));
629: for (j = 0; j < bs; j++) sums[j] = 0.0;
630: for (i = 0; i < n; i += bs) {
631: for (j = 0; j < bs; j++) sums[j] += x[i + j];
632: }
633: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, sums, ibs, MPIU_SCALAR, MPIU_SUM, comm));
635: PetscCall(VecRestoreArrayRead(v, &x));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: /*@
640: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
641: separate vectors.
643: Collective
645: Input Parameters:
646: + v - the vector
647: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
649: Output Parameter:
650: . s - the location where the subvectors are stored
652: Level: advanced
654: Notes:
655: One must call `VecSetBlockSize()` before this routine to set the stride
656: information, or use a vector created from a multicomponent `DMDA`.
658: If x is the array representing the vector x then this gathers
659: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
660: for start=0,1,2,...bs-1
662: The parallel layout of the vector and the subvector must be the same;
663: i.e., nlocal of v = stride*(nlocal of s)
665: Not optimized; could be easily
667: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
668: `VecStrideScatterAll()`
669: @*/
670: PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
671: {
672: PetscInt i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
673: PetscScalar **y;
674: const PetscScalar *x;
676: PetscFunctionBegin;
678: PetscAssertPointer(s, 2);
680: PetscCall(VecGetLocalSize(v, &n));
681: PetscCall(VecGetLocalSize(s[0], &n2));
682: PetscCall(VecGetArrayRead(v, &x));
683: PetscCall(VecGetBlockSize(v, &bs));
684: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");
686: PetscCall(PetscMalloc2(bs, &y, bs, &bss));
687: nv = 0;
688: nvc = 0;
689: for (i = 0; i < bs; i++) {
690: PetscCall(VecGetBlockSize(s[i], &bss[i]));
691: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
692: PetscCall(VecGetArray(s[i], &y[i]));
693: nvc += bss[i];
694: nv++;
695: PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
696: if (nvc == bs) break;
697: }
699: n = n / bs;
701: jj = 0;
702: if (addv == INSERT_VALUES) {
703: for (j = 0; j < nv; j++) {
704: for (k = 0; k < bss[j]; k++) {
705: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
706: }
707: jj += bss[j];
708: }
709: } else if (addv == ADD_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] += x[bs * i + jj + k];
713: }
714: jj += bss[j];
715: }
716: #if !PetscDefined(USE_COMPLEX)
717: } else if (addv == MAX_VALUES) {
718: for (j = 0; j < nv; j++) {
719: for (k = 0; k < bss[j]; k++) {
720: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
721: }
722: jj += bss[j];
723: }
724: #endif
725: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
727: PetscCall(VecRestoreArrayRead(v, &x));
728: for (i = 0; i < nv; i++) PetscCall(VecRestoreArray(s[i], &y[i]));
730: PetscCall(PetscFree2(y, bss));
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: /*@
735: VecStrideScatterAll - Scatters all the single components from separate vectors into
736: a multi-component vector.
738: Collective
740: Input Parameters:
741: + s - the location where the subvectors are stored
742: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
744: Output Parameter:
745: . v - the multicomponent vector
747: Level: advanced
749: Notes:
750: One must call `VecSetBlockSize()` before this routine to set the stride
751: information, or use a vector created from a multicomponent `DMDA`.
753: The parallel layout of the vector and the subvector must be the same;
754: i.e., nlocal of v = stride*(nlocal of s)
756: Not optimized; could be easily
758: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`
759: @*/
760: PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
761: {
762: PetscInt i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
763: PetscScalar *x;
764: PetscScalar const **y;
766: PetscFunctionBegin;
768: PetscAssertPointer(s, 1);
770: PetscCall(VecGetLocalSize(v, &n));
771: PetscCall(VecGetLocalSize(s[0], &n2));
772: PetscCall(VecGetArray(v, &x));
773: PetscCall(VecGetBlockSize(v, &bs));
774: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");
776: PetscCall(PetscMalloc2(bs, (PetscScalar ***)&y, bs, &bss));
777: nv = 0;
778: nvc = 0;
779: for (i = 0; i < bs; i++) {
780: PetscCall(VecGetBlockSize(s[i], &bss[i]));
781: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
782: PetscCall(VecGetArrayRead(s[i], &y[i]));
783: nvc += bss[i];
784: nv++;
785: PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
786: if (nvc == bs) break;
787: }
789: n = n / bs;
791: jj = 0;
792: if (addv == INSERT_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: } else if (addv == ADD_VALUES) {
800: for (j = 0; j < nv; j++) {
801: for (k = 0; k < bss[j]; k++) {
802: for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
803: }
804: jj += bss[j];
805: }
806: #if !PetscDefined(USE_COMPLEX)
807: } else if (addv == MAX_VALUES) {
808: for (j = 0; j < nv; j++) {
809: for (k = 0; k < bss[j]; k++) {
810: for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
811: }
812: jj += bss[j];
813: }
814: #endif
815: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
817: PetscCall(VecRestoreArray(v, &x));
818: for (i = 0; i < nv; i++) PetscCall(VecRestoreArrayRead(s[i], &y[i]));
819: PetscCall(PetscFree2(*(PetscScalar ***)&y, bss));
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
823: /*@
824: VecStrideGather - Gathers a single component from a multi-component vector into
825: another vector.
827: Collective
829: Input Parameters:
830: + v - the vector
831: . start - starting point of the subvector (defined by a stride)
832: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
834: Output Parameter:
835: . s - the location where the subvector is stored
837: Level: advanced
839: Notes:
840: One must call `VecSetBlockSize()` before this routine to set the stride
841: information, or use a vector created from a multicomponent `DMDA`.
843: If x is the array representing the vector x then this gathers
844: the array (x[start],x[start+stride],x[start+2*stride], ....)
846: The parallel layout of the vector and the subvector must be the same;
847: i.e., nlocal of v = stride*(nlocal of s)
849: Not optimized; could be easily
851: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
852: `VecStrideScatterAll()`
853: @*/
854: PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
855: {
856: PetscFunctionBegin;
860: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
861: PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
862: v->map->bs);
863: PetscUseTypeMethod(v, stridegather, start, s, addv);
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
870: Collective
872: Input Parameters:
873: + s - the single-component vector
874: . start - starting point of the subvector (defined by a stride)
875: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
877: Output Parameter:
878: . v - the location where the subvector is scattered (the multi-component vector)
880: Level: advanced
882: Notes:
883: One must call `VecSetBlockSize()` on the multi-component vector before this
884: routine to set the stride information, or use a vector created from a multicomponent `DMDA`.
886: The parallel layout of the vector and the subvector must be the same;
887: i.e., nlocal of v = stride*(nlocal of s)
889: Not optimized; could be easily
891: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
892: `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
893: @*/
894: PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
895: {
896: PetscFunctionBegin;
900: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
901: PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
902: v->map->bs);
903: PetscCall((*v->ops->stridescatter)(s, start, v, addv));
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: /*@
908: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
909: another vector.
911: Collective
913: Input Parameters:
914: + v - the vector
915: . nidx - the number of indices
916: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
917: . 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`
918: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
920: Output Parameter:
921: . s - the location where the subvector is stored
923: Level: advanced
925: Notes:
926: One must call `VecSetBlockSize()` on both vectors before this routine to set the stride
927: information, or use a vector created from a multicomponent `DMDA`.
929: The parallel layout of the vector and the subvector must be the same;
931: Not optimized; could be easily
933: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
934: `VecStrideScatterAll()`
935: @*/
936: PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
937: {
938: PetscFunctionBegin;
941: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
942: PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: /*@
947: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
949: Collective
951: Input Parameters:
952: + s - the smaller-component vector
953: . nidx - the number of indices in idx
954: . 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`
955: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
956: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
958: Output Parameter:
959: . v - the location where the subvector is into scattered (the multi-component vector)
961: Level: advanced
963: Notes:
964: One must call `VecSetBlockSize()` on the vectors before this
965: routine to set the stride information, or use a vector created from a multicomponent `DMDA`.
967: The parallel layout of the vector and the subvector must be the same;
969: Not optimized; could be easily
971: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
972: `VecStrideScatterAll()`
973: @*/
974: PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
975: {
976: PetscFunctionBegin;
979: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
980: PetscCall((*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv));
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
985: {
986: PetscInt i, n, bs, ns;
987: const PetscScalar *x;
988: PetscScalar *y;
990: PetscFunctionBegin;
991: PetscCall(VecGetLocalSize(v, &n));
992: PetscCall(VecGetLocalSize(s, &ns));
993: PetscCall(VecGetArrayRead(v, &x));
994: PetscCall(VecGetArray(s, &y));
996: bs = v->map->bs;
997: 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);
998: x += start;
999: n = n / bs;
1001: if (addv == INSERT_VALUES) {
1002: for (i = 0; i < n; i++) y[i] = x[bs * i];
1003: } else if (addv == ADD_VALUES) {
1004: for (i = 0; i < n; i++) y[i] += x[bs * i];
1005: #if !PetscDefined(USE_COMPLEX)
1006: } else if (addv == MAX_VALUES) {
1007: for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
1008: #endif
1009: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1011: PetscCall(VecRestoreArrayRead(v, &x));
1012: PetscCall(VecRestoreArray(s, &y));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
1017: {
1018: PetscInt i, n, bs, ns;
1019: PetscScalar *x;
1020: const PetscScalar *y;
1022: PetscFunctionBegin;
1023: PetscCall(VecGetLocalSize(v, &n));
1024: PetscCall(VecGetLocalSize(s, &ns));
1025: PetscCall(VecGetArray(v, &x));
1026: PetscCall(VecGetArrayRead(s, &y));
1028: bs = v->map->bs;
1029: 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);
1030: x += start;
1031: n = n / bs;
1033: if (addv == INSERT_VALUES) {
1034: for (i = 0; i < n; i++) x[bs * i] = y[i];
1035: } else if (addv == ADD_VALUES) {
1036: for (i = 0; i < n; i++) x[bs * i] += y[i];
1037: #if !PetscDefined(USE_COMPLEX)
1038: } else if (addv == MAX_VALUES) {
1039: for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
1040: #endif
1041: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1043: PetscCall(VecRestoreArray(v, &x));
1044: PetscCall(VecRestoreArrayRead(s, &y));
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
1049: {
1050: PetscInt i, j, n, bs, bss, ns;
1051: const PetscScalar *x;
1052: PetscScalar *y;
1054: PetscFunctionBegin;
1055: PetscCall(VecGetLocalSize(v, &n));
1056: PetscCall(VecGetLocalSize(s, &ns));
1057: PetscCall(VecGetArrayRead(v, &x));
1058: PetscCall(VecGetArray(s, &y));
1060: bs = v->map->bs;
1061: bss = s->map->bs;
1062: n = n / bs;
1064: if (PetscDefined(USE_DEBUG)) {
1065: PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1066: for (j = 0; j < nidx; j++) {
1067: PetscCheck(idxv[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxv[j]);
1068: 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);
1069: }
1070: PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not gathering into all locations");
1071: }
1073: if (addv == INSERT_VALUES) {
1074: if (!idxs) {
1075: for (i = 0; i < n; i++) {
1076: for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
1077: }
1078: } else {
1079: for (i = 0; i < n; i++) {
1080: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
1081: }
1082: }
1083: } else if (addv == ADD_VALUES) {
1084: if (!idxs) {
1085: for (i = 0; i < n; i++) {
1086: for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
1087: }
1088: } else {
1089: for (i = 0; i < n; i++) {
1090: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
1091: }
1092: }
1093: #if !PetscDefined(USE_COMPLEX)
1094: } else if (addv == MAX_VALUES) {
1095: if (!idxs) {
1096: for (i = 0; i < n; i++) {
1097: for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1098: }
1099: } else {
1100: for (i = 0; i < n; i++) {
1101: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1102: }
1103: }
1104: #endif
1105: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1107: PetscCall(VecRestoreArrayRead(v, &x));
1108: PetscCall(VecRestoreArray(s, &y));
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
1113: {
1114: PetscInt j, i, n, bs, ns, bss;
1115: PetscScalar *x;
1116: const PetscScalar *y;
1118: PetscFunctionBegin;
1119: PetscCall(VecGetLocalSize(v, &n));
1120: PetscCall(VecGetLocalSize(s, &ns));
1121: PetscCall(VecGetArray(v, &x));
1122: PetscCall(VecGetArrayRead(s, &y));
1124: bs = v->map->bs;
1125: bss = s->map->bs;
1126: n = n / bs;
1128: if (PetscDefined(USE_DEBUG)) {
1129: PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1130: for (j = 0; j < bss; j++) {
1131: if (idxs) {
1132: PetscCheck(idxs[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxs[j]);
1133: 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);
1134: }
1135: }
1136: PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not scattering from all locations");
1137: }
1139: if (addv == INSERT_VALUES) {
1140: if (!idxs) {
1141: for (i = 0; i < n; i++) {
1142: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
1143: }
1144: } else {
1145: for (i = 0; i < n; i++) {
1146: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
1147: }
1148: }
1149: } else if (addv == ADD_VALUES) {
1150: if (!idxs) {
1151: for (i = 0; i < n; i++) {
1152: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
1153: }
1154: } else {
1155: for (i = 0; i < n; i++) {
1156: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
1157: }
1158: }
1159: #if !PetscDefined(USE_COMPLEX)
1160: } else if (addv == MAX_VALUES) {
1161: if (!idxs) {
1162: for (i = 0; i < n; i++) {
1163: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1164: }
1165: } else {
1166: for (i = 0; i < n; i++) {
1167: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1168: }
1169: }
1170: #endif
1171: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
1173: PetscCall(VecRestoreArray(v, &x));
1174: PetscCall(VecRestoreArrayRead(s, &y));
1175: PetscFunctionReturn(PETSC_SUCCESS);
1176: }
1178: static PetscErrorCode VecApplyUnary_Private(Vec v, PetscDeviceContext dctx, const char async_name[], PetscErrorCode (*unary_op)(Vec), PetscScalar (*UnaryFunc)(PetscScalar))
1179: {
1180: PetscFunctionBegin;
1182: PetscCall(VecSetErrorIfLocked(v, 1));
1183: if (dctx) {
1184: PetscErrorCode (*unary_op_async)(Vec, PetscDeviceContext);
1186: PetscCall(PetscObjectQueryFunction((PetscObject)v, async_name, &unary_op_async));
1187: if (unary_op_async) {
1188: PetscCall((*unary_op_async)(v, dctx));
1189: PetscFunctionReturn(PETSC_SUCCESS);
1190: }
1191: }
1192: if (unary_op) {
1194: PetscCall((*unary_op)(v));
1195: } else {
1196: PetscInt n;
1197: PetscScalar *x;
1200: PetscCall(VecGetLocalSize(v, &n));
1201: PetscCall(VecGetArray(v, &x));
1202: for (PetscInt i = 0; i < n; ++i) x[i] = UnaryFunc(x[i]);
1203: PetscCall(VecRestoreArray(v, &x));
1204: }
1205: PetscFunctionReturn(PETSC_SUCCESS);
1206: }
1208: static PetscScalar ScalarReciprocal_Function(PetscScalar x)
1209: {
1210: const PetscScalar zero = 0.0;
1212: return x == zero ? zero : ((PetscScalar)1.0) / x;
1213: }
1215: PetscErrorCode VecReciprocalAsync_Private(Vec v, PetscDeviceContext dctx)
1216: {
1217: PetscFunctionBegin;
1218: PetscCall(PetscLogEventBegin(VEC_Reciprocal, v, NULL, NULL, NULL));
1219: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Reciprocal), v->ops->reciprocal, ScalarReciprocal_Function));
1220: PetscCall(PetscLogEventEnd(VEC_Reciprocal, v, NULL, NULL, NULL));
1221: PetscFunctionReturn(PETSC_SUCCESS);
1222: }
1224: PetscErrorCode VecReciprocal_Default(Vec v)
1225: {
1226: PetscFunctionBegin;
1227: PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarReciprocal_Function));
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: static PetscScalar ScalarExp_Function(PetscScalar x)
1232: {
1233: return PetscExpScalar(x);
1234: }
1236: PetscErrorCode VecExpAsync_Private(Vec v, PetscDeviceContext dctx)
1237: {
1238: PetscFunctionBegin;
1240: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Exp), v->ops->exp, ScalarExp_Function));
1241: PetscFunctionReturn(PETSC_SUCCESS);
1242: }
1244: /*@
1245: VecExp - Replaces each component of a vector by e^x_i
1247: Not Collective
1249: Input Parameter:
1250: . v - The vector
1252: Output Parameter:
1253: . v - The vector of exponents
1255: Level: beginner
1257: .seealso: `Vec`, `VecLog()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1258: @*/
1259: PetscErrorCode VecExp(Vec v)
1260: {
1261: PetscFunctionBegin;
1262: PetscCall(VecExpAsync_Private(v, NULL));
1263: PetscFunctionReturn(PETSC_SUCCESS);
1264: }
1266: static PetscScalar ScalarLog_Function(PetscScalar x)
1267: {
1268: return PetscLogScalar(x);
1269: }
1271: PetscErrorCode VecLogAsync_Private(Vec v, PetscDeviceContext dctx)
1272: {
1273: PetscFunctionBegin;
1275: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Log), v->ops->log, ScalarLog_Function));
1276: PetscFunctionReturn(PETSC_SUCCESS);
1277: }
1279: /*@
1280: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1282: Not Collective
1284: Input Parameter:
1285: . v - The vector
1287: Output Parameter:
1288: . v - The vector of logs
1290: Level: beginner
1292: .seealso: `Vec`, `VecExp()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1293: @*/
1294: PetscErrorCode VecLog(Vec v)
1295: {
1296: PetscFunctionBegin;
1297: PetscCall(VecLogAsync_Private(v, NULL));
1298: PetscFunctionReturn(PETSC_SUCCESS);
1299: }
1301: static PetscScalar ScalarAbs_Function(PetscScalar x)
1302: {
1303: return PetscAbsScalar(x);
1304: }
1306: PetscErrorCode VecAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1307: {
1308: PetscFunctionBegin;
1310: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Abs), v->ops->abs, ScalarAbs_Function));
1311: PetscFunctionReturn(PETSC_SUCCESS);
1312: }
1314: /*@
1315: VecAbs - Replaces every element in a vector with its absolute value.
1317: Logically Collective
1319: Input Parameter:
1320: . v - the vector
1322: Level: intermediate
1324: .seealso: `Vec`, `VecExp()`, `VecSqrtAbs()`, `VecReciprocal()`, `VecLog()`, `VecPointwiseSign()`
1325: @*/
1326: PetscErrorCode VecAbs(Vec v)
1327: {
1328: PetscFunctionBegin;
1329: PetscCall(VecAbsAsync_Private(v, NULL));
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }
1333: static PetscScalar ScalarConjugate_Function(PetscScalar x)
1334: {
1335: return PetscConj(x);
1336: }
1338: PetscErrorCode VecConjugateAsync_Private(Vec v, PetscDeviceContext dctx)
1339: {
1340: PetscFunctionBegin;
1342: if (PetscDefined(USE_COMPLEX)) PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Conjugate), v->ops->conjugate, ScalarConjugate_Function));
1343: PetscFunctionReturn(PETSC_SUCCESS);
1344: }
1346: /*@
1347: VecConjugate - Conjugates a vector. That is, replace every entry in a vector with its complex conjugate
1349: Logically Collective
1351: Input Parameter:
1352: . x - the vector
1354: Level: intermediate
1356: .seealso: [](ch_vectors), `Vec`, `VecSet()`
1357: @*/
1358: PetscErrorCode VecConjugate(Vec x)
1359: {
1360: PetscFunctionBegin;
1361: PetscCall(VecConjugateAsync_Private(x, NULL));
1362: PetscFunctionReturn(PETSC_SUCCESS);
1363: }
1365: static PetscScalar ScalarSqrtAbs_Function(PetscScalar x)
1366: {
1367: return PetscSqrtScalar(ScalarAbs_Function(x));
1368: }
1370: PetscErrorCode VecSqrtAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1371: {
1372: PetscFunctionBegin;
1374: PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(SqrtAbs), v->ops->sqrt, ScalarSqrtAbs_Function));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: /*@
1379: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1381: Not Collective
1383: Input Parameter:
1384: . v - The vector
1386: Level: beginner
1388: Note:
1389: The actual function is sqrt(|x_i|)
1391: .seealso: `Vec`, `VecLog()`, `VecExp()`, `VecReciprocal()`, `VecAbs()`
1392: @*/
1393: PetscErrorCode VecSqrtAbs(Vec v)
1394: {
1395: PetscFunctionBegin;
1396: PetscCall(VecSqrtAbsAsync_Private(v, NULL));
1397: PetscFunctionReturn(PETSC_SUCCESS);
1398: }
1400: #if PetscDefined(USE_COMPLEX)
1401: static PetscScalar ScalarImaginaryPart_Function(PetscScalar x)
1402: {
1403: const PetscReal imag = PetscImaginaryPart(x);
1405: return PetscCMPLX(imag, 0.0);
1406: }
1407: #endif
1409: /*@
1410: VecImaginaryPart - Replaces a complex vector with its imaginary part
1412: Collective
1414: Input Parameter:
1415: . v - the vector
1417: Level: beginner
1419: .seealso: `Vec`, `VecNorm()`, `VecRealPart()`
1420: @*/
1421: PetscErrorCode VecImaginaryPart(Vec v)
1422: {
1423: PetscFunctionBegin;
1425: #if PetscDefined(USE_COMPLEX)
1426: PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarImaginaryPart_Function));
1427: #else
1428: PetscCall(VecZeroEntries(v));
1429: #endif
1430: PetscFunctionReturn(PETSC_SUCCESS);
1431: }
1433: #if PetscDefined(USE_COMPLEX)
1434: static PetscScalar ScalarRealPart_Function(PetscScalar x)
1435: {
1436: const PetscReal real = PetscRealPart(x);
1438: return PetscCMPLX(real, 0.0);
1439: }
1440: #endif
1442: /*@
1443: VecRealPart - Replaces a complex vector with its real part
1445: Collective
1447: Input Parameter:
1448: . v - the vector
1450: Level: beginner
1452: .seealso: `Vec`, `VecNorm()`, `VecImaginaryPart()`
1453: @*/
1454: PetscErrorCode VecRealPart(Vec v)
1455: {
1456: PetscFunctionBegin;
1458: #if PetscDefined(USE_COMPLEX)
1459: PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarRealPart_Function));
1460: #endif
1461: PetscFunctionReturn(PETSC_SUCCESS);
1462: }
1464: /*@
1465: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1467: Collective
1469: Input Parameters:
1470: + s - first vector
1471: - t - second vector
1473: Output Parameters:
1474: + dp - s'conj(t)
1475: - nm - t'conj(t)
1477: Level: advanced
1479: Note:
1480: conj(x) is the complex conjugate of x when x is complex
1482: .seealso: `Vec`, `VecDot()`, `VecNorm()`, `VecDotBegin()`, `VecNormBegin()`, `VecDotEnd()`, `VecNormEnd()`
1483: @*/
1484: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1485: {
1486: PetscScalar work[] = {0.0, 0.0};
1488: PetscFunctionBegin;
1491: PetscAssertPointer(dp, 3);
1492: PetscAssertPointer(nm, 4);
1495: PetscCheckSameTypeAndComm(s, 1, t, 2);
1496: PetscCheck(s->map->N == t->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector global lengths");
1497: PetscCheck(s->map->n == t->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector local lengths");
1499: PetscCall(PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0));
1500: if (s->ops->dotnorm2) {
1501: PetscUseTypeMethod(s, dotnorm2, t, work, work + 1);
1502: } else {
1503: const PetscScalar *sx, *tx;
1504: PetscInt n;
1506: PetscCall(VecGetLocalSize(s, &n));
1507: PetscCall(VecGetArrayRead(s, &sx));
1508: PetscCall(VecGetArrayRead(t, &tx));
1509: for (PetscInt i = 0; i < n; ++i) {
1510: const PetscScalar txconj = PetscConj(tx[i]);
1512: work[0] += sx[i] * txconj;
1513: work[1] += tx[i] * txconj;
1514: }
1515: PetscCall(VecRestoreArrayRead(t, &tx));
1516: PetscCall(VecRestoreArrayRead(s, &sx));
1517: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, work, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
1518: PetscCall(PetscLogFlops(4.0 * n));
1519: }
1520: PetscCall(PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0));
1521: *dp = work[0];
1522: *nm = PetscRealPart(work[1]);
1523: PetscFunctionReturn(PETSC_SUCCESS);
1524: }
1526: /*@
1527: VecSum - Computes the sum of all the components of a vector.
1529: Collective
1531: Input Parameter:
1532: . v - the vector
1534: Output Parameter:
1535: . sum - the result
1537: Level: beginner
1539: .seealso: `Vec`, `VecMean()`, `VecNorm()`
1540: @*/
1541: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1542: {
1543: PetscScalar tmp = 0.0;
1545: PetscFunctionBegin;
1547: PetscAssertPointer(sum, 2);
1548: if (v->ops->sum) PetscUseTypeMethod(v, sum, &tmp);
1549: else {
1550: const PetscScalar *x;
1551: PetscInt n;
1553: PetscCall(VecGetLocalSize(v, &n));
1554: PetscCall(VecGetArrayRead(v, &x));
1555: for (PetscInt i = 0; i < n; ++i) tmp += x[i];
1556: PetscCall(VecRestoreArrayRead(v, &x));
1557: }
1558: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &tmp, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
1559: *sum = tmp;
1560: PetscFunctionReturn(PETSC_SUCCESS);
1561: }
1563: /*@
1564: VecMean - Computes the arithmetic mean of all the components of a vector.
1566: Collective
1568: Input Parameter:
1569: . v - the vector
1571: Output Parameter:
1572: . mean - the result
1574: Level: beginner
1576: .seealso: `Vec`, `VecSum()`, `VecNorm()`
1577: @*/
1578: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1579: {
1580: PetscInt n;
1582: PetscFunctionBegin;
1584: PetscAssertPointer(mean, 2);
1585: PetscCall(VecGetSize(v, &n));
1586: PetscCall(VecSum(v, mean));
1587: *mean /= n;
1588: PetscFunctionReturn(PETSC_SUCCESS);
1589: }
1591: PetscErrorCode VecShiftAsync_Private(Vec v, PetscScalar shift, PetscDeviceContext dctx)
1592: {
1593: PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext) = NULL;
1595: PetscFunctionBegin;
1596: if (dctx) {
1597: PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext);
1599: PetscCall(PetscObjectQueryFunction((PetscObject)v, VecAsyncFnName(Shift), &shift_async));
1600: }
1601: if (shift_async) PetscCall((*shift_async)(v, shift, dctx));
1602: else if (v->ops->shift) PetscUseTypeMethod(v, shift, shift);
1603: else {
1604: PetscInt n;
1605: PetscScalar *x;
1607: PetscCall(VecGetLocalSize(v, &n));
1608: PetscCall(VecGetArray(v, &x));
1609: for (PetscInt i = 0; i < n; ++i) x[i] += shift;
1610: PetscCall(VecRestoreArray(v, &x));
1611: PetscCall(PetscLogFlops(n));
1612: }
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: /*@
1617: VecShift - Shifts all of the components of a vector by computing
1618: `x[i] = x[i] + shift`.
1620: Logically Collective
1622: Input Parameters:
1623: + v - the vector
1624: - shift - the shift
1626: Level: intermediate
1628: .seealso: `Vec`, `VecISShift()`
1629: @*/
1630: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1631: {
1632: PetscFunctionBegin;
1635: PetscCall(VecSetErrorIfLocked(v, 1));
1636: if (shift == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
1637: PetscCall(PetscLogEventBegin(VEC_Shift, v, 0, 0, 0));
1638: PetscCall(VecShiftAsync_Private(v, shift, NULL));
1639: PetscCall(PetscLogEventEnd(VEC_Shift, v, 0, 0, 0));
1640: PetscFunctionReturn(PETSC_SUCCESS);
1641: }
1643: /*@
1644: VecPermute - Permutes a vector in place using the given ordering.
1646: Input Parameters:
1647: + x - The vector
1648: . row - The ordering
1649: - inv - The flag for inverting the permutation
1651: Level: beginner
1653: Note:
1654: This function does not yet support parallel Index Sets with non-local permutations
1656: .seealso: `Vec`, `MatPermute()`
1657: @*/
1658: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1659: {
1660: PetscScalar *array, *newArray;
1661: const PetscInt *idx;
1662: PetscInt i, rstart, rend;
1664: PetscFunctionBegin;
1667: PetscCall(VecSetErrorIfLocked(x, 1));
1668: PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1669: PetscCall(ISGetIndices(row, &idx));
1670: PetscCall(VecGetArray(x, &array));
1671: PetscCall(PetscMalloc1(x->map->n, &newArray));
1672: PetscCall(PetscArraycpy(newArray, array, x->map->n));
1673: if (PetscDefined(USE_DEBUG)) {
1674: 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]);
1675: }
1676: if (!inv) {
1677: for (i = 0; i < x->map->n; i++) array[i] = newArray[idx[i] - rstart];
1678: } else {
1679: for (i = 0; i < x->map->n; i++) array[idx[i] - rstart] = newArray[i];
1680: }
1681: PetscCall(VecRestoreArray(x, &array));
1682: PetscCall(ISRestoreIndices(row, &idx));
1683: PetscCall(PetscFree(newArray));
1684: PetscFunctionReturn(PETSC_SUCCESS);
1685: }
1687: /*@
1688: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1689: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1690: Does NOT take round-off errors into account.
1692: Collective
1694: Input Parameters:
1695: + vec1 - the first vector
1696: - vec2 - the second vector
1698: Output Parameter:
1699: . flg - `PETSC_TRUE` if the vectors are equal; `PETSC_FALSE` otherwise.
1701: Level: intermediate
1703: .seealso: `Vec`
1704: @*/
1705: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1706: {
1707: const PetscScalar *v1, *v2;
1708: PetscInt n1, n2, N1, N2;
1710: PetscFunctionBegin;
1713: PetscAssertPointer(flg, 3);
1714: if (vec1 == vec2) *flg = PETSC_TRUE;
1715: else {
1716: PetscCall(VecGetSize(vec1, &N1));
1717: PetscCall(VecGetSize(vec2, &N2));
1718: if (N1 != N2) *flg = PETSC_FALSE;
1719: else {
1720: PetscCall(VecGetLocalSize(vec1, &n1));
1721: PetscCall(VecGetLocalSize(vec2, &n2));
1722: if (n1 != n2) *flg = PETSC_FALSE;
1723: else {
1724: PetscCall(VecGetArrayRead(vec1, &v1));
1725: PetscCall(VecGetArrayRead(vec2, &v2));
1726: PetscCall(PetscArraycmp(v1, v2, n1, flg));
1727: PetscCall(VecRestoreArrayRead(vec1, &v1));
1728: PetscCall(VecRestoreArrayRead(vec2, &v2));
1729: }
1730: }
1731: /* combine results from all processors */
1732: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)vec1)));
1733: }
1734: PetscFunctionReturn(PETSC_SUCCESS);
1735: }
1737: /*@
1738: VecUniqueEntries - Compute the number of unique entries, and those entries
1740: Collective
1742: Input Parameter:
1743: . vec - the vector
1745: Output Parameters:
1746: + n - The number of unique entries
1747: - e - The entries, each MPI process receives all the unique entries
1749: Level: intermediate
1751: .seealso: `Vec`
1752: @*/
1753: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar *e[])
1754: {
1755: const PetscScalar *v;
1756: PetscScalar *tmp, *vals;
1757: PetscMPIInt *N, *displs, l;
1758: PetscInt ng, m, i, j, p;
1759: PetscMPIInt size;
1761: PetscFunctionBegin;
1763: PetscAssertPointer(n, 2);
1764: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1765: PetscCall(VecGetLocalSize(vec, &m));
1766: PetscCall(VecGetArrayRead(vec, &v));
1767: PetscCall(PetscMalloc2(m, &tmp, size, &N));
1768: for (i = 0, l = 0; i < m; ++i) {
1769: /* Can speed this up with sorting */
1770: for (j = 0; j < l; ++j) {
1771: if (v[i] == tmp[j]) break;
1772: }
1773: if (j == l) {
1774: tmp[j] = v[i];
1775: ++l;
1776: }
1777: }
1778: PetscCall(VecRestoreArrayRead(vec, &v));
1779: /* Gather serial results */
1780: PetscCallMPI(MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec)));
1781: for (p = 0, ng = 0; p < size; ++p) ng += N[p];
1782: PetscCall(PetscMalloc2(ng, &vals, size + 1, &displs));
1783: for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
1784: PetscCallMPI(MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec)));
1785: /* Find unique entries */
1786: #if PetscDefined(USE_COMPLEX)
1787: SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1788: #else
1789: *n = displs[size];
1790: PetscCall(PetscSortRemoveDupsReal(n, vals));
1791: if (e) {
1792: PetscAssertPointer(e, 3);
1793: PetscCall(PetscMalloc1(*n, e));
1794: for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
1795: }
1796: PetscCall(PetscFree2(vals, displs));
1797: PetscCall(PetscFree2(tmp, N));
1798: PetscFunctionReturn(PETSC_SUCCESS);
1799: #endif
1800: }