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:     PetscReal work2[2],z2[2];
     95:     PetscInt  rstart;
     96:     rstart   = xin->map->rstart;
     97:     work2[0] = work;
     98:     work2[1] = *idx + rstart;
     99:     MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)xin));
    100:     *z       = z2[0];
    101:     *idx     = (PetscInt)z2[1];
    102:   }
    103:   return(0);
    104: }

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

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

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

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