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) PetscMalloc1(nv, &work);
 13:   VecMDot_Seq(xin, nv, y, work);
 14:   MPIU_Allreduce(work, z, nv, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin));
 15:   if (nv > 128) PetscFree(work);
 16:   return 0;
 17: }

 19: PetscErrorCode VecMTDot_MPI(Vec xin, PetscInt nv, const Vec y[], PetscScalar *z)
 20: {
 21:   PetscScalar awork[128], *work = awork;

 23:   if (nv > 128) PetscMalloc1(nv, &work);
 24:   VecMTDot_Seq(xin, nv, y, work);
 25:   MPIU_Allreduce(work, z, nv, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin));
 26:   if (nv > 128) PetscFree(work);
 27:   return 0;
 28: }

 30: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
 31: PetscErrorCode VecNorm_MPI(Vec xin, NormType type, PetscReal *z)
 32: {
 33:   PetscReal          sum, work = 0.0;
 34:   const PetscScalar *xx;
 35:   PetscInt           n   = xin->map->n;
 36:   PetscBLASInt       one = 1, bn = 0;

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

 67: PetscErrorCode VecMax_MPI(Vec xin, PetscInt *idx, PetscReal *z)
 68: {
 69:   PetscReal work;

 71:   /* Find the local max */
 72:   VecMax_Seq(xin, idx, &work);
 73: #if defined(PETSC_HAVE_MPIUNI)
 74:   *z = work;
 75: #else
 76:   /* Find the global max */
 77:   if (!idx) {
 78:     MPIU_Allreduce(&work, z, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)xin));
 79:   } else {
 80:     struct {
 81:       PetscReal v;
 82:       PetscInt  i;
 83:     } in, out;

 85:     in.v = work;
 86:     in.i = *idx + xin->map->rstart;
 87:     MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)xin));
 88:     *z   = out.v;
 89:     *idx = out.i;
 90:   }
 91: #endif
 92:   return 0;
 93: }

 95: PetscErrorCode VecMin_MPI(Vec xin, PetscInt *idx, PetscReal *z)
 96: {
 97:   PetscReal work;

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

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