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