1 #define PETSCVEC_DLL 2 /* 3 Code for some of the parallel vector primatives. 4 */ 5 #include "../src/vec/vec/impls/mpi/pvecimpl.h" 6 #include "petscblaslapack.h" 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "VecMDot_MPI" 10 PetscErrorCode VecMDot_MPI(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z) 11 { 12 PetscScalar awork[128],*work = awork; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 if (nv > 128) { 17 ierr = PetscMalloc(nv*sizeof(PetscScalar),&work);CHKERRQ(ierr); 18 } 19 ierr = VecMDot_Seq(xin,nv,y,work);CHKERRQ(ierr); 20 ierr = MPI_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 21 if (nv > 128) { 22 ierr = PetscFree(work);CHKERRQ(ierr); 23 } 24 PetscFunctionReturn(0); 25 } 26 27 #undef __FUNCT__ 28 #define __FUNCT__ "VecMTDot_MPI" 29 PetscErrorCode VecMTDot_MPI(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z) 30 { 31 PetscScalar awork[128],*work = awork; 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 if (nv > 128) { 36 ierr = PetscMalloc(nv*sizeof(PetscScalar),&work);CHKERRQ(ierr); 37 } 38 ierr = VecMTDot_Seq(xin,nv,y,work);CHKERRQ(ierr); 39 ierr = MPI_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 40 if (nv > 128) { 41 ierr = PetscFree(work);CHKERRQ(ierr); 42 } 43 PetscFunctionReturn(0); 44 } 45 46 #include "../src/vec/vec/impls/seq/ftn-kernels/fnorm.h" 47 #undef __FUNCT__ 48 #define __FUNCT__ "VecNorm_MPI" 49 PetscErrorCode VecNorm_MPI(Vec xin,NormType type,PetscReal *z) 50 { 51 PetscReal sum,work = 0.0; 52 PetscScalar *xx = *(PetscScalar**)xin->data; 53 PetscErrorCode ierr; 54 PetscInt n = xin->map->n; 55 PetscBLASInt one = 1,bn = PetscBLASIntCast(n); 56 57 PetscFunctionBegin; 58 if (type == NORM_2 || type == NORM_FROBENIUS) { 59 work = BLASnrm2_(&bn,xx,&one); 60 work *= work; 61 ierr = MPI_Allreduce(&work,&sum,1,MPIU_REAL,MPI_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 62 *z = sqrt(sum); 63 ierr = PetscLogFlops(2.0*xin->map->n);CHKERRQ(ierr); 64 } else if (type == NORM_1) { 65 /* Find the local part */ 66 ierr = VecNorm_Seq(xin,NORM_1,&work);CHKERRQ(ierr); 67 /* Find the global max */ 68 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 69 } else if (type == NORM_INFINITY) { 70 /* Find the local max */ 71 ierr = VecNorm_Seq(xin,NORM_INFINITY,&work);CHKERRQ(ierr); 72 /* Find the global max */ 73 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_MAX,((PetscObject)xin)->comm);CHKERRQ(ierr); 74 } else if (type == NORM_1_AND_2) { 75 PetscReal temp[2]; 76 ierr = VecNorm_Seq(xin,NORM_1,temp);CHKERRQ(ierr); 77 ierr = VecNorm_Seq(xin,NORM_2,temp+1);CHKERRQ(ierr); 78 temp[1] = temp[1]*temp[1]; 79 ierr = MPI_Allreduce(temp,z,2,MPIU_REAL,MPI_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 80 z[1] = sqrt(z[1]); 81 } 82 PetscFunctionReturn(0); 83 } 84 85 /* 86 These two functions are the MPI reduction operation used for max and min with index 87 The call below to MPI_Op_create() converts the function Vec[Max,Min]_Local() to the 88 MPI operator Vec[Max,Min]_Local_Op. 89 */ 90 MPI_Op VecMax_Local_Op = 0; 91 MPI_Op VecMin_Local_Op = 0; 92 93 EXTERN_C_BEGIN 94 #undef __FUNCT__ 95 #define __FUNCT__ "VecMax_Local" 96 void PETSCVEC_DLLEXPORT MPIAPI VecMax_Local(void *in,void *out,PetscMPIInt *cnt,MPI_Datatype *datatype) 97 { 98 PetscReal *xin = (PetscReal *)in,*xout = (PetscReal*)out; 99 100 PetscFunctionBegin; 101 if (*datatype != MPIU_REAL) { 102 (*PetscErrorPrintf)("Can only handle MPIU_REAL data types"); 103 MPI_Abort(MPI_COMM_WORLD,1); 104 } 105 if (xin[0] > xout[0]) { 106 xout[0] = xin[0]; 107 xout[1] = xin[1]; 108 } else if (xin[0] == xout[0]) { 109 xout[1] = PetscMin(xin[1],xout[1]); 110 } 111 PetscFunctionReturnVoid(); /* cannot return a value */ 112 } 113 EXTERN_C_END 114 115 EXTERN_C_BEGIN 116 #undef __FUNCT__ 117 #define __FUNCT__ "VecMin_Local" 118 void PETSCVEC_DLLEXPORT MPIAPI VecMin_Local(void *in,void *out,PetscMPIInt *cnt,MPI_Datatype *datatype) 119 { 120 PetscReal *xin = (PetscReal *)in,*xout = (PetscReal*)out; 121 122 PetscFunctionBegin; 123 if (*datatype != MPIU_REAL) { 124 (*PetscErrorPrintf)("Can only handle MPIU_REAL data types"); 125 MPI_Abort(MPI_COMM_WORLD,1); 126 } 127 if (xin[0] < xout[0]) { 128 xout[0] = xin[0]; 129 xout[1] = xin[1]; 130 } else if (xin[0] == xout[0]) { 131 xout[1] = PetscMin(xin[1],xout[1]); 132 } 133 PetscFunctionReturnVoid(); 134 } 135 EXTERN_C_END 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "VecMax_MPI" 139 PetscErrorCode VecMax_MPI(Vec xin,PetscInt *idx,PetscReal *z) 140 { 141 PetscErrorCode ierr; 142 PetscReal work; 143 144 PetscFunctionBegin; 145 /* Find the local max */ 146 ierr = VecMax_Seq(xin,idx,&work);CHKERRQ(ierr); 147 148 /* Find the global max */ 149 if (!idx) { 150 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_MAX,((PetscObject)xin)->comm);CHKERRQ(ierr); 151 } else { 152 PetscReal work2[2],z2[2]; 153 PetscInt rstart; 154 rstart = xin->map->rstart; 155 work2[0] = work; 156 work2[1] = *idx + rstart; 157 ierr = MPI_Allreduce(work2,z2,2,MPIU_REAL,VecMax_Local_Op,((PetscObject)xin)->comm);CHKERRQ(ierr); 158 *z = z2[0]; 159 *idx = (PetscInt)z2[1]; 160 } 161 PetscFunctionReturn(0); 162 } 163 164 #undef __FUNCT__ 165 #define __FUNCT__ "VecMin_MPI" 166 PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z) 167 { 168 PetscErrorCode ierr; 169 PetscReal work; 170 171 PetscFunctionBegin; 172 /* Find the local Min */ 173 ierr = VecMin_Seq(xin,idx,&work);CHKERRQ(ierr); 174 175 /* Find the global Min */ 176 if (!idx) { 177 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_MIN,((PetscObject)xin)->comm);CHKERRQ(ierr); 178 } else { 179 PetscReal work2[2],z2[2]; 180 PetscInt rstart; 181 182 ierr = VecGetOwnershipRange(xin,&rstart,PETSC_NULL);CHKERRQ(ierr); 183 work2[0] = work; 184 work2[1] = *idx + rstart; 185 ierr = MPI_Allreduce(work2,z2,2,MPIU_REAL,VecMin_Local_Op,((PetscObject)xin)->comm);CHKERRQ(ierr); 186 *z = z2[0]; 187 *idx = (PetscInt)z2[1]; 188 } 189 PetscFunctionReturn(0); 190 } 191 192 193 194 195 196 197 198 199