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);
 89: #if defined(PETSC_HAVE_MPIUNI)
 90:   *z = work;
 91: #else
 92:   /* Find the global max */
 93:   if (!idx) {
 94:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
 95:   } else {
 96:     struct { PetscReal v; PetscInt i; } in,out;

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

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

114:   /* Find the local Min */
115:   VecMin_Seq(xin,idx,&work);
116: #if defined(PETSC_HAVE_MPIUNI)
117:   *z = work;
118: #else
119:   /* Find the global Min */
120:   if (!idx) {
121:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));
122:   } else {
123:     struct { PetscReal v; PetscInt i; } in,out;

125:     in.v  = work;
126:     in.i  = *idx + xin->map->rstart;
127:     MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)xin));
128:     *z    = out.v;
129:     *idx  = out.i;
130:   }
131: #endif
132:   return(0);
133: }