Actual source code: pvec2.c


  2: /*
  3:      Code for some of the parallel vector primitives.
  4: */
  5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
  6: #include <petscblaslapack.h>

  8: PetscErrorCode VecMDot_MPI(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z)
  9: {
 10:   PetscScalar    awork[128],*work = awork;

 12:   if (nv > 128) {
 13:     PetscMalloc1(nv,&work);
 14:   }
 15:   VecMDot_Seq(xin,nv,y,work);
 16:   MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 17:   if (nv > 128) {
 18:     PetscFree(work);
 19:   }
 20:   return 0;
 21: }

 23: PetscErrorCode VecMTDot_MPI(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z)
 24: {
 25:   PetscScalar    awork[128],*work = awork;

 27:   if (nv > 128) {
 28:     PetscMalloc1(nv,&work);
 29:   }
 30:   VecMTDot_Seq(xin,nv,y,work);
 31:   MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 32:   if (nv > 128) {
 33:     PetscFree(work);
 34:   }
 35:   return 0;
 36: }

 38: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
 39: PetscErrorCode VecNorm_MPI(Vec xin,NormType type,PetscReal *z)
 40: {
 41:   PetscReal         sum,work = 0.0;
 42:   const PetscScalar *xx;
 43:   PetscInt          n   = xin->map->n;
 44:   PetscBLASInt      one = 1,bn = 0;

 46:   PetscBLASIntCast(n,&bn);
 47:   if (type == NORM_2 || type == NORM_FROBENIUS) {
 48:     VecGetArrayRead(xin,&xx);
 49:     work = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
 50:     VecRestoreArrayRead(xin,&xx);
 51:     MPIU_Allreduce(&work,&sum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 52:     *z   = PetscSqrtReal(sum);
 53:     PetscLogFlops(2.0*xin->map->n);
 54:   } else if (type == NORM_1) {
 55:     /* Find the local part */
 56:     VecNorm_Seq(xin,NORM_1,&work);
 57:     /* Find the global max */
 58:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 59:   } else if (type == NORM_INFINITY) {
 60:     /* Find the local max */
 61:     VecNorm_Seq(xin,NORM_INFINITY,&work);
 62:     /* Find the global max */
 63:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
 64:   } else if (type == NORM_1_AND_2) {
 65:     PetscReal temp[2];
 66:     VecNorm_Seq(xin,NORM_1,temp);
 67:     VecNorm_Seq(xin,NORM_2,temp+1);
 68:     temp[1] = temp[1]*temp[1];
 69:     MPIU_Allreduce(temp,z,2,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 70:     z[1] = PetscSqrtReal(z[1]);
 71:   }
 72:   return 0;
 73: }

 75: PetscErrorCode VecMax_MPI(Vec xin,PetscInt *idx,PetscReal *z)
 76: {
 77:   PetscReal      work;

 79:   /* Find the local max */
 80:   VecMax_Seq(xin,idx,&work);
 81: #if defined(PETSC_HAVE_MPIUNI)
 82:   *z = work;
 83: #else
 84:   /* Find the global max */
 85:   if (!idx) {
 86:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
 87:   } else {
 88:     struct { PetscReal v; PetscInt i; } in,out;

 90:     in.v  = work;
 91:     in.i  = *idx + xin->map->rstart;
 92:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)xin));
 93:     *z    = out.v;
 94:     *idx  = out.i;
 95:   }
 96: #endif
 97:   return 0;
 98: }

100: PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z)
101: {
102:   PetscReal      work;

104:   /* Find the local Min */
105:   VecMin_Seq(xin,idx,&work);
106: #if defined(PETSC_HAVE_MPIUNI)
107:   *z = work;
108: #else
109:   /* Find the global Min */
110:   if (!idx) {
111:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));
112:   } else {
113:     struct { PetscReal v; PetscInt i; } in,out;

115:     in.v  = work;
116:     in.i  = *idx + xin->map->rstart;
117:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)xin));
118:     *z    = out.v;
119:     *idx  = out.i;
120:   }
121: #endif
122:   return 0;
123: }