1 2 /* 3 Code for some of the parallel vector primitives. 4 */ 5 #include <../src/vec/vec/impls/mpi/pvecimpl.h> 6 #include <petscblaslapack.h> 7 8 PetscErrorCode VecMDot_MPI(Vec xin, PetscInt nv, const Vec y[], PetscScalar *z) { 9 PetscScalar awork[128], *work = awork; 10 11 PetscFunctionBegin; 12 if (nv > 128) PetscCall(PetscMalloc1(nv, &work)); 13 PetscCall(VecMDot_Seq(xin, nv, y, work)); 14 PetscCall(MPIU_Allreduce(work, z, nv, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin))); 15 if (nv > 128) PetscCall(PetscFree(work)); 16 PetscFunctionReturn(0); 17 } 18 19 PetscErrorCode VecMTDot_MPI(Vec xin, PetscInt nv, const Vec y[], PetscScalar *z) { 20 PetscScalar awork[128], *work = awork; 21 22 PetscFunctionBegin; 23 if (nv > 128) PetscCall(PetscMalloc1(nv, &work)); 24 PetscCall(VecMTDot_Seq(xin, nv, y, work)); 25 PetscCall(MPIU_Allreduce(work, z, nv, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin))); 26 if (nv > 128) PetscCall(PetscFree(work)); 27 PetscFunctionReturn(0); 28 } 29 30 #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h> 31 PetscErrorCode VecNorm_MPI(Vec xin, NormType type, PetscReal *z) { 32 PetscReal sum, work = 0.0; 33 const PetscScalar *xx; 34 PetscInt n = xin->map->n; 35 PetscBLASInt one = 1, bn = 0; 36 37 PetscFunctionBegin; 38 PetscCall(PetscBLASIntCast(n, &bn)); 39 if (type == NORM_2 || type == NORM_FROBENIUS) { 40 PetscCall(VecGetArrayRead(xin, &xx)); 41 work = PetscRealPart(BLASdot_(&bn, xx, &one, xx, &one)); 42 PetscCall(VecRestoreArrayRead(xin, &xx)); 43 PetscCall(MPIU_Allreduce(&work, &sum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin))); 44 *z = PetscSqrtReal(sum); 45 PetscCall(PetscLogFlops(2.0 * xin->map->n)); 46 } else if (type == NORM_1) { 47 /* Find the local part */ 48 PetscCall(VecNorm_Seq(xin, NORM_1, &work)); 49 /* Find the global max */ 50 PetscCall(MPIU_Allreduce(&work, z, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin))); 51 } else if (type == NORM_INFINITY) { 52 /* Find the local max */ 53 PetscCall(VecNorm_Seq(xin, NORM_INFINITY, &work)); 54 /* Find the global max */ 55 PetscCall(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 PetscCall(VecNorm_Seq(xin, NORM_1, temp)); 59 PetscCall(VecNorm_Seq(xin, NORM_2, temp + 1)); 60 temp[1] = temp[1] * temp[1]; 61 PetscCall(MPIU_Allreduce(temp, z, 2, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin))); 62 z[1] = PetscSqrtReal(z[1]); 63 } 64 PetscFunctionReturn(0); 65 } 66 67 PetscErrorCode VecMax_MPI(Vec xin, PetscInt *idx, PetscReal *z) { 68 PetscReal work; 69 70 PetscFunctionBegin; 71 /* Find the local max */ 72 PetscCall(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 PetscCall(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; 84 85 in.v = work; 86 in.i = *idx + xin->map->rstart; 87 PetscCall(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 PetscFunctionReturn(0); 93 } 94 95 PetscErrorCode VecMin_MPI(Vec xin, PetscInt *idx, PetscReal *z) { 96 PetscReal work; 97 98 PetscFunctionBegin; 99 /* Find the local Min */ 100 PetscCall(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 PetscCall(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; 112 113 in.v = work; 114 in.i = *idx + xin->map->rstart; 115 PetscCall(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 PetscFunctionReturn(0); 121 } 122