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