Actual source code: pvec2.c


  2: /*
  3:      Code for some of the parallel vector primatives.
  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;

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

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

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

 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;

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

 81: PetscErrorCode VecMax_MPI(Vec xin,PetscInt *idx,PetscReal *z)
 82: {
 84:   PetscReal      work;

 87:   /* Find the local max */
 88:   VecMax_Seq(xin,idx,&work);

 90:   /* Find the global max */
 91:   if (!idx) {
 92:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
 93:   } else {
 94:     struct { PetscReal v; PetscInt i; } in,out;

 96:     in.v  = work;
 97:     in.i  = *idx + xin->map->rstart;
 98:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)xin));
 99:     *z    = out.v;
100:     *idx  = out.i;
101:   }
102:   return(0);
103: }

105: PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z)
106: {
108:   PetscReal      work;

111:   /* Find the local Min */
112:   VecMin_Seq(xin,idx,&work);

114:   /* Find the global Min */
115:   if (!idx) {
116:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));
117:   } else {
118:     struct { PetscReal v; PetscInt i; } in,out;

120:     in.v  = work;
121:     in.i  = *idx + xin->map->rstart;
122:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)xin));
123:     *z    = out.v;
124:     *idx  = out.i;
125:   }
126:   return(0);
127: }