Actual source code: vinv.c
2: /*
3: Some useful vector utility functions.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
7: /*@
8: VecStrideSet - Sets a subvector of a vector defined
9: by a starting point and a stride with a given value
11: Logically Collective
13: Input Parameters:
14: + v - the vector
15: . start - starting point of the subvector (defined by a stride)
16: - s - value to set for each entry in that subvector
18: Level: advanced
20: Notes:
21: One must call `VecSetBlockSize()` before this routine to set the stride
22: information, or use a vector created from a multicomponent `DMDA`.
24: This will only work if the desire subvector is a stride subvector
26: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideScale()`
27: @*/
28: PetscErrorCode VecStrideSet(Vec v, PetscInt start, PetscScalar s)
29: {
30: PetscInt i, n, bs;
31: PetscScalar *x;
33: 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\n 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()`, `VecStrideScale()`
65: @*/
66: PetscErrorCode VecStrideScale(Vec v, PetscInt start, PetscScalar scale)
67: {
68: PetscInt i, n, bs;
69: PetscScalar *x;
71: PetscFunctionBegin;
73: PetscCall(VecGetLocalSize(v, &n));
74: PetscCall(VecGetBlockSize(v, &bs));
75: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
76: 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);
77: PetscCall(VecGetArray(v, &x));
78: for (i = start; i < n; i += bs) x[i] *= scale;
79: PetscCall(VecRestoreArray(v, &x));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: /*@
84: VecStrideNorm - Computes the norm of subvector of a vector defined
85: by a starting point and a stride.
87: Collective
89: Input Parameters:
90: + v - the vector
91: . start - starting point of the subvector (defined by a stride)
92: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`
94: Output Parameter:
95: . norm - the norm
97: Level: advanced
99: Notes:
100: One must call `VecSetBlockSize()` before this routine to set the stride
101: information, or use a vector created from a multicomponent `DMDA`.
103: If x is the array representing the vector x then this computes the norm
104: of the array (x[start],x[start+stride],x[start+2*stride], ....)
106: This is useful for computing, say the norm of the pressure variable when
107: the pressure is stored (interlaced) with other variables, say density etc.
109: This will only work if the desire subvector is a stride subvector
111: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
112: @*/
113: PetscErrorCode VecStrideNorm(Vec v, PetscInt start, NormType ntype, PetscReal *nrm)
114: {
115: PetscInt i, n, bs;
116: const PetscScalar *x;
117: PetscReal tnorm;
119: PetscFunctionBegin;
123: PetscCall(VecGetLocalSize(v, &n));
124: PetscCall(VecGetBlockSize(v, &bs));
125: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
126: 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);
127: PetscCall(VecGetArrayRead(v, &x));
128: if (ntype == NORM_2) {
129: PetscScalar sum = 0.0;
130: for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
131: tnorm = PetscRealPart(sum);
132: PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
133: *nrm = PetscSqrtReal(*nrm);
134: } else if (ntype == NORM_1) {
135: tnorm = 0.0;
136: for (i = start; i < n; i += bs) tnorm += PetscAbsScalar(x[i]);
137: PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
138: } else if (ntype == NORM_INFINITY) {
139: tnorm = 0.0;
140: for (i = start; i < n; i += bs) {
141: if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
142: }
143: PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
144: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
145: PetscCall(VecRestoreArrayRead(v, &x));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: VecStrideMax - Computes the maximum of subvector of a vector defined
151: by a starting point and a stride and optionally its location.
153: Collective
155: Input Parameters:
156: + v - the vector
157: - start - starting point of the subvector (defined by a stride)
159: Output Parameters:
160: + idex - the location where the maximum occurred (pass `NULL` if not required)
161: - nrm - the maximum value in the subvector
163: Level: advanced
165: Notes:
166: One must call `VecSetBlockSize()` before this routine to set the stride
167: information, or use a vector created from a multicomponent `DMDA`.
169: If xa is the array representing the vector x, then this computes the max
170: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
172: This is useful for computing, say the maximum of the pressure variable when
173: the pressure is stored (interlaced) with other variables, e.g., density, etc.
174: This will only work if the desire subvector is a stride subvector.
176: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
177: @*/
178: PetscErrorCode VecStrideMax(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
179: {
180: PetscInt i, n, bs, id = -1;
181: const PetscScalar *x;
182: PetscReal max = PETSC_MIN_REAL;
184: PetscFunctionBegin;
187: PetscCall(VecGetLocalSize(v, &n));
188: PetscCall(VecGetBlockSize(v, &bs));
189: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
190: 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);
191: PetscCall(VecGetArrayRead(v, &x));
192: for (i = start; i < n; i += bs) {
193: if (PetscRealPart(x[i]) > max) {
194: max = PetscRealPart(x[i]);
195: id = i;
196: }
197: }
198: PetscCall(VecRestoreArrayRead(v, &x));
199: #if defined(PETSC_HAVE_MPIUNI)
200: *nrm = max;
201: if (idex) *idex = id;
202: #else
203: if (!idex) {
204: PetscCall(MPIU_Allreduce(&max, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
205: } else {
206: struct {
207: PetscReal v;
208: PetscInt i;
209: } in, out;
210: PetscInt rstart;
212: PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
213: in.v = max;
214: in.i = rstart + id;
215: PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v)));
216: *nrm = out.v;
217: *idex = out.i;
218: }
219: #endif
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*@
224: VecStrideMin - Computes the minimum of subvector of a vector defined
225: by a starting point and a stride and optionally its location.
227: Collective
229: Input Parameters:
230: + v - the vector
231: - start - starting point of the subvector (defined by a stride)
233: Output Parameters:
234: + idex - the location where the minimum occurred. (pass `NULL` if not required)
235: - nrm - the minimum value in the subvector
237: Level: advanced
239: Notes:
240: One must call `VecSetBlockSize()` before this routine to set the stride
241: information, or use a vector created from a multicomponent `DMDA`.
243: If xa is the array representing the vector x, then this computes the min
244: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
246: This is useful for computing, say the minimum of the pressure variable when
247: the pressure is stored (interlaced) with other variables, e.g., density, etc.
248: This will only work if the desire subvector is a stride subvector.
250: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
251: @*/
252: PetscErrorCode VecStrideMin(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
253: {
254: PetscInt i, n, bs, id = -1;
255: const PetscScalar *x;
256: PetscReal min = PETSC_MAX_REAL;
258: PetscFunctionBegin;
261: PetscCall(VecGetLocalSize(v, &n));
262: PetscCall(VecGetBlockSize(v, &bs));
263: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
264: 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);
265: PetscCall(VecGetArrayRead(v, &x));
266: for (i = start; i < n; i += bs) {
267: if (PetscRealPart(x[i]) < min) {
268: min = PetscRealPart(x[i]);
269: id = i;
270: }
271: }
272: PetscCall(VecRestoreArrayRead(v, &x));
273: #if defined(PETSC_HAVE_MPIUNI)
274: *nrm = min;
275: if (idex) *idex = id;
276: #else
277: if (!idex) {
278: PetscCall(MPIU_Allreduce(&min, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v)));
279: } else {
280: struct {
281: PetscReal v;
282: PetscInt i;
283: } in, out;
284: PetscInt rstart;
286: PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
287: in.v = min;
288: in.i = rstart + id;
289: PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v)));
290: *nrm = out.v;
291: *idex = out.i;
292: }
293: #endif
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@
298: VecStrideSum - Computes the sum of subvector of a vector defined
299: by a starting point and a stride.
301: Collective
303: Input Parameters:
304: + v - the vector
305: - start - starting point of the subvector (defined by a stride)
307: Output Parameter:
308: . sum - the sum
310: Level: advanced
312: Notes:
313: One must call `VecSetBlockSize()` before this routine to set the stride
314: information, or use a vector created from a multicomponent `DMDA`.
316: If x is the array representing the vector x then this computes the sum
317: of the array (x[start],x[start+stride],x[start+2*stride], ....)
319: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
320: @*/
321: PetscErrorCode VecStrideSum(Vec v, PetscInt start, PetscScalar *sum)
322: {
323: PetscInt i, n, bs;
324: const PetscScalar *x;
325: PetscScalar local_sum = 0.0;
327: PetscFunctionBegin;
330: PetscCall(VecGetLocalSize(v, &n));
331: PetscCall(VecGetBlockSize(v, &bs));
332: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
333: 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);
334: PetscCall(VecGetArrayRead(v, &x));
335: for (i = start; i < n; i += bs) local_sum += x[i];
336: PetscCall(MPIU_Allreduce(&local_sum, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
337: PetscCall(VecRestoreArrayRead(v, &x));
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: VecStrideScaleAll - Scales the subvectors of a vector defined
343: by a starting point and a stride.
345: Logically Collective
347: Input Parameters:
348: + v - the vector
349: - scales - values to multiply each subvector entry by
351: Level: advanced
353: Notes:
354: One must call `VecSetBlockSize()` before this routine to set the stride
355: information, or use a vector created from a multicomponent `DMDA`.
357: The dimension of scales must be the same as the vector block size
359: .seealso: `Vec`, `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
360: @*/
361: PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
362: {
363: PetscInt i, j, n, bs;
364: PetscScalar *x;
366: PetscFunctionBegin;
369: PetscCall(VecGetLocalSize(v, &n));
370: PetscCall(VecGetBlockSize(v, &bs));
371: PetscCall(VecGetArray(v, &x));
372: /* need to provide optimized code for each bs */
373: for (i = 0; i < n; i += bs) {
374: for (j = 0; j < bs; j++) x[i + j] *= scales[j];
375: }
376: PetscCall(VecRestoreArray(v, &x));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*@
381: VecStrideNormAll - Computes the norms of subvectors of a vector defined
382: by a starting point and a stride.
384: Collective
386: Input Parameters:
387: + v - the vector
388: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`
390: Output Parameter:
391: . nrm - the norms
393: Level: advanced
395: Notes:
396: One must call `VecSetBlockSize()` before this routine to set the stride
397: information, or use a vector created from a multicomponent `DMDA`.
399: If x is the array representing the vector x then this computes the norm
400: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
402: The dimension of nrm must be the same as the vector block size
404: This will only work if the desire subvector is a stride subvector
406: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
407: @*/
408: PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
409: {
410: PetscInt i, j, n, bs;
411: const PetscScalar *x;
412: PetscReal tnorm[128];
413: MPI_Comm comm;
415: PetscFunctionBegin;
419: PetscCall(VecGetLocalSize(v, &n));
420: PetscCall(VecGetArrayRead(v, &x));
421: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
423: PetscCall(VecGetBlockSize(v, &bs));
424: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
426: if (ntype == NORM_2) {
427: PetscScalar sum[128];
428: for (j = 0; j < bs; j++) sum[j] = 0.0;
429: for (i = 0; i < n; i += bs) {
430: for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
431: }
432: for (j = 0; j < bs; j++) tnorm[j] = PetscRealPart(sum[j]);
434: PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
435: for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
436: } else if (ntype == NORM_1) {
437: for (j = 0; j < bs; j++) tnorm[j] = 0.0;
439: for (i = 0; i < n; i += bs) {
440: for (j = 0; j < bs; j++) tnorm[j] += PetscAbsScalar(x[i + j]);
441: }
443: PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
444: } else if (ntype == NORM_INFINITY) {
445: PetscReal tmp;
446: for (j = 0; j < bs; j++) tnorm[j] = 0.0;
448: for (i = 0; i < n; i += bs) {
449: for (j = 0; j < bs; j++) {
450: if ((tmp = PetscAbsScalar(x[i + j])) > tnorm[j]) tnorm[j] = tmp;
451: /* check special case of tmp == NaN */
452: if (tmp != tmp) {
453: tnorm[j] = tmp;
454: break;
455: }
456: }
457: }
458: PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_MAX, comm));
459: } else SETERRQ(comm, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
460: PetscCall(VecRestoreArrayRead(v, &x));
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
466: by a starting point and a stride and optionally its location.
468: Collective
470: Input Parameter:
471: . v - the vector
473: Output Parameters:
474: + index - the location where the maximum occurred (not supported, pass `NULL`,
475: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
476: - nrm - the maximum values of each subvector
478: Level: advanced
480: Notes:
481: One must call `VecSetBlockSize()` before this routine to set the stride
482: information, or use a vector created from a multicomponent `DMDA`.
484: The dimension of nrm must be the same as the vector block size
486: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
487: @*/
488: PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
489: {
490: PetscInt i, j, n, bs;
491: const PetscScalar *x;
492: PetscReal max[128], tmp;
493: MPI_Comm comm;
495: PetscFunctionBegin;
498: 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");
499: PetscCall(VecGetLocalSize(v, &n));
500: PetscCall(VecGetArrayRead(v, &x));
501: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
503: PetscCall(VecGetBlockSize(v, &bs));
504: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
506: if (!n) {
507: for (j = 0; j < bs; j++) max[j] = PETSC_MIN_REAL;
508: } else {
509: for (j = 0; j < bs; j++) max[j] = PetscRealPart(x[j]);
511: for (i = bs; i < n; i += bs) {
512: for (j = 0; j < bs; j++) {
513: if ((tmp = PetscRealPart(x[i + j])) > max[j]) max[j] = tmp;
514: }
515: }
516: }
517: PetscCall(MPIU_Allreduce(max, nrm, bs, MPIU_REAL, MPIU_MAX, comm));
519: PetscCall(VecRestoreArrayRead(v, &x));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: /*@
524: VecStrideMinAll - Computes the minimum of subvector of a vector defined
525: by a starting point and a stride and optionally its location.
527: Collective
529: Input Parameter:
530: . v - the vector
532: Output Parameters:
533: + idex - the location where the minimum occurred (not supported, pass `NULL`,
534: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
535: - nrm - the minimums of each subvector
537: Level: advanced
539: Notes:
540: One must call `VecSetBlockSize()` before this routine to set the stride
541: information, or use a vector created from a multicomponent `DMDA`.
543: The dimension of `nrm` must be the same as the vector block size
545: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
546: @*/
547: PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
548: {
549: PetscInt i, n, bs, j;
550: const PetscScalar *x;
551: PetscReal min[128], tmp;
552: MPI_Comm comm;
554: PetscFunctionBegin;
557: 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");
558: PetscCall(VecGetLocalSize(v, &n));
559: PetscCall(VecGetArrayRead(v, &x));
560: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
562: PetscCall(VecGetBlockSize(v, &bs));
563: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
565: if (!n) {
566: for (j = 0; j < bs; j++) min[j] = PETSC_MAX_REAL;
567: } else {
568: for (j = 0; j < bs; j++) min[j] = PetscRealPart(x[j]);
570: for (i = bs; i < n; i += bs) {
571: for (j = 0; j < bs; j++) {
572: if ((tmp = PetscRealPart(x[i + j])) < min[j]) min[j] = tmp;
573: }
574: }
575: }
576: PetscCall(MPIU_Allreduce(min, nrm, bs, MPIU_REAL, MPIU_MIN, comm));
578: PetscCall(VecRestoreArrayRead(v, &x));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*@
583: VecStrideSumAll - Computes the sums of subvectors of a vector defined by a stride.
585: Collective
587: Input Parameter:
588: . v - the vector
590: Output Parameter:
591: . sums - the sums
593: Level: advanced
595: Notes:
596: One must call `VecSetBlockSize()` before this routine to set the stride
597: information, or use a vector created from a multicomponent `DMDA`.
599: If x is the array representing the vector x then this computes the sum
600: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
602: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
603: @*/
604: PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
605: {
606: PetscInt i, j, n, bs;
607: const PetscScalar *x;
608: PetscScalar local_sums[128];
609: MPI_Comm comm;
611: PetscFunctionBegin;
614: PetscCall(VecGetLocalSize(v, &n));
615: PetscCall(VecGetArrayRead(v, &x));
616: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
618: PetscCall(VecGetBlockSize(v, &bs));
619: PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");
621: for (j = 0; j < bs; j++) local_sums[j] = 0.0;
622: for (i = 0; i < n; i += bs) {
623: for (j = 0; j < bs; j++) local_sums[j] += x[i + j];
624: }
625: PetscCall(MPIU_Allreduce(local_sums, sums, bs, MPIU_SCALAR, MPIU_SUM, comm));
627: PetscCall(VecRestoreArrayRead(v, &x));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: /*----------------------------------------------------------------------------------------------*/
632: /*@
633: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
634: separate vectors.
636: Collective
638: Input Parameters:
639: + v - the vector
640: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
642: Output Parameter:
643: . s - the location where the subvectors are stored
645: Level: advanced
647: Notes:
648: One must call `VecSetBlockSize()` before this routine to set the stride
649: information, or use a vector created from a multicomponent `DMDA`.
651: If x is the array representing the vector x then this gathers
652: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
653: for start=0,1,2,...bs-1
655: The parallel layout of the vector and the subvector must be the same;
656: i.e., nlocal of v = stride*(nlocal of s)
658: Not optimized; could be easily
660: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
661: `VecStrideScatterAll()`
662: @*/
663: PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
664: {
665: PetscInt i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
666: PetscScalar **y;
667: const PetscScalar *x;
669: PetscFunctionBegin;
673: PetscCall(VecGetLocalSize(v, &n));
674: PetscCall(VecGetLocalSize(s[0], &n2));
675: PetscCall(VecGetArrayRead(v, &x));
676: PetscCall(VecGetBlockSize(v, &bs));
677: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");
679: PetscCall(PetscMalloc2(bs, &y, bs, &bss));
680: nv = 0;
681: nvc = 0;
682: for (i = 0; i < bs; i++) {
683: PetscCall(VecGetBlockSize(s[i], &bss[i]));
684: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
685: PetscCall(VecGetArray(s[i], &y[i]));
686: nvc += bss[i];
687: nv++;
688: PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
689: if (nvc == bs) break;
690: }
692: n = n / bs;
694: jj = 0;
695: if (addv == INSERT_VALUES) {
696: for (j = 0; j < nv; j++) {
697: for (k = 0; k < bss[j]; k++) {
698: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
699: }
700: jj += bss[j];
701: }
702: } else if (addv == ADD_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: #if !defined(PETSC_USE_COMPLEX)
710: } else if (addv == MAX_VALUES) {
711: for (j = 0; j < nv; j++) {
712: for (k = 0; k < bss[j]; k++) {
713: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
714: }
715: jj += bss[j];
716: }
717: #endif
718: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
720: PetscCall(VecRestoreArrayRead(v, &x));
721: for (i = 0; i < nv; i++) PetscCall(VecRestoreArray(s[i], &y[i]));
723: PetscCall(PetscFree2(y, bss));
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: /*@
728: VecStrideScatterAll - Scatters all the single components from separate vectors into
729: a multi-component vector.
731: Collective
733: Input Parameters:
734: + s - the location where the subvectors are stored
735: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
737: Output Parameter:
738: . v - the multicomponent vector
740: Level: advanced
742: Notes:
743: One must call `VecSetBlockSize()` before this routine to set the stride
744: information, or use a vector created from a multicomponent `DMDA`.
746: The parallel layout of the vector and the subvector must be the same;
747: i.e., nlocal of v = stride*(nlocal of s)
749: Not optimized; could be easily
751: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
752: `VecStrideScatterAll()`
753: @*/
754: PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
755: {
756: PetscInt i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
757: PetscScalar *x;
758: PetscScalar const **y;
760: PetscFunctionBegin;
764: PetscCall(VecGetLocalSize(v, &n));
765: PetscCall(VecGetLocalSize(s[0], &n2));
766: PetscCall(VecGetArray(v, &x));
767: PetscCall(VecGetBlockSize(v, &bs));
768: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");
770: PetscCall(PetscMalloc2(bs, (PetscScalar ***)&y, bs, &bss));
771: nv = 0;
772: nvc = 0;
773: for (i = 0; i < bs; i++) {
774: PetscCall(VecGetBlockSize(s[i], &bss[i]));
775: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
776: PetscCall(VecGetArrayRead(s[i], &y[i]));
777: nvc += bss[i];
778: nv++;
779: PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
780: if (nvc == bs) break;
781: }
783: n = n / bs;
785: jj = 0;
786: if (addv == INSERT_VALUES) {
787: for (j = 0; j < nv; j++) {
788: for (k = 0; k < bss[j]; k++) {
789: for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
790: }
791: jj += bss[j];
792: }
793: } else if (addv == ADD_VALUES) {
794: for (j = 0; j < nv; j++) {
795: for (k = 0; k < bss[j]; k++) {
796: for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
797: }
798: jj += bss[j];
799: }
800: #if !defined(PETSC_USE_COMPLEX)
801: } else if (addv == MAX_VALUES) {
802: for (j = 0; j < nv; j++) {
803: for (k = 0; k < bss[j]; k++) {
804: for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
805: }
806: jj += bss[j];
807: }
808: #endif
809: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
811: PetscCall(VecRestoreArray(v, &x));
812: for (i = 0; i < nv; i++) PetscCall(VecRestoreArrayRead(s[i], &y[i]));
813: PetscCall(PetscFree2(*(PetscScalar ***)&y, bss));
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: /*@
818: VecStrideGather - Gathers a single component from a multi-component vector into
819: another vector.
821: Collective
823: Input Parameters:
824: + v - the vector
825: . start - starting point of the subvector (defined by a stride)
826: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
828: Output Parameter:
829: . s - the location where the subvector is stored
831: Level: advanced
833: Notes:
834: One must call `VecSetBlockSize()` before this routine to set the stride
835: information, or use a vector created from a multicomponent `DMDA`.
837: If x is the array representing the vector x then this gathers
838: the array (x[start],x[start+stride],x[start+2*stride], ....)
840: The parallel layout of the vector and the subvector must be the same;
841: i.e., nlocal of v = stride*(nlocal of s)
843: Not optimized; could be easily
845: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
846: `VecStrideScatterAll()`
847: @*/
848: PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
849: {
850: PetscFunctionBegin;
853: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
854: 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,
855: v->map->bs);
856: PetscUseTypeMethod(v, stridegather, start, s, addv);
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*@
861: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
863: Collective
865: Input Parameters:
866: + s - the single-component vector
867: . start - starting point of the subvector (defined by a stride)
868: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
870: Output Parameter:
871: . v - the location where the subvector is scattered (the multi-component vector)
873: Level: advanced
875: Notes:
876: One must call `VecSetBlockSize()` on the multi-component vector before this
877: routine to set the stride information, or use a vector created from a multicomponent `DMDA`.
879: The parallel layout of the vector and the subvector must be the same;
880: i.e., nlocal of v = stride*(nlocal of s)
882: Not optimized; could be easily
884: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
885: `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
886: @*/
887: PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
888: {
889: PetscFunctionBegin;
892: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
893: 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,
894: v->map->bs);
895: PetscCall((*v->ops->stridescatter)(s, start, v, addv));
896: PetscFunctionReturn(PETSC_SUCCESS);
897: }
899: /*@
900: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
901: another vector.
903: Collective
905: Input Parameters:
906: + v - the vector
907: . nidx - the number of indices
908: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
909: . 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`
910: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
912: Output Parameter:
913: . s - the location where the subvector is stored
915: Level: advanced
917: Notes:
918: One must call `VecSetBlockSize()` on both vectors before this routine to set the stride
919: information, or use a vector created from a multicomponent `DMDA`.
921: The parallel layout of the vector and the subvector must be the same;
923: Not optimized; could be easily
925: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
926: `VecStrideScatterAll()`
927: @*/
928: PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
929: {
930: PetscFunctionBegin;
933: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
934: PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: /*@
939: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
941: Collective
943: Input Parameters:
944: + s - the smaller-component vector
945: . nidx - the number of indices in idx
946: . 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`
947: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
948: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`
950: Output Parameter:
951: . v - the location where the subvector is into scattered (the multi-component vector)
953: Level: advanced
955: Notes:
956: One must call `VecSetBlockSize()` on the vectors before this
957: routine to set the stride information, or use a vector created from a multicomponent `DMDA`.
959: The parallel layout of the vector and the subvector must be the same;
961: Not optimized; could be easily
963: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
964: `VecStrideScatterAll()`
965: @*/
966: PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
967: {
968: PetscFunctionBegin;
971: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
972: PetscCall((*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
977: {
978: PetscInt i, n, bs, ns;
979: const PetscScalar *x;
980: PetscScalar *y;
982: PetscFunctionBegin;
983: PetscCall(VecGetLocalSize(v, &n));
984: PetscCall(VecGetLocalSize(s, &ns));
985: PetscCall(VecGetArrayRead(v, &x));
986: PetscCall(VecGetArray(s, &y));
988: bs = v->map->bs;
989: 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);
990: x += start;
991: n = n / bs;
993: if (addv == INSERT_VALUES) {
994: for (i = 0; i < n; i++) y[i] = x[bs * i];
995: } else if (addv == ADD_VALUES) {
996: for (i = 0; i < n; i++) y[i] += x[bs * i];
997: #if !defined(PETSC_USE_COMPLEX)
998: } else if (addv == MAX_VALUES) {
999: for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
1000: #endif
1001: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
1003: PetscCall(VecRestoreArrayRead(v, &x));
1004: PetscCall(VecRestoreArray(s, &y));
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
1009: {
1010: PetscInt i, n, bs, ns;
1011: PetscScalar *x;
1012: const PetscScalar *y;
1014: PetscFunctionBegin;
1015: PetscCall(VecGetLocalSize(v, &n));
1016: PetscCall(VecGetLocalSize(s, &ns));
1017: PetscCall(VecGetArray(v, &x));
1018: PetscCall(VecGetArrayRead(s, &y));
1020: bs = v->map->bs;
1021: 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);
1022: x += start;
1023: n = n / bs;
1025: if (addv == INSERT_VALUES) {
1026: for (i = 0; i < n; i++) x[bs * i] = y[i];
1027: } else if (addv == ADD_VALUES) {
1028: for (i = 0; i < n; i++) x[bs * i] += y[i];
1029: #if !defined(PETSC_USE_COMPLEX)
1030: } else if (addv == MAX_VALUES) {
1031: for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
1032: #endif
1033: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
1035: PetscCall(VecRestoreArray(v, &x));
1036: PetscCall(VecRestoreArrayRead(s, &y));
1037: PetscFunctionReturn(PETSC_SUCCESS);
1038: }
1040: PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
1041: {
1042: PetscInt i, j, n, bs, bss, ns;
1043: const PetscScalar *x;
1044: PetscScalar *y;
1046: PetscFunctionBegin;
1047: PetscCall(VecGetLocalSize(v, &n));
1048: PetscCall(VecGetLocalSize(s, &ns));
1049: PetscCall(VecGetArrayRead(v, &x));
1050: PetscCall(VecGetArray(s, &y));
1052: bs = v->map->bs;
1053: bss = s->map->bs;
1054: n = n / bs;
1056: if (PetscDefined(USE_DEBUG)) {
1057: PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1058: for (j = 0; j < nidx; j++) {
1059: PetscCheck(idxv[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxv[j]);
1060: 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);
1061: }
1062: PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not gathering into all locations");
1063: }
1065: if (addv == INSERT_VALUES) {
1066: if (!idxs) {
1067: for (i = 0; i < n; i++) {
1068: for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
1069: }
1070: } else {
1071: for (i = 0; i < n; i++) {
1072: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
1073: }
1074: }
1075: } else if (addv == ADD_VALUES) {
1076: if (!idxs) {
1077: for (i = 0; i < n; i++) {
1078: for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
1079: }
1080: } else {
1081: for (i = 0; i < n; i++) {
1082: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
1083: }
1084: }
1085: #if !defined(PETSC_USE_COMPLEX)
1086: } else if (addv == MAX_VALUES) {
1087: if (!idxs) {
1088: for (i = 0; i < n; i++) {
1089: for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1090: }
1091: } else {
1092: for (i = 0; i < n; i++) {
1093: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1094: }
1095: }
1096: #endif
1097: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
1099: PetscCall(VecRestoreArrayRead(v, &x));
1100: PetscCall(VecRestoreArray(s, &y));
1101: PetscFunctionReturn(PETSC_SUCCESS);
1102: }
1104: PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
1105: {
1106: PetscInt j, i, n, bs, ns, bss;
1107: PetscScalar *x;
1108: const PetscScalar *y;
1110: PetscFunctionBegin;
1111: PetscCall(VecGetLocalSize(v, &n));
1112: PetscCall(VecGetLocalSize(s, &ns));
1113: PetscCall(VecGetArray(v, &x));
1114: PetscCall(VecGetArrayRead(s, &y));
1116: bs = v->map->bs;
1117: bss = s->map->bs;
1118: n = n / bs;
1120: if (PetscDefined(USE_DEBUG)) {
1121: PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1122: for (j = 0; j < bss; j++) {
1123: if (idxs) {
1124: PetscCheck(idxs[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxs[j]);
1125: 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);
1126: }
1127: }
1128: PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not scattering from all locations");
1129: }
1131: if (addv == INSERT_VALUES) {
1132: if (!idxs) {
1133: for (i = 0; i < n; i++) {
1134: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
1135: }
1136: } else {
1137: for (i = 0; i < n; i++) {
1138: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
1139: }
1140: }
1141: } else if (addv == ADD_VALUES) {
1142: if (!idxs) {
1143: for (i = 0; i < n; i++) {
1144: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
1145: }
1146: } else {
1147: for (i = 0; i < n; i++) {
1148: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
1149: }
1150: }
1151: #if !defined(PETSC_USE_COMPLEX)
1152: } else if (addv == MAX_VALUES) {
1153: if (!idxs) {
1154: for (i = 0; i < n; i++) {
1155: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1156: }
1157: } else {
1158: for (i = 0; i < n; i++) {
1159: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1160: }
1161: }
1162: #endif
1163: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
1165: PetscCall(VecRestoreArray(v, &x));
1166: PetscCall(VecRestoreArrayRead(s, &y));
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: static PetscErrorCode VecApplyUnary_Private(Vec v, PetscErrorCode (*unary_op)(Vec), PetscScalar (*UnaryFunc)(PetscScalar))
1171: {
1172: PetscFunctionBegin;
1174: PetscCall(VecSetErrorIfLocked(v, 1));
1175: if (unary_op) {
1177: PetscCall((*unary_op)(v));
1178: } else {
1179: PetscInt n;
1180: PetscScalar *x;
1183: PetscCall(VecGetLocalSize(v, &n));
1184: PetscCall(VecGetArray(v, &x));
1185: for (PetscInt i = 0; i < n; ++i) x[i] = UnaryFunc(x[i]);
1186: PetscCall(VecRestoreArray(v, &x));
1187: }
1188: PetscFunctionReturn(PETSC_SUCCESS);
1189: }
1191: static PetscScalar ScalarReciprocal_Fn(PetscScalar x)
1192: {
1193: const PetscScalar zero = 0.0;
1195: return x == zero ? zero : ((PetscScalar)1.0) / x;
1196: }
1198: PetscErrorCode VecReciprocal_Default(Vec v)
1199: {
1200: PetscFunctionBegin;
1201: PetscCall(VecApplyUnary_Private(v, NULL, ScalarReciprocal_Fn));
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: static PetscScalar ScalarExp_Fn(PetscScalar x)
1206: {
1207: return PetscExpScalar(x);
1208: }
1210: /*@
1211: VecExp - Replaces each component of a vector by e^x_i
1213: Not Collective
1215: Input Parameter:
1216: . v - The vector
1218: Output Parameter:
1219: . v - The vector of exponents
1221: Level: beginner
1223: .seealso: `Vec`, `VecLog()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1225: @*/
1226: PetscErrorCode VecExp(Vec v)
1227: {
1228: PetscFunctionBegin;
1230: PetscCall(VecApplyUnary_Private(v, v->ops->exp, ScalarExp_Fn));
1231: PetscFunctionReturn(PETSC_SUCCESS);
1232: }
1234: static PetscScalar ScalarLog_Fn(PetscScalar x)
1235: {
1236: return PetscLogScalar(x);
1237: }
1239: /*@
1240: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1242: Not Collective
1244: Input Parameter:
1245: . v - The vector
1247: Output Parameter:
1248: . v - The vector of logs
1250: Level: beginner
1252: .seealso: `Vec`, `VecExp()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1254: @*/
1255: PetscErrorCode VecLog(Vec v)
1256: {
1257: PetscFunctionBegin;
1259: PetscCall(VecApplyUnary_Private(v, v->ops->log, ScalarLog_Fn));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: static PetscScalar ScalarAbs_Fn(PetscScalar x)
1264: {
1265: return PetscAbsScalar(x);
1266: }
1268: /*@
1269: VecAbs - Replaces every element in a vector with its absolute value.
1271: Logically Collective
1273: Input Parameter:
1274: . v - the vector
1276: Level: intermediate
1278: @*/
1279: PetscErrorCode VecAbs(Vec v)
1280: {
1281: PetscFunctionBegin;
1283: PetscCall(VecApplyUnary_Private(v, v->ops->abs, ScalarAbs_Fn));
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: static PetscScalar ScalarSqrtAbs_Fn(PetscScalar x)
1288: {
1289: return PetscSqrtScalar(ScalarAbs_Fn(x));
1290: }
1292: /*@
1293: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1295: Not Collective
1297: Input Parameter:
1298: . v - The vector
1300: Level: beginner
1302: Note:
1303: The actual function is sqrt(|x_i|)
1305: .seealso: `Vec`, `VecLog()`, `VecExp()`, `VecReciprocal()`, `VecAbs()`
1307: @*/
1308: PetscErrorCode VecSqrtAbs(Vec v)
1309: {
1310: PetscFunctionBegin;
1312: PetscCall(VecApplyUnary_Private(v, v->ops->sqrt, ScalarSqrtAbs_Fn));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: static PetscScalar ScalarImaginaryPart_Fn(PetscScalar x)
1317: {
1318: const PetscReal imag = PetscImaginaryPart(x);
1320: #if PetscDefined(USE_COMPLEX)
1321: return PetscCMPLX(imag, 0.0);
1322: #else
1323: return imag;
1324: #endif
1325: }
1327: /*@
1328: VecImaginaryPart - Replaces a complex vector with its imginary part
1330: Collective
1332: Input Parameter:
1333: . v - the vector
1335: Level: beginner
1337: .seealso: `Vec`, `VecNorm()`, `VecRealPart()`
1338: @*/
1339: PetscErrorCode VecImaginaryPart(Vec v)
1340: {
1341: PetscFunctionBegin;
1343: PetscCall(VecApplyUnary_Private(v, NULL, ScalarImaginaryPart_Fn));
1344: PetscFunctionReturn(PETSC_SUCCESS);
1345: }
1347: static PetscScalar ScalarRealPart_Fn(PetscScalar x)
1348: {
1349: const PetscReal real = PetscRealPart(x);
1351: #if PetscDefined(USE_COMPLEX)
1352: return PetscCMPLX(real, 0.0);
1353: #else
1354: return real;
1355: #endif
1356: }
1358: /*@
1359: VecRealPart - Replaces a complex vector with its real part
1361: Collective
1363: Input Parameter:
1364: . v - the vector
1366: Level: beginner
1368: .seealso: `Vec`, `VecNorm()`, `VecImaginaryPart()`
1369: @*/
1370: PetscErrorCode VecRealPart(Vec v)
1371: {
1372: PetscFunctionBegin;
1374: PetscCall(VecApplyUnary_Private(v, NULL, ScalarRealPart_Fn));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: /*@
1379: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1381: Collective
1383: Input Parameters:
1384: + s - first vector
1385: - t - second vector
1387: Output Parameters:
1388: + dp - s'conj(t)
1389: - nm - t'conj(t)
1391: Level: advanced
1393: Note:
1394: conj(x) is the complex conjugate of x when x is complex
1396: .seealso: `Vec`, `VecDot()`, `VecNorm()`, `VecDotBegin()`, `VecNormBegin()`, `VecDotEnd()`, `VecNormEnd()`
1398: @*/
1399: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1400: {
1401: PetscScalar work[] = {0.0, 0.0};
1403: PetscFunctionBegin;
1410: PetscCheckSameTypeAndComm(s, 1, t, 2);
1411: PetscCheck(s->map->N == t->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector global lengths");
1412: PetscCheck(s->map->n == t->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector local lengths");
1414: PetscCall(PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0));
1415: if (s->ops->dotnorm2) {
1416: PetscUseTypeMethod(s, dotnorm2, t, work, work + 1);
1417: } else {
1418: const PetscScalar *sx, *tx;
1419: PetscInt n;
1421: PetscCall(VecGetLocalSize(s, &n));
1422: PetscCall(VecGetArrayRead(s, &sx));
1423: PetscCall(VecGetArrayRead(t, &tx));
1424: for (PetscInt i = 0; i < n; ++i) {
1425: const PetscScalar txconj = PetscConj(tx[i]);
1427: work[0] += sx[i] * txconj;
1428: work[1] += tx[i] * txconj;
1429: }
1430: PetscCall(VecRestoreArrayRead(t, &tx));
1431: PetscCall(VecRestoreArrayRead(s, &sx));
1432: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, work, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
1433: PetscCall(PetscLogFlops(4.0 * n));
1434: }
1435: PetscCall(PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0));
1436: *dp = work[0];
1437: *nm = PetscRealPart(work[1]);
1438: PetscFunctionReturn(PETSC_SUCCESS);
1439: }
1441: /*@
1442: VecSum - Computes the sum of all the components of a vector.
1444: Collective
1446: Input Parameter:
1447: . v - the vector
1449: Output Parameter:
1450: . sum - the result
1452: Level: beginner
1454: .seealso: `Vec`, `VecMean()`, `VecNorm()`
1455: @*/
1456: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1457: {
1458: PetscScalar tmp = 0.0;
1460: PetscFunctionBegin;
1463: if (v->ops->sum) {
1464: PetscUseTypeMethod(v, sum, &tmp);
1465: } else {
1466: const PetscScalar *x;
1467: PetscInt n;
1469: PetscCall(VecGetLocalSize(v, &n));
1470: PetscCall(VecGetArrayRead(v, &x));
1471: for (PetscInt i = 0; i < n; ++i) tmp += x[i];
1472: PetscCall(VecRestoreArrayRead(v, &x));
1473: }
1474: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &tmp, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
1475: *sum = tmp;
1476: PetscFunctionReturn(PETSC_SUCCESS);
1477: }
1479: /*@
1480: VecMean - Computes the arithmetic mean of all the components of a vector.
1482: Collective
1484: Input Parameter:
1485: . v - the vector
1487: Output Parameter:
1488: . mean - the result
1490: Level: beginner
1492: .seealso: `Vec`, `VecSum()`, `VecNorm()`
1493: @*/
1494: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1495: {
1496: PetscInt n;
1498: PetscFunctionBegin;
1501: PetscCall(VecGetSize(v, &n));
1502: PetscCall(VecSum(v, mean));
1503: *mean /= n;
1504: PetscFunctionReturn(PETSC_SUCCESS);
1505: }
1507: /*@
1508: VecShift - Shifts all of the components of a vector by computing
1509: `x[i] = x[i] + shift`.
1511: Logically Collective
1513: Input Parameters:
1514: + v - the vector
1515: - shift - the shift
1517: Level: intermediate
1519: .seealso: `Vec`
1520: @*/
1521: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1522: {
1523: PetscFunctionBegin;
1526: PetscCall(VecSetErrorIfLocked(v, 1));
1527: if (shift == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
1529: if (v->ops->shift) {
1530: PetscUseTypeMethod(v, shift, shift);
1531: } else {
1532: PetscInt n;
1533: PetscScalar *x;
1535: PetscCall(VecGetLocalSize(v, &n));
1536: PetscCall(VecGetArray(v, &x));
1537: for (PetscInt i = 0; i < n; ++i) x[i] += shift;
1538: PetscCall(VecRestoreArray(v, &x));
1539: }
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: /*@
1544: VecPermute - Permutes a vector in place using the given ordering.
1546: Input Parameters:
1547: + vec - The vector
1548: . order - The ordering
1549: - inv - The flag for inverting the permutation
1551: Level: beginner
1553: Note:
1554: This function does not yet support parallel Index Sets with non-local permutations
1556: .seealso: `Vec`, `MatPermute()`
1557: @*/
1558: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1559: {
1560: const PetscScalar *array;
1561: PetscScalar *newArray;
1562: const PetscInt *idx;
1563: PetscInt i, rstart, rend;
1565: PetscFunctionBegin;
1568: PetscCall(VecSetErrorIfLocked(x, 1));
1569: PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1570: PetscCall(ISGetIndices(row, &idx));
1571: PetscCall(VecGetArrayRead(x, &array));
1572: PetscCall(PetscMalloc1(x->map->n, &newArray));
1573: if (PetscDefined(USE_DEBUG)) {
1574: 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]);
1575: }
1576: if (!inv) {
1577: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i] - rstart];
1578: } else {
1579: for (i = 0; i < x->map->n; i++) newArray[idx[i] - rstart] = array[i];
1580: }
1581: PetscCall(VecRestoreArrayRead(x, &array));
1582: PetscCall(ISRestoreIndices(row, &idx));
1583: PetscCall(VecReplaceArray(x, newArray));
1584: PetscFunctionReturn(PETSC_SUCCESS);
1585: }
1587: /*@
1588: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1589: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1590: Does NOT take round-off errors into account.
1592: Collective
1594: Input Parameters:
1595: + vec1 - the first vector
1596: - vec2 - the second vector
1598: Output Parameter:
1599: . flg - `PETSC_TRUE` if the vectors are equal; `PETSC_FALSE` otherwise.
1601: Level: intermediate
1603: .seealso: `Vec`
1604: @*/
1605: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1606: {
1607: const PetscScalar *v1, *v2;
1608: PetscInt n1, n2, N1, N2;
1609: PetscBool flg1;
1611: PetscFunctionBegin;
1615: if (vec1 == vec2) *flg = PETSC_TRUE;
1616: else {
1617: PetscCall(VecGetSize(vec1, &N1));
1618: PetscCall(VecGetSize(vec2, &N2));
1619: if (N1 != N2) flg1 = PETSC_FALSE;
1620: else {
1621: PetscCall(VecGetLocalSize(vec1, &n1));
1622: PetscCall(VecGetLocalSize(vec2, &n2));
1623: if (n1 != n2) flg1 = PETSC_FALSE;
1624: else {
1625: PetscCall(VecGetArrayRead(vec1, &v1));
1626: PetscCall(VecGetArrayRead(vec2, &v2));
1627: PetscCall(PetscArraycmp(v1, v2, n1, &flg1));
1628: PetscCall(VecRestoreArrayRead(vec1, &v1));
1629: PetscCall(VecRestoreArrayRead(vec2, &v2));
1630: }
1631: }
1632: /* combine results from all processors */
1633: PetscCall(MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1)));
1634: }
1635: PetscFunctionReturn(PETSC_SUCCESS);
1636: }
1638: /*@
1639: VecUniqueEntries - Compute the number of unique entries, and those entries
1641: Collective
1643: Input Parameter:
1644: . vec - the vector
1646: Output Parameters:
1647: + n - The number of unique entries
1648: - e - The entries
1650: Level: intermediate
1652: .seealso: `Vec`
1653: @*/
1654: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1655: {
1656: const PetscScalar *v;
1657: PetscScalar *tmp, *vals;
1658: PetscMPIInt *N, *displs, l;
1659: PetscInt ng, m, i, j, p;
1660: PetscMPIInt size;
1662: PetscFunctionBegin;
1665: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1666: PetscCall(VecGetLocalSize(vec, &m));
1667: PetscCall(VecGetArrayRead(vec, &v));
1668: PetscCall(PetscMalloc2(m, &tmp, size, &N));
1669: for (i = 0, j = 0, l = 0; i < m; ++i) {
1670: /* Can speed this up with sorting */
1671: for (j = 0; j < l; ++j) {
1672: if (v[i] == tmp[j]) break;
1673: }
1674: if (j == l) {
1675: tmp[j] = v[i];
1676: ++l;
1677: }
1678: }
1679: PetscCall(VecRestoreArrayRead(vec, &v));
1680: /* Gather serial results */
1681: PetscCallMPI(MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec)));
1682: for (p = 0, ng = 0; p < size; ++p) ng += N[p];
1683: PetscCall(PetscMalloc2(ng, &vals, size + 1, &displs));
1684: for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
1685: PetscCallMPI(MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec)));
1686: /* Find unique entries */
1687: #ifdef PETSC_USE_COMPLEX
1688: SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1689: #else
1690: *n = displs[size];
1691: PetscCall(PetscSortRemoveDupsReal(n, (PetscReal *)vals));
1692: if (e) {
1694: PetscCall(PetscMalloc1(*n, e));
1695: for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
1696: }
1697: PetscCall(PetscFree2(vals, displs));
1698: PetscCall(PetscFree2(tmp, N));
1699: PetscFunctionReturn(PETSC_SUCCESS);
1700: #endif
1701: }