Actual source code: pvec2.c
2: /*
3: Code for some of the parallel vector primitives.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
6: #include <petscblaslapack.h>
8: PetscErrorCode VecMDot_MPI(Vec xin, PetscInt nv, const Vec y[], PetscScalar *z)
9: {
10: PetscScalar awork[128], *work = awork;
12: if (nv > 128) PetscMalloc1(nv, &work);
13: VecMDot_Seq(xin, nv, y, work);
14: MPIU_Allreduce(work, z, nv, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin));
15: if (nv > 128) PetscFree(work);
16: return 0;
17: }
19: PetscErrorCode VecMTDot_MPI(Vec xin, PetscInt nv, const Vec y[], PetscScalar *z)
20: {
21: PetscScalar awork[128], *work = awork;
23: if (nv > 128) PetscMalloc1(nv, &work);
24: VecMTDot_Seq(xin, nv, y, work);
25: MPIU_Allreduce(work, z, nv, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin));
26: if (nv > 128) PetscFree(work);
27: return 0;
28: }
30: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
31: PetscErrorCode VecNorm_MPI(Vec xin, NormType type, PetscReal *z)
32: {
33: PetscReal sum, work = 0.0;
34: const PetscScalar *xx;
35: PetscInt n = xin->map->n;
36: PetscBLASInt one = 1, bn = 0;
38: PetscBLASIntCast(n, &bn);
39: if (type == NORM_2 || type == NORM_FROBENIUS) {
40: VecGetArrayRead(xin, &xx);
41: work = PetscRealPart(BLASdot_(&bn, xx, &one, xx, &one));
42: VecRestoreArrayRead(xin, &xx);
43: MPIU_Allreduce(&work, &sum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin));
44: *z = PetscSqrtReal(sum);
45: PetscLogFlops(2.0 * xin->map->n);
46: } else if (type == NORM_1) {
47: /* Find the local part */
48: VecNorm_Seq(xin, NORM_1, &work);
49: /* Find the global max */
50: MPIU_Allreduce(&work, z, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin));
51: } else if (type == NORM_INFINITY) {
52: /* Find the local max */
53: VecNorm_Seq(xin, NORM_INFINITY, &work);
54: /* Find the global max */
55: MPIU_Allreduce(&work, z, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)xin));
56: } else if (type == NORM_1_AND_2) {
57: PetscReal temp[2];
58: VecNorm_Seq(xin, NORM_1, temp);
59: VecNorm_Seq(xin, NORM_2, temp + 1);
60: temp[1] = temp[1] * temp[1];
61: MPIU_Allreduce(temp, z, 2, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin));
62: z[1] = PetscSqrtReal(z[1]);
63: }
64: return 0;
65: }
67: PetscErrorCode VecMax_MPI(Vec xin, PetscInt *idx, PetscReal *z)
68: {
69: PetscReal work;
71: /* Find the local max */
72: VecMax_Seq(xin, idx, &work);
73: #if defined(PETSC_HAVE_MPIUNI)
74: *z = work;
75: #else
76: /* Find the global max */
77: if (!idx) {
78: MPIU_Allreduce(&work, z, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)xin));
79: } else {
80: struct {
81: PetscReal v;
82: PetscInt i;
83: } in, out;
85: in.v = work;
86: in.i = *idx + xin->map->rstart;
87: MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)xin));
88: *z = out.v;
89: *idx = out.i;
90: }
91: #endif
92: return 0;
93: }
95: PetscErrorCode VecMin_MPI(Vec xin, PetscInt *idx, PetscReal *z)
96: {
97: PetscReal work;
99: /* Find the local Min */
100: VecMin_Seq(xin, idx, &work);
101: #if defined(PETSC_HAVE_MPIUNI)
102: *z = work;
103: #else
104: /* Find the global Min */
105: if (!idx) {
106: MPIU_Allreduce(&work, z, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)xin));
107: } else {
108: struct {
109: PetscReal v;
110: PetscInt i;
111: } in, out;
113: in.v = work;
114: in.i = *idx + xin->map->rstart;
115: MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)xin));
116: *z = out.v;
117: *idx = out.i;
118: }
119: #endif
120: return 0;
121: }