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