1 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 = PetscMalloc1(nv,&work);CHKERRQ(ierr); 18 } 19 ierr = VecMDot_Seq(xin,nv,y,work);CHKERRQ(ierr); 20 ierr = MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));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 = PetscMalloc1(nv,&work);CHKERRQ(ierr); 37 } 38 ierr = VecMTDot_Seq(xin,nv,y,work);CHKERRQ(ierr); 39 ierr = MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));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 const PetscScalar *xx; 53 PetscErrorCode ierr; 54 PetscInt n = xin->map->n; 55 PetscBLASInt one = 1,bn; 56 57 PetscFunctionBegin; 58 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 59 if (type == NORM_2 || type == NORM_FROBENIUS) { 60 ierr = VecGetArrayRead(xin,&xx);CHKERRQ(ierr); 61 work = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one)); 62 ierr = VecRestoreArrayRead(xin,&xx);CHKERRQ(ierr); 63 ierr = MPIU_Allreduce(&work,&sum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 64 *z = PetscSqrtReal(sum); 65 ierr = PetscLogFlops(2.0*xin->map->n);CHKERRQ(ierr); 66 } else if (type == NORM_1) { 67 /* Find the local part */ 68 ierr = VecNorm_Seq(xin,NORM_1,&work);CHKERRQ(ierr); 69 /* Find the global max */ 70 ierr = MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 71 } else if (type == NORM_INFINITY) { 72 /* Find the local max */ 73 ierr = VecNorm_Seq(xin,NORM_INFINITY,&work);CHKERRQ(ierr); 74 /* Find the global max */ 75 ierr = MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 76 } else if (type == NORM_1_AND_2) { 77 PetscReal temp[2]; 78 ierr = VecNorm_Seq(xin,NORM_1,temp);CHKERRQ(ierr); 79 ierr = VecNorm_Seq(xin,NORM_2,temp+1);CHKERRQ(ierr); 80 temp[1] = temp[1]*temp[1]; 81 ierr = MPIU_Allreduce(temp,z,2,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 82 z[1] = PetscSqrtReal(z[1]); 83 } 84 PetscFunctionReturn(0); 85 } 86 87 extern MPI_Op MPIU_MAXINDEX_OP, MPIU_MININDEX_OP; 88 89 #undef __FUNCT__ 90 #define __FUNCT__ "VecMax_MPI" 91 PetscErrorCode VecMax_MPI(Vec xin,PetscInt *idx,PetscReal *z) 92 { 93 PetscErrorCode ierr; 94 PetscReal work; 95 96 PetscFunctionBegin; 97 /* Find the local max */ 98 ierr = VecMax_Seq(xin,idx,&work);CHKERRQ(ierr); 99 100 /* Find the global max */ 101 if (!idx) { 102 ierr = MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 103 } else { 104 PetscReal work2[2],z2[2]; 105 PetscInt rstart; 106 rstart = xin->map->rstart; 107 work2[0] = work; 108 work2[1] = *idx + rstart; 109 ierr = MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 110 *z = z2[0]; 111 *idx = (PetscInt)z2[1]; 112 } 113 PetscFunctionReturn(0); 114 } 115 116 #undef __FUNCT__ 117 #define __FUNCT__ "VecMin_MPI" 118 PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z) 119 { 120 PetscErrorCode ierr; 121 PetscReal work; 122 123 PetscFunctionBegin; 124 /* Find the local Min */ 125 ierr = VecMin_Seq(xin,idx,&work);CHKERRQ(ierr); 126 127 /* Find the global Min */ 128 if (!idx) { 129 ierr = MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 130 } else { 131 PetscReal work2[2],z2[2]; 132 PetscInt rstart; 133 134 ierr = VecGetOwnershipRange(xin,&rstart,NULL);CHKERRQ(ierr); 135 work2[0] = work; 136 work2[1] = *idx + rstart; 137 ierr = MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr); 138 *z = z2[0]; 139 *idx = (PetscInt)z2[1]; 140 } 141 PetscFunctionReturn(0); 142 } 143 144 145 146 147 148 149 150 151