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;

 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: extern MPI_Op MPIU_MAXINDEX_OP, MPIU_MININDEX_OP;

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

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

 92:   /* Find the global max */
 93:   if (!idx) {
 94:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
 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:     MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)xin));
102:     *z       = z2[0];
103:     *idx     = (PetscInt)z2[1];
104:   }
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);

117:   /* Find the global Min */
118:   if (!idx) {
119:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));
120:   } else {
121:     PetscReal work2[2],z2[2];
122:     PetscInt  rstart;

124:     VecGetOwnershipRange(xin,&rstart,NULL);
125:     work2[0] = work;
126:     work2[1] = *idx + rstart;
127:     MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)xin));
128:     *z   = z2[0];
129:     *idx = (PetscInt)z2[1];
130:   }
131:   return(0);
132: }