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; 53 PetscErrorCode ierr; 54 PetscInt n = xin->map->n; 55 56 PetscFunctionBegin; 57 ierr = VecGetArray(xin,&xx);CHKERRQ(ierr); 58 if (type == NORM_2 || type == NORM_FROBENIUS) { 59 60 #if defined(PETSC_HAVE_SLOW_BLAS_NORM2) 61 #if defined(PETSC_USE_FORTRAN_KERNEL_NORM) 62 fortrannormsqr_(xx,&n,&work); 63 #elif defined(PETSC_USE_UNROLLED_NORM) 64 switch (n & 0x3) { 65 case 3: work += PetscRealPart(xx[0]*PetscConj(xx[0])); xx++; 66 case 2: work += PetscRealPart(xx[0]*PetscConj(xx[0])); xx++; 67 case 1: work += PetscRealPart(xx[0]*PetscConj(xx[0])); xx++; n -= 4; 68 } 69 while (n>0) { 70 work += PetscRealPart(xx[0]*PetscConj(xx[0])+xx[1]*PetscConj(xx[1])+ 71 xx[2]*PetscConj(xx[2])+xx[3]*PetscConj(xx[3])); 72 xx += 4; n -= 4; 73 } 74 #else 75 {PetscInt i; for (i=0; i<n; i++) work += PetscRealPart((xx[i])*(PetscConj(xx[i])));} 76 #endif 77 #else 78 {PetscBLASInt one = 1,bn = PetscBLASIntCast(n); 79 work = BLASnrm2_(&bn,xx,&one); 80 work *= work; 81 } 82 #endif 83 ierr = MPI_Allreduce(&work,&sum,1,MPIU_REAL,MPI_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 84 *z = sqrt(sum); 85 ierr = PetscLogFlops(2.0*xin->map->n);CHKERRQ(ierr); 86 } else if (type == NORM_1) { 87 /* Find the local part */ 88 ierr = VecNorm_Seq(xin,NORM_1,&work);CHKERRQ(ierr); 89 /* Find the global max */ 90 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 91 } else if (type == NORM_INFINITY) { 92 /* Find the local max */ 93 ierr = VecNorm_Seq(xin,NORM_INFINITY,&work);CHKERRQ(ierr); 94 /* Find the global max */ 95 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_MAX,((PetscObject)xin)->comm);CHKERRQ(ierr); 96 } else if (type == NORM_1_AND_2) { 97 PetscReal temp[2]; 98 ierr = VecNorm_Seq(xin,NORM_1,temp);CHKERRQ(ierr); 99 ierr = VecNorm_Seq(xin,NORM_2,temp+1);CHKERRQ(ierr); 100 temp[1] = temp[1]*temp[1]; 101 ierr = MPI_Allreduce(temp,z,2,MPIU_REAL,MPI_SUM,((PetscObject)xin)->comm);CHKERRQ(ierr); 102 z[1] = sqrt(z[1]); 103 } 104 ierr = VecRestoreArray(xin,&xx);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 108 /* 109 These two functions are the MPI reduction operation used for max and min with index 110 The call below to MPI_Op_create() converts the function Vec[Max,Min]_Local() to the 111 MPI operator Vec[Max,Min]_Local_Op. 112 */ 113 MPI_Op VecMax_Local_Op = 0; 114 MPI_Op VecMin_Local_Op = 0; 115 116 EXTERN_C_BEGIN 117 #undef __FUNCT__ 118 #define __FUNCT__ "VecMax_Local" 119 void PETSCVEC_DLLEXPORT MPIAPI VecMax_Local(void *in,void *out,PetscMPIInt *cnt,MPI_Datatype *datatype) 120 { 121 PetscReal *xin = (PetscReal *)in,*xout = (PetscReal*)out; 122 123 PetscFunctionBegin; 124 if (*datatype != MPIU_REAL) { 125 (*PetscErrorPrintf)("Can only handle MPIU_REAL data types"); 126 MPI_Abort(MPI_COMM_WORLD,1); 127 } 128 if (xin[0] > xout[0]) { 129 xout[0] = xin[0]; 130 xout[1] = xin[1]; 131 } else if (xin[0] == xout[0]) { 132 xout[1] = PetscMin(xin[1],xout[1]); 133 } 134 PetscFunctionReturnVoid(); /* cannot return a value */ 135 } 136 EXTERN_C_END 137 138 EXTERN_C_BEGIN 139 #undef __FUNCT__ 140 #define __FUNCT__ "VecMin_Local" 141 void PETSCVEC_DLLEXPORT MPIAPI VecMin_Local(void *in,void *out,PetscMPIInt *cnt,MPI_Datatype *datatype) 142 { 143 PetscReal *xin = (PetscReal *)in,*xout = (PetscReal*)out; 144 145 PetscFunctionBegin; 146 if (*datatype != MPIU_REAL) { 147 (*PetscErrorPrintf)("Can only handle MPIU_REAL data types"); 148 MPI_Abort(MPI_COMM_WORLD,1); 149 } 150 if (xin[0] < xout[0]) { 151 xout[0] = xin[0]; 152 xout[1] = xin[1]; 153 } else if (xin[0] == xout[0]) { 154 xout[1] = PetscMin(xin[1],xout[1]); 155 } 156 PetscFunctionReturnVoid(); 157 } 158 EXTERN_C_END 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "VecMax_MPI" 162 PetscErrorCode VecMax_MPI(Vec xin,PetscInt *idx,PetscReal *z) 163 { 164 PetscErrorCode ierr; 165 PetscReal work; 166 167 PetscFunctionBegin; 168 /* Find the local max */ 169 ierr = VecMax_Seq(xin,idx,&work);CHKERRQ(ierr); 170 171 /* Find the global max */ 172 if (!idx) { 173 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_MAX,((PetscObject)xin)->comm);CHKERRQ(ierr); 174 } else { 175 PetscReal work2[2],z2[2]; 176 PetscInt rstart; 177 rstart = xin->map->rstart; 178 work2[0] = work; 179 work2[1] = *idx + rstart; 180 ierr = MPI_Allreduce(work2,z2,2,MPIU_REAL,VecMax_Local_Op,((PetscObject)xin)->comm);CHKERRQ(ierr); 181 *z = z2[0]; 182 *idx = (PetscInt)z2[1]; 183 } 184 PetscFunctionReturn(0); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "VecMin_MPI" 189 PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z) 190 { 191 PetscErrorCode ierr; 192 PetscReal work; 193 194 PetscFunctionBegin; 195 /* Find the local Min */ 196 ierr = VecMin_Seq(xin,idx,&work);CHKERRQ(ierr); 197 198 /* Find the global Min */ 199 if (!idx) { 200 ierr = MPI_Allreduce(&work,z,1,MPIU_REAL,MPI_MIN,((PetscObject)xin)->comm);CHKERRQ(ierr); 201 } else { 202 PetscReal work2[2],z2[2]; 203 PetscInt rstart; 204 205 ierr = VecGetOwnershipRange(xin,&rstart,PETSC_NULL);CHKERRQ(ierr); 206 work2[0] = work; 207 work2[1] = *idx + rstart; 208 ierr = MPI_Allreduce(work2,z2,2,MPIU_REAL,VecMin_Local_Op,((PetscObject)xin)->comm);CHKERRQ(ierr); 209 *z = z2[0]; 210 *idx = (PetscInt)z2[1]; 211 } 212 PetscFunctionReturn(0); 213 } 214 215 216 217 218 219 220 221 222