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: Notes:
19: One must call VecSetBlockSize() before this routine to set the stride
20: information, or use a vector created from a multicomponent DMDA.
22: This will only work if the desire subvector is a stride subvector
24: Level: advanced
26: .seealso: `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;
34: VecGetLocalSize(v, &n);
35: VecGetBlockSize(v, &bs);
38: VecGetArray(v, &x);
39: for (i = start; i < n; i += bs) x[i] = s;
40: VecRestoreArray(v, &x);
41: return 0;
42: }
44: /*@
45: VecStrideScale - Scales a subvector of a vector defined
46: by a starting point and a stride.
48: Logically Collective
50: Input Parameters:
51: + v - the vector
52: . start - starting point of the subvector (defined by a stride)
53: - scale - value to multiply each subvector entry by
55: Notes:
56: One must call VecSetBlockSize() before this routine to set the stride
57: information, or use a vector created from a multicomponent DMDA.
59: This will only work if the desire subvector is a stride subvector
61: Level: advanced
63: .seealso: `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideScale()`
64: @*/
65: PetscErrorCode VecStrideScale(Vec v, PetscInt start, PetscScalar scale)
66: {
67: PetscInt i, n, bs;
68: PetscScalar *x;
71: VecGetLocalSize(v, &n);
72: VecGetBlockSize(v, &bs);
75: VecGetArray(v, &x);
76: for (i = start; i < n; i += bs) x[i] *= scale;
77: VecRestoreArray(v, &x);
78: return 0;
79: }
81: /*@
82: VecStrideNorm - Computes the norm of subvector of a vector defined
83: by a starting point and a stride.
85: Collective
87: Input Parameters:
88: + v - the vector
89: . start - starting point of the subvector (defined by a stride)
90: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
92: Output Parameter:
93: . norm - the norm
95: Notes:
96: One must call VecSetBlockSize() before this routine to set the stride
97: information, or use a vector created from a multicomponent DMDA.
99: If x is the array representing the vector x then this computes the norm
100: of the array (x[start],x[start+stride],x[start+2*stride], ....)
102: This is useful for computing, say the norm of the pressure variable when
103: the pressure is stored (interlaced) with other variables, say density etc.
105: This will only work if the desire subvector is a stride subvector
107: Level: advanced
109: .seealso: `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
110: @*/
111: PetscErrorCode VecStrideNorm(Vec v, PetscInt start, NormType ntype, PetscReal *nrm)
112: {
113: PetscInt i, n, bs;
114: const PetscScalar *x;
115: PetscReal tnorm;
120: VecGetLocalSize(v, &n);
121: VecGetBlockSize(v, &bs);
124: VecGetArrayRead(v, &x);
125: if (ntype == NORM_2) {
126: PetscScalar sum = 0.0;
127: for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
128: tnorm = PetscRealPart(sum);
129: MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v));
130: *nrm = PetscSqrtReal(*nrm);
131: } else if (ntype == NORM_1) {
132: tnorm = 0.0;
133: for (i = start; i < n; i += bs) tnorm += PetscAbsScalar(x[i]);
134: MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v));
135: } else if (ntype == NORM_INFINITY) {
136: tnorm = 0.0;
137: for (i = start; i < n; i += bs) {
138: if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
139: }
140: MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v));
141: } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
142: VecRestoreArrayRead(v, &x);
143: return 0;
144: }
146: /*@
147: VecStrideMax - Computes the maximum of subvector of a vector defined
148: by a starting point and a stride and optionally its location.
150: Collective
152: Input Parameters:
153: + v - the vector
154: - start - starting point of the subvector (defined by a stride)
156: Output Parameters:
157: + idex - the location where the maximum occurred (pass NULL if not required)
158: - nrm - the maximum value in the subvector
160: Notes:
161: One must call VecSetBlockSize() before this routine to set the stride
162: information, or use a vector created from a multicomponent DMDA.
164: If xa is the array representing the vector x, then this computes the max
165: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
167: This is useful for computing, say the maximum of the pressure variable when
168: the pressure is stored (interlaced) with other variables, e.g., density, etc.
169: This will only work if the desire subvector is a stride subvector.
171: Level: advanced
173: .seealso: `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
174: @*/
175: PetscErrorCode VecStrideMax(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
176: {
177: PetscInt i, n, bs, id = -1;
178: const PetscScalar *x;
179: PetscReal max = PETSC_MIN_REAL;
183: VecGetLocalSize(v, &n);
184: VecGetBlockSize(v, &bs);
187: VecGetArrayRead(v, &x);
188: for (i = start; i < n; i += bs) {
189: if (PetscRealPart(x[i]) > max) {
190: max = PetscRealPart(x[i]);
191: id = i;
192: }
193: }
194: VecRestoreArrayRead(v, &x);
195: #if defined(PETSC_HAVE_MPIUNI)
196: *nrm = max;
197: if (idex) *idex = id;
198: #else
199: if (!idex) {
200: MPIU_Allreduce(&max, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v));
201: } else {
202: struct {
203: PetscReal v;
204: PetscInt i;
205: } in, out;
206: PetscInt rstart;
208: VecGetOwnershipRange(v, &rstart, NULL);
209: in.v = max;
210: in.i = rstart + id;
211: MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v));
212: *nrm = out.v;
213: *idex = out.i;
214: }
215: #endif
216: return 0;
217: }
219: /*@
220: VecStrideMin - Computes the minimum of subvector of a vector defined
221: by a starting point and a stride and optionally its location.
223: Collective
225: Input Parameters:
226: + v - the vector
227: - start - starting point of the subvector (defined by a stride)
229: Output Parameters:
230: + idex - the location where the minimum occurred. (pass NULL if not required)
231: - nrm - the minimum value in the subvector
233: Level: advanced
235: Notes:
236: One must call VecSetBlockSize() before this routine to set the stride
237: information, or use a vector created from a multicomponent DMDA.
239: If xa is the array representing the vector x, then this computes the min
240: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
242: This is useful for computing, say the minimum of the pressure variable when
243: the pressure is stored (interlaced) with other variables, e.g., density, etc.
244: This will only work if the desire subvector is a stride subvector.
246: .seealso: `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
247: @*/
248: PetscErrorCode VecStrideMin(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
249: {
250: PetscInt i, n, bs, id = -1;
251: const PetscScalar *x;
252: PetscReal min = PETSC_MAX_REAL;
256: VecGetLocalSize(v, &n);
257: VecGetBlockSize(v, &bs);
260: VecGetArrayRead(v, &x);
261: for (i = start; i < n; i += bs) {
262: if (PetscRealPart(x[i]) < min) {
263: min = PetscRealPart(x[i]);
264: id = i;
265: }
266: }
267: VecRestoreArrayRead(v, &x);
268: #if defined(PETSC_HAVE_MPIUNI)
269: *nrm = min;
270: if (idex) *idex = id;
271: #else
272: if (!idex) {
273: MPIU_Allreduce(&min, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v));
274: } else {
275: struct {
276: PetscReal v;
277: PetscInt i;
278: } in, out;
279: PetscInt rstart;
281: VecGetOwnershipRange(v, &rstart, NULL);
282: in.v = min;
283: in.i = rstart + id;
284: MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v));
285: *nrm = out.v;
286: *idex = out.i;
287: }
288: #endif
289: return 0;
290: }
292: /*@
293: VecStrideSum - Computes the sum of subvector of a vector defined
294: by a starting point and a stride.
296: Collective
298: Input Parameters:
299: + v - the vector
300: - start - starting point of the subvector (defined by a stride)
302: Output Parameter:
303: . sum - the sum
305: Notes:
306: One must call `VecSetBlockSize()` before this routine to set the stride
307: information, or use a vector created from a multicomponent `DMDA`.
309: If x is the array representing the vector x then this computes the sum
310: of the array (x[start],x[start+stride],x[start+2*stride], ....)
312: Level: advanced
314: .seealso: `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
315: @*/
316: PetscErrorCode VecStrideSum(Vec v, PetscInt start, PetscScalar *sum)
317: {
318: PetscInt i, n, bs;
319: const PetscScalar *x;
320: PetscScalar local_sum = 0.0;
324: VecGetLocalSize(v, &n);
325: VecGetBlockSize(v, &bs);
328: VecGetArrayRead(v, &x);
329: for (i = start; i < n; i += bs) local_sum += x[i];
330: MPIU_Allreduce(&local_sum, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v));
331: VecRestoreArrayRead(v, &x);
332: return 0;
333: }
335: /*@
336: VecStrideScaleAll - Scales the subvectors of a vector defined
337: by a starting point and a stride.
339: Logically Collective
341: Input Parameters:
342: + v - the vector
343: - scales - values to multiply each subvector entry by
345: Notes:
346: One must call VecSetBlockSize() before this routine to set the stride
347: information, or use a vector created from a multicomponent DMDA.
349: The dimension of scales must be the same as the vector block size
351: Level: advanced
353: .seealso: `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
354: @*/
355: PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
356: {
357: PetscInt i, j, n, bs;
358: PetscScalar *x;
362: VecGetLocalSize(v, &n);
363: VecGetBlockSize(v, &bs);
364: VecGetArray(v, &x);
365: /* need to provide optimized code for each bs */
366: for (i = 0; i < n; i += bs) {
367: for (j = 0; j < bs; j++) x[i + j] *= scales[j];
368: }
369: VecRestoreArray(v, &x);
370: return 0;
371: }
373: /*@
374: VecStrideNormAll - Computes the norms of subvectors of a vector defined
375: by a starting point and a stride.
377: Collective
379: Input Parameters:
380: + v - the vector
381: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
383: Output Parameter:
384: . nrm - the norms
386: Notes:
387: One must call VecSetBlockSize() before this routine to set the stride
388: information, or use a vector created from a multicomponent DMDA.
390: If x is the array representing the vector x then this computes the norm
391: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
393: The dimension of nrm must be the same as the vector block size
395: This will only work if the desire subvector is a stride subvector
397: Level: advanced
399: .seealso: `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
400: @*/
401: PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
402: {
403: PetscInt i, j, n, bs;
404: const PetscScalar *x;
405: PetscReal tnorm[128];
406: MPI_Comm comm;
411: VecGetLocalSize(v, &n);
412: VecGetArrayRead(v, &x);
413: PetscObjectGetComm((PetscObject)v, &comm);
415: VecGetBlockSize(v, &bs);
418: if (ntype == NORM_2) {
419: PetscScalar sum[128];
420: for (j = 0; j < bs; j++) sum[j] = 0.0;
421: for (i = 0; i < n; i += bs) {
422: for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
423: }
424: for (j = 0; j < bs; j++) tnorm[j] = PetscRealPart(sum[j]);
426: MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm);
427: for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
428: } else if (ntype == NORM_1) {
429: for (j = 0; j < bs; j++) tnorm[j] = 0.0;
431: for (i = 0; i < n; i += bs) {
432: for (j = 0; j < bs; j++) tnorm[j] += PetscAbsScalar(x[i + j]);
433: }
435: MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm);
436: } else if (ntype == NORM_INFINITY) {
437: PetscReal tmp;
438: for (j = 0; j < bs; j++) tnorm[j] = 0.0;
440: for (i = 0; i < n; i += bs) {
441: for (j = 0; j < bs; j++) {
442: if ((tmp = PetscAbsScalar(x[i + j])) > tnorm[j]) tnorm[j] = tmp;
443: /* check special case of tmp == NaN */
444: if (tmp != tmp) {
445: tnorm[j] = tmp;
446: break;
447: }
448: }
449: }
450: MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_MAX, comm);
451: } else SETERRQ(comm, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
452: VecRestoreArrayRead(v, &x);
453: return 0;
454: }
456: /*@
457: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
458: by a starting point and a stride and optionally its location.
460: Collective
462: Input Parameter:
463: . v - the vector
465: Output Parameters:
466: + index - the location where the maximum occurred (not supported, pass NULL,
467: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
468: - nrm - the maximum values of each subvector
470: Notes:
471: One must call VecSetBlockSize() before this routine to set the stride
472: information, or use a vector created from a multicomponent DMDA.
474: The dimension of nrm must be the same as the vector block size
476: Level: advanced
478: .seealso: `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
479: @*/
480: PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
481: {
482: PetscInt i, j, n, bs;
483: const PetscScalar *x;
484: PetscReal max[128], tmp;
485: MPI_Comm comm;
490: VecGetLocalSize(v, &n);
491: VecGetArrayRead(v, &x);
492: PetscObjectGetComm((PetscObject)v, &comm);
494: VecGetBlockSize(v, &bs);
497: if (!n) {
498: for (j = 0; j < bs; j++) max[j] = PETSC_MIN_REAL;
499: } else {
500: for (j = 0; j < bs; j++) max[j] = PetscRealPart(x[j]);
502: for (i = bs; i < n; i += bs) {
503: for (j = 0; j < bs; j++) {
504: if ((tmp = PetscRealPart(x[i + j])) > max[j]) max[j] = tmp;
505: }
506: }
507: }
508: MPIU_Allreduce(max, nrm, bs, MPIU_REAL, MPIU_MAX, comm);
510: VecRestoreArrayRead(v, &x);
511: return 0;
512: }
514: /*@
515: VecStrideMinAll - Computes the minimum of subvector of a vector defined
516: by a starting point and a stride and optionally its location.
518: Collective
520: Input Parameter:
521: . v - the vector
523: Output Parameters:
524: + idex - the location where the minimum occurred (not supported, pass NULL,
525: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
526: - nrm - the minimums of each subvector
528: Level: advanced
530: Notes:
531: One must call VecSetBlockSize() before this routine to set the stride
532: information, or use a vector created from a multicomponent DMDA.
534: The dimension of nrm must be the same as the vector block size
536: .seealso: `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
537: @*/
538: PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
539: {
540: PetscInt i, n, bs, j;
541: const PetscScalar *x;
542: PetscReal min[128], tmp;
543: MPI_Comm comm;
548: VecGetLocalSize(v, &n);
549: VecGetArrayRead(v, &x);
550: PetscObjectGetComm((PetscObject)v, &comm);
552: VecGetBlockSize(v, &bs);
555: if (!n) {
556: for (j = 0; j < bs; j++) min[j] = PETSC_MAX_REAL;
557: } else {
558: for (j = 0; j < bs; j++) min[j] = PetscRealPart(x[j]);
560: for (i = bs; i < n; i += bs) {
561: for (j = 0; j < bs; j++) {
562: if ((tmp = PetscRealPart(x[i + j])) < min[j]) min[j] = tmp;
563: }
564: }
565: }
566: MPIU_Allreduce(min, nrm, bs, MPIU_REAL, MPIU_MIN, comm);
568: VecRestoreArrayRead(v, &x);
569: return 0;
570: }
572: /*@
573: VecStrideSumAll - Computes the sums of subvectors of a vector defined by a stride.
575: Collective
577: Input Parameter:
578: . v - the vector
580: Output Parameter:
581: . sums - the sums
583: Notes:
584: One must call `VecSetBlockSize()` before this routine to set the stride
585: information, or use a vector created from a multicomponent `DMDA`.
587: If x is the array representing the vector x then this computes the sum
588: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
590: Level: advanced
592: .seealso: `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
593: @*/
594: PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
595: {
596: PetscInt i, j, n, bs;
597: const PetscScalar *x;
598: PetscScalar local_sums[128];
599: MPI_Comm comm;
603: VecGetLocalSize(v, &n);
604: VecGetArrayRead(v, &x);
605: PetscObjectGetComm((PetscObject)v, &comm);
607: VecGetBlockSize(v, &bs);
610: for (j = 0; j < bs; j++) local_sums[j] = 0.0;
611: for (i = 0; i < n; i += bs) {
612: for (j = 0; j < bs; j++) local_sums[j] += x[i + j];
613: }
614: MPIU_Allreduce(local_sums, sums, bs, MPIU_SCALAR, MPIU_SUM, comm);
616: VecRestoreArrayRead(v, &x);
617: return 0;
618: }
620: /*----------------------------------------------------------------------------------------------*/
621: /*@
622: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
623: separate vectors.
625: Collective
627: Input Parameters:
628: + v - the vector
629: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
631: Output Parameter:
632: . s - the location where the subvectors are stored
634: Notes:
635: One must call VecSetBlockSize() before this routine to set the stride
636: information, or use a vector created from a multicomponent DMDA.
638: If x is the array representing the vector x then this gathers
639: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
640: for start=0,1,2,...bs-1
642: The parallel layout of the vector and the subvector must be the same;
643: i.e., nlocal of v = stride*(nlocal of s)
645: Not optimized; could be easily
647: Level: advanced
649: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
650: `VecStrideScatterAll()`
651: @*/
652: PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
653: {
654: PetscInt i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
655: PetscScalar **y;
656: const PetscScalar *x;
661: VecGetLocalSize(v, &n);
662: VecGetLocalSize(s[0], &n2);
663: VecGetArrayRead(v, &x);
664: VecGetBlockSize(v, &bs);
667: PetscMalloc2(bs, &y, bs, &bss);
668: nv = 0;
669: nvc = 0;
670: for (i = 0; i < bs; i++) {
671: VecGetBlockSize(s[i], &bss[i]);
672: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
673: VecGetArray(s[i], &y[i]);
674: nvc += bss[i];
675: nv++;
677: if (nvc == bs) break;
678: }
680: n = n / bs;
682: jj = 0;
683: if (addv == INSERT_VALUES) {
684: for (j = 0; j < nv; j++) {
685: for (k = 0; k < bss[j]; k++) {
686: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
687: }
688: jj += bss[j];
689: }
690: } else if (addv == ADD_VALUES) {
691: for (j = 0; j < nv; j++) {
692: for (k = 0; k < bss[j]; k++) {
693: for (i = 0; i < n; i++) y[j][i * bss[j] + k] += x[bs * i + jj + k];
694: }
695: jj += bss[j];
696: }
697: #if !defined(PETSC_USE_COMPLEX)
698: } else if (addv == MAX_VALUES) {
699: for (j = 0; j < nv; j++) {
700: for (k = 0; k < bss[j]; k++) {
701: for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
702: }
703: jj += bss[j];
704: }
705: #endif
706: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
708: VecRestoreArrayRead(v, &x);
709: for (i = 0; i < nv; i++) VecRestoreArray(s[i], &y[i]);
711: PetscFree2(y, bss);
712: return 0;
713: }
715: /*@
716: VecStrideScatterAll - Scatters all the single components from separate vectors into
717: a multi-component vector.
719: Collective
721: Input Parameters:
722: + s - the location where the subvectors are stored
723: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
725: Output Parameter:
726: . v - the multicomponent vector
728: Notes:
729: One must call VecSetBlockSize() before this routine to set the stride
730: information, or use a vector created from a multicomponent DMDA.
732: The parallel layout of the vector and the subvector must be the same;
733: i.e., nlocal of v = stride*(nlocal of s)
735: Not optimized; could be easily
737: Level: advanced
739: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
740: `VecStrideScatterAll()`
741: @*/
742: PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
743: {
744: PetscInt i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
745: PetscScalar *x;
746: PetscScalar const **y;
751: VecGetLocalSize(v, &n);
752: VecGetLocalSize(s[0], &n2);
753: VecGetArray(v, &x);
754: VecGetBlockSize(v, &bs);
757: PetscMalloc2(bs, (PetscScalar ***)&y, bs, &bss);
758: nv = 0;
759: nvc = 0;
760: for (i = 0; i < bs; i++) {
761: VecGetBlockSize(s[i], &bss[i]);
762: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
763: VecGetArrayRead(s[i], &y[i]);
764: nvc += bss[i];
765: nv++;
767: if (nvc == bs) break;
768: }
770: n = n / bs;
772: jj = 0;
773: if (addv == INSERT_VALUES) {
774: for (j = 0; j < nv; j++) {
775: for (k = 0; k < bss[j]; k++) {
776: for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
777: }
778: jj += bss[j];
779: }
780: } else if (addv == ADD_VALUES) {
781: for (j = 0; j < nv; j++) {
782: for (k = 0; k < bss[j]; k++) {
783: for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
784: }
785: jj += bss[j];
786: }
787: #if !defined(PETSC_USE_COMPLEX)
788: } else if (addv == MAX_VALUES) {
789: for (j = 0; j < nv; j++) {
790: for (k = 0; k < bss[j]; k++) {
791: for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
792: }
793: jj += bss[j];
794: }
795: #endif
796: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
798: VecRestoreArray(v, &x);
799: for (i = 0; i < nv; i++) VecRestoreArrayRead(s[i], &y[i]);
800: PetscFree2(*(PetscScalar ***)&y, bss);
801: return 0;
802: }
804: /*@
805: VecStrideGather - Gathers a single component from a multi-component vector into
806: another vector.
808: Collective
810: Input Parameters:
811: + v - the vector
812: . start - starting point of the subvector (defined by a stride)
813: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
815: Output Parameter:
816: . s - the location where the subvector is stored
818: Notes:
819: One must call VecSetBlockSize() before this routine to set the stride
820: information, or use a vector created from a multicomponent DMDA.
822: If x is the array representing the vector x then this gathers
823: the array (x[start],x[start+stride],x[start+2*stride], ....)
825: The parallel layout of the vector and the subvector must be the same;
826: i.e., nlocal of v = stride*(nlocal of s)
828: Not optimized; could be easily
830: Level: advanced
832: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
833: `VecStrideScatterAll()`
834: @*/
835: PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
836: {
841: v->map->bs);
842: PetscUseTypeMethod(v, stridegather, start, s, addv);
843: return 0;
844: }
846: /*@
847: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
849: Collective
851: Input Parameters:
852: + s - the single-component vector
853: . start - starting point of the subvector (defined by a stride)
854: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
856: Output Parameter:
857: . v - the location where the subvector is scattered (the multi-component vector)
859: Notes:
860: One must call VecSetBlockSize() on the multi-component vector before this
861: routine to set the stride information, or use a vector created from a multicomponent DMDA.
863: The parallel layout of the vector and the subvector must be the same;
864: i.e., nlocal of v = stride*(nlocal of s)
866: Not optimized; could be easily
868: Level: advanced
870: .seealso: `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
871: `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
872: @*/
873: PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
874: {
879: v->map->bs);
880: (*v->ops->stridescatter)(s, start, v, addv);
881: return 0;
882: }
884: /*@
885: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
886: another vector.
888: Collective
890: Input Parameters:
891: + v - the vector
892: . nidx - the number of indices
893: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
894: . 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
895: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
897: Output Parameter:
898: . s - the location where the subvector is stored
900: Notes:
901: One must call VecSetBlockSize() on both vectors before this routine to set the stride
902: information, or use a vector created from a multicomponent DMDA.
904: The parallel layout of the vector and the subvector must be the same;
906: Not optimized; could be easily
908: Level: advanced
910: .seealso: `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
911: `VecStrideScatterAll()`
912: @*/
913: PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
914: {
917: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
918: PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
919: return 0;
920: }
922: /*@
923: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
925: Collective
927: Input Parameters:
928: + s - the smaller-component vector
929: . nidx - the number of indices in idx
930: . 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
931: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
932: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
934: Output Parameter:
935: . v - the location where the subvector is into scattered (the multi-component vector)
937: Notes:
938: One must call VecSetBlockSize() on the vectors before this
939: routine to set the stride information, or use a vector created from a multicomponent DMDA.
941: The parallel layout of the vector and the subvector must be the same;
943: Not optimized; could be easily
945: Level: advanced
947: .seealso: `VecStrideNorm()`, `VecStrideGather()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
948: `VecStrideScatterAll()`
949: @*/
950: PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
951: {
954: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
955: (*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv);
956: return 0;
957: }
959: PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
960: {
961: PetscInt i, n, bs, ns;
962: const PetscScalar *x;
963: PetscScalar *y;
965: VecGetLocalSize(v, &n);
966: VecGetLocalSize(s, &ns);
967: VecGetArrayRead(v, &x);
968: VecGetArray(s, &y);
970: bs = v->map->bs;
972: x += start;
973: n = n / bs;
975: if (addv == INSERT_VALUES) {
976: for (i = 0; i < n; i++) y[i] = x[bs * i];
977: } else if (addv == ADD_VALUES) {
978: for (i = 0; i < n; i++) y[i] += x[bs * i];
979: #if !defined(PETSC_USE_COMPLEX)
980: } else if (addv == MAX_VALUES) {
981: for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
982: #endif
983: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
985: VecRestoreArrayRead(v, &x);
986: VecRestoreArray(s, &y);
987: return 0;
988: }
990: PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
991: {
992: PetscInt i, n, bs, ns;
993: PetscScalar *x;
994: const PetscScalar *y;
996: VecGetLocalSize(v, &n);
997: VecGetLocalSize(s, &ns);
998: VecGetArray(v, &x);
999: VecGetArrayRead(s, &y);
1001: bs = v->map->bs;
1003: x += start;
1004: n = n / bs;
1006: if (addv == INSERT_VALUES) {
1007: for (i = 0; i < n; i++) x[bs * i] = y[i];
1008: } else if (addv == ADD_VALUES) {
1009: for (i = 0; i < n; i++) x[bs * i] += y[i];
1010: #if !defined(PETSC_USE_COMPLEX)
1011: } else if (addv == MAX_VALUES) {
1012: for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
1013: #endif
1014: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
1016: VecRestoreArray(v, &x);
1017: VecRestoreArrayRead(s, &y);
1018: return 0;
1019: }
1021: PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
1022: {
1023: PetscInt i, j, n, bs, bss, ns;
1024: const PetscScalar *x;
1025: PetscScalar *y;
1027: VecGetLocalSize(v, &n);
1028: VecGetLocalSize(s, &ns);
1029: VecGetArrayRead(v, &x);
1030: VecGetArray(s, &y);
1032: bs = v->map->bs;
1033: bss = s->map->bs;
1034: n = n / bs;
1036: if (PetscDefined(USE_DEBUG)) {
1038: for (j = 0; j < nidx; j++) {
1041: }
1043: }
1045: if (addv == INSERT_VALUES) {
1046: if (!idxs) {
1047: for (i = 0; i < n; i++) {
1048: for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
1049: }
1050: } else {
1051: for (i = 0; i < n; i++) {
1052: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
1053: }
1054: }
1055: } else if (addv == ADD_VALUES) {
1056: if (!idxs) {
1057: for (i = 0; i < n; i++) {
1058: for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
1059: }
1060: } else {
1061: for (i = 0; i < n; i++) {
1062: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
1063: }
1064: }
1065: #if !defined(PETSC_USE_COMPLEX)
1066: } else if (addv == MAX_VALUES) {
1067: if (!idxs) {
1068: for (i = 0; i < n; i++) {
1069: for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1070: }
1071: } else {
1072: for (i = 0; i < n; i++) {
1073: for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1074: }
1075: }
1076: #endif
1077: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
1079: VecRestoreArrayRead(v, &x);
1080: VecRestoreArray(s, &y);
1081: return 0;
1082: }
1084: PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
1085: {
1086: PetscInt j, i, n, bs, ns, bss;
1087: PetscScalar *x;
1088: const PetscScalar *y;
1090: VecGetLocalSize(v, &n);
1091: VecGetLocalSize(s, &ns);
1092: VecGetArray(v, &x);
1093: VecGetArrayRead(s, &y);
1095: bs = v->map->bs;
1096: bss = s->map->bs;
1097: n = n / bs;
1099: if (PetscDefined(USE_DEBUG)) {
1101: for (j = 0; j < bss; j++) {
1102: if (idxs) {
1105: }
1106: }
1108: }
1110: if (addv == INSERT_VALUES) {
1111: if (!idxs) {
1112: for (i = 0; i < n; i++) {
1113: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
1114: }
1115: } else {
1116: for (i = 0; i < n; i++) {
1117: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
1118: }
1119: }
1120: } else if (addv == ADD_VALUES) {
1121: if (!idxs) {
1122: for (i = 0; i < n; i++) {
1123: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
1124: }
1125: } else {
1126: for (i = 0; i < n; i++) {
1127: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
1128: }
1129: }
1130: #if !defined(PETSC_USE_COMPLEX)
1131: } else if (addv == MAX_VALUES) {
1132: if (!idxs) {
1133: for (i = 0; i < n; i++) {
1134: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1135: }
1136: } else {
1137: for (i = 0; i < n; i++) {
1138: for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1139: }
1140: }
1141: #endif
1142: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");
1144: VecRestoreArray(v, &x);
1145: VecRestoreArrayRead(s, &y);
1146: return 0;
1147: }
1149: PetscErrorCode VecReciprocal_Default(Vec v)
1150: {
1151: PetscInt i, n;
1152: PetscScalar *x;
1154: VecGetLocalSize(v, &n);
1155: VecGetArray(v, &x);
1156: for (i = 0; i < n; i++) {
1157: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0 / x[i];
1158: }
1159: VecRestoreArray(v, &x);
1160: return 0;
1161: }
1163: /*@
1164: VecExp - Replaces each component of a vector by e^x_i
1166: Not collective
1168: Input Parameter:
1169: . v - The vector
1171: Output Parameter:
1172: . v - The vector of exponents
1174: Level: beginner
1176: .seealso: `VecLog()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1178: @*/
1179: PetscErrorCode VecExp(Vec v)
1180: {
1181: PetscScalar *x;
1182: PetscInt i, n;
1185: if (v->ops->exp) PetscUseTypeMethod(v, exp);
1186: else {
1187: VecGetLocalSize(v, &n);
1188: VecGetArray(v, &x);
1189: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1190: VecRestoreArray(v, &x);
1191: }
1192: return 0;
1193: }
1195: /*@
1196: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1198: Not collective
1200: Input Parameter:
1201: . v - The vector
1203: Output Parameter:
1204: . v - The vector of logs
1206: Level: beginner
1208: .seealso: `VecExp()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`
1210: @*/
1211: PetscErrorCode VecLog(Vec v)
1212: {
1213: PetscScalar *x;
1214: PetscInt i, n;
1217: if (v->ops->log) PetscUseTypeMethod(v, log);
1218: else {
1219: VecGetLocalSize(v, &n);
1220: VecGetArray(v, &x);
1221: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1222: VecRestoreArray(v, &x);
1223: }
1224: return 0;
1225: }
1227: /*@
1228: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1230: Not collective
1232: Input Parameter:
1233: . v - The vector
1235: Output Parameter:
1236: . v - The vector square root
1238: Level: beginner
1240: Note: The actual function is sqrt(|x_i|)
1242: .seealso: `VecLog()`, `VecExp()`, `VecReciprocal()`, `VecAbs()`
1244: @*/
1245: PetscErrorCode VecSqrtAbs(Vec v)
1246: {
1247: PetscScalar *x;
1248: PetscInt i, n;
1251: if (v->ops->sqrt) PetscUseTypeMethod(v, sqrt);
1252: else {
1253: VecGetLocalSize(v, &n);
1254: VecGetArray(v, &x);
1255: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1256: VecRestoreArray(v, &x);
1257: }
1258: return 0;
1259: }
1261: /*@
1262: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1264: Collective
1266: Input Parameters:
1267: + s - first vector
1268: - t - second vector
1270: Output Parameters:
1271: + dp - s'conj(t)
1272: - nm - t'conj(t)
1274: Level: advanced
1276: Notes:
1277: conj(x) is the complex conjugate of x when x is complex
1279: .seealso: `VecDot()`, `VecNorm()`, `VecDotBegin()`, `VecNormBegin()`, `VecDotEnd()`, `VecNormEnd()`
1281: @*/
1282: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1283: {
1284: const PetscScalar *sx, *tx;
1285: PetscScalar dpx = 0.0, nmx = 0.0, work[2], sum[2];
1286: PetscInt i, n;
1298: PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0);
1299: if (s->ops->dotnorm2) {
1300: PetscUseTypeMethod(s, dotnorm2, t, dp, &dpx);
1301: *nm = PetscRealPart(dpx);
1302: } else {
1303: VecGetLocalSize(s, &n);
1304: VecGetArrayRead(s, &sx);
1305: VecGetArrayRead(t, &tx);
1307: for (i = 0; i < n; i++) {
1308: dpx += sx[i] * PetscConj(tx[i]);
1309: nmx += tx[i] * PetscConj(tx[i]);
1310: }
1311: work[0] = dpx;
1312: work[1] = nmx;
1314: MPIU_Allreduce(work, sum, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s));
1315: *dp = sum[0];
1316: *nm = PetscRealPart(sum[1]);
1318: VecRestoreArrayRead(t, &tx);
1319: VecRestoreArrayRead(s, &sx);
1320: PetscLogFlops(4.0 * n);
1321: }
1322: PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0);
1323: return 0;
1324: }
1326: /*@
1327: VecSum - Computes the sum of all the components of a vector.
1329: Collective
1331: Input Parameter:
1332: . v - the vector
1334: Output Parameter:
1335: . sum - the result
1337: Level: beginner
1339: .seealso: `VecMean()`, `VecNorm()`
1340: @*/
1341: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1342: {
1343: PetscInt i, n;
1344: const PetscScalar *x;
1348: *sum = 0.0;
1349: if (v->ops->sum) {
1350: PetscUseTypeMethod(v, sum, sum);
1351: } else {
1352: VecGetLocalSize(v, &n);
1353: VecGetArrayRead(v, &x);
1354: for (i = 0; i < n; i++) *sum += x[i];
1355: VecRestoreArrayRead(v, &x);
1356: }
1357: MPIU_Allreduce(MPI_IN_PLACE, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v));
1358: return 0;
1359: }
1361: /*@
1362: VecMean - Computes the arithmetic mean of all the components of a vector.
1364: Collective
1366: Input Parameter:
1367: . v - the vector
1369: Output Parameter:
1370: . mean - the result
1372: Level: beginner
1374: .seealso: `VecSum()`, `VecNorm()`
1375: @*/
1376: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1377: {
1378: PetscInt n;
1382: VecGetSize(v, &n);
1383: VecSum(v, mean);
1384: *mean /= n;
1385: return 0;
1386: }
1388: /*@
1389: VecImaginaryPart - Replaces a complex vector with its imginary part
1391: Collective
1393: Input Parameter:
1394: . v - the vector
1396: Level: beginner
1398: .seealso: `VecNorm()`, `VecRealPart()`
1399: @*/
1400: PetscErrorCode VecImaginaryPart(Vec v)
1401: {
1402: PetscInt i, n;
1403: PetscScalar *x;
1406: VecGetLocalSize(v, &n);
1407: VecGetArray(v, &x);
1408: for (i = 0; i < n; i++) x[i] = PetscImaginaryPart(x[i]);
1409: VecRestoreArray(v, &x);
1410: return 0;
1411: }
1413: /*@
1414: VecRealPart - Replaces a complex vector with its real part
1416: Collective
1418: Input Parameter:
1419: . v - the vector
1421: Level: beginner
1423: .seealso: `VecNorm()`, `VecImaginaryPart()`
1424: @*/
1425: PetscErrorCode VecRealPart(Vec v)
1426: {
1427: PetscInt i, n;
1428: PetscScalar *x;
1431: VecGetLocalSize(v, &n);
1432: VecGetArray(v, &x);
1433: for (i = 0; i < n; i++) x[i] = PetscRealPart(x[i]);
1434: VecRestoreArray(v, &x);
1435: return 0;
1436: }
1438: /*@
1439: VecShift - Shifts all of the components of a vector by computing
1440: `x[i] = x[i] + shift`.
1442: Logically Collective
1444: Input Parameters:
1445: + v - the vector
1446: - shift - the shift
1448: Level: intermediate
1450: @*/
1451: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1452: {
1453: PetscInt i, n;
1454: PetscScalar *x;
1458: VecSetErrorIfLocked(v, 1);
1459: if (shift == 0.0) return 0;
1461: if (v->ops->shift) PetscUseTypeMethod(v, shift, shift);
1462: else {
1463: VecGetLocalSize(v, &n);
1464: VecGetArray(v, &x);
1465: for (i = 0; i < n; i++) x[i] += shift;
1466: VecRestoreArray(v, &x);
1467: }
1468: return 0;
1469: }
1471: /*@
1472: VecAbs - Replaces every element in a vector with its absolute value.
1474: Logically Collective
1476: Input Parameters:
1477: . v - the vector
1479: Level: intermediate
1481: @*/
1482: PetscErrorCode VecAbs(Vec v)
1483: {
1484: PetscInt i, n;
1485: PetscScalar *x;
1488: VecSetErrorIfLocked(v, 1);
1490: if (v->ops->abs) {
1491: PetscUseTypeMethod(v, abs);
1492: } else {
1493: VecGetLocalSize(v, &n);
1494: VecGetArray(v, &x);
1495: for (i = 0; i < n; i++) x[i] = PetscAbsScalar(x[i]);
1496: VecRestoreArray(v, &x);
1497: }
1498: return 0;
1499: }
1501: /*@
1502: VecPermute - Permutes a vector in place using the given ordering.
1504: Input Parameters:
1505: + vec - The vector
1506: . order - The ordering
1507: - inv - The flag for inverting the permutation
1509: Level: beginner
1511: Note: This function does not yet support parallel Index Sets with non-local permutations
1513: .seealso: `MatPermute()`
1514: @*/
1515: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1516: {
1517: const PetscScalar *array;
1518: PetscScalar *newArray;
1519: const PetscInt *idx;
1520: PetscInt i, rstart, rend;
1524: VecSetErrorIfLocked(x, 1);
1525: VecGetOwnershipRange(x, &rstart, &rend);
1526: ISGetIndices(row, &idx);
1527: VecGetArrayRead(x, &array);
1528: PetscMalloc1(x->map->n, &newArray);
1529: if (PetscDefined(USE_DEBUG)) {
1531: }
1532: if (!inv) {
1533: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i] - rstart];
1534: } else {
1535: for (i = 0; i < x->map->n; i++) newArray[idx[i] - rstart] = array[i];
1536: }
1537: VecRestoreArrayRead(x, &array);
1538: ISRestoreIndices(row, &idx);
1539: VecReplaceArray(x, newArray);
1540: return 0;
1541: }
1543: /*@
1544: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1545: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1546: Does NOT take round-off errors into account.
1548: Collective
1550: Input Parameters:
1551: + vec1 - the first vector
1552: - vec2 - the second vector
1554: Output Parameter:
1555: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1557: Level: intermediate
1558: @*/
1559: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1560: {
1561: const PetscScalar *v1, *v2;
1562: PetscInt n1, n2, N1, N2;
1563: PetscBool flg1;
1568: if (vec1 == vec2) *flg = PETSC_TRUE;
1569: else {
1570: VecGetSize(vec1, &N1);
1571: VecGetSize(vec2, &N2);
1572: if (N1 != N2) flg1 = PETSC_FALSE;
1573: else {
1574: VecGetLocalSize(vec1, &n1);
1575: VecGetLocalSize(vec2, &n2);
1576: if (n1 != n2) flg1 = PETSC_FALSE;
1577: else {
1578: VecGetArrayRead(vec1, &v1);
1579: VecGetArrayRead(vec2, &v2);
1580: PetscArraycmp(v1, v2, n1, &flg1);
1581: VecRestoreArrayRead(vec1, &v1);
1582: VecRestoreArrayRead(vec2, &v2);
1583: }
1584: }
1585: /* combine results from all processors */
1586: MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1));
1587: }
1588: return 0;
1589: }
1591: /*@
1592: VecUniqueEntries - Compute the number of unique entries, and those entries
1594: Collective
1596: Input Parameter:
1597: . vec - the vector
1599: Output Parameters:
1600: + n - The number of unique entries
1601: - e - The entries
1603: Level: intermediate
1605: @*/
1606: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1607: {
1608: const PetscScalar *v;
1609: PetscScalar *tmp, *vals;
1610: PetscMPIInt *N, *displs, l;
1611: PetscInt ng, m, i, j, p;
1612: PetscMPIInt size;
1616: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1617: VecGetLocalSize(vec, &m);
1618: VecGetArrayRead(vec, &v);
1619: PetscMalloc2(m, &tmp, size, &N);
1620: for (i = 0, j = 0, l = 0; i < m; ++i) {
1621: /* Can speed this up with sorting */
1622: for (j = 0; j < l; ++j) {
1623: if (v[i] == tmp[j]) break;
1624: }
1625: if (j == l) {
1626: tmp[j] = v[i];
1627: ++l;
1628: }
1629: }
1630: VecRestoreArrayRead(vec, &v);
1631: /* Gather serial results */
1632: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec));
1633: for (p = 0, ng = 0; p < size; ++p) ng += N[p];
1634: PetscMalloc2(ng, &vals, size + 1, &displs);
1635: for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
1636: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec));
1637: /* Find unique entries */
1638: #ifdef PETSC_USE_COMPLEX
1639: SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1640: #else
1641: *n = displs[size];
1642: PetscSortRemoveDupsReal(n, (PetscReal *)vals);
1643: if (e) {
1645: PetscMalloc1(*n, e);
1646: for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
1647: }
1648: PetscFree2(vals, displs);
1649: PetscFree2(tmp, N);
1650: return 0;
1651: #endif
1652: }