vec/utils/vinv.c

/*
     Some useful vector utility functions.
*/
#include <../src/vec/vec/impls/mpi/pvecimpl.h> /*I "petscvec.h" I*/

/*@
   VecStrideSet - Sets a subvector of a vector defined
   by a starting point and a stride with a given value

   Logically Collective

   Input Parameters:
+  v - the vector
.  start - starting point of the subvector (defined by a stride)
-  s - value to set for each entry in that subvector

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   This will only work if the desire subvector is a stride subvector

.seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideScale()`
@*/
PetscErrorCode VecStrideSet(Vec v, PetscInt start, PetscScalar s)
{
  PetscInt     i, n, bs;
  PetscScalar *x;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n  Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
  PetscCall(VecGetArray(v, &x));
  for (i = start; i < n; i += bs) x[i] = s;
  PetscCall(VecRestoreArray(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideScale - Scales a subvector of a vector defined
   by a starting point and a stride.

   Logically Collective

   Input Parameters:
+  v - the vector
.  start - starting point of the subvector (defined by a stride)
-  scale - value to multiply each subvector entry by

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   This will only work if the desire subvector is a stride subvector

.seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideScale()`
@*/
PetscErrorCode VecStrideScale(Vec v, PetscInt start, PetscScalar scale)
{
  PetscInt     i, n, bs;
  PetscScalar *x;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n  Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
  PetscCall(VecGetArray(v, &x));
  for (i = start; i < n; i += bs) x[i] *= scale;
  PetscCall(VecRestoreArray(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideNorm - Computes the norm of subvector of a vector defined
   by a starting point and a stride.

   Collective

   Input Parameters:
+  v - the vector
.  start - starting point of the subvector (defined by a stride)
-  ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`

   Output Parameter:
.  norm - the norm

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If x is the array representing the vector x then this computes the norm
   of the array (x[start],x[start+stride],x[start+2*stride], ....)

   This is useful for computing, say the norm of the pressure variable when
   the pressure is stored (interlaced) with other variables, say density etc.

   This will only work if the desire subvector is a stride subvector

.seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
@*/
PetscErrorCode VecStrideNorm(Vec v, PetscInt start, NormType ntype, PetscReal *nrm)
{
  PetscInt           i, n, bs;
  const PetscScalar *x;
  PetscReal          tnorm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidLogicalCollectiveEnum(v, ntype, 3);
  PetscValidRealPointer(nrm, 4);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
  PetscCall(VecGetArrayRead(v, &x));
  if (ntype == NORM_2) {
    PetscScalar sum = 0.0;
    for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
    tnorm = PetscRealPart(sum);
    PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
    *nrm = PetscSqrtReal(*nrm);
  } else if (ntype == NORM_1) {
    tnorm = 0.0;
    for (i = start; i < n; i += bs) tnorm += PetscAbsScalar(x[i]);
    PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
  } else if (ntype == NORM_INFINITY) {
    tnorm = 0.0;
    for (i = start; i < n; i += bs) {
      if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
    }
    PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
  } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
  PetscCall(VecRestoreArrayRead(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideMax - Computes the maximum of subvector of a vector defined
   by a starting point and a stride and optionally its location.

   Collective

   Input Parameters:
+  v - the vector
-  start - starting point of the subvector (defined by a stride)

   Output Parameters:
+  idex - the location where the maximum occurred  (pass `NULL` if not required)
-  nrm - the maximum value in the subvector

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If xa is the array representing the vector x, then this computes the max
   of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)

   This is useful for computing, say the maximum of the pressure variable when
   the pressure is stored (interlaced) with other variables, e.g., density, etc.
   This will only work if the desire subvector is a stride subvector.

.seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
@*/
PetscErrorCode VecStrideMax(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
{
  PetscInt           i, n, bs, id = -1;
  const PetscScalar *x;
  PetscReal          max = PETSC_MIN_REAL;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidRealPointer(nrm, 4);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
  PetscCall(VecGetArrayRead(v, &x));
  for (i = start; i < n; i += bs) {
    if (PetscRealPart(x[i]) > max) {
      max = PetscRealPart(x[i]);
      id  = i;
    }
  }
  PetscCall(VecRestoreArrayRead(v, &x));
#if defined(PETSC_HAVE_MPIUNI)
  *nrm = max;
  if (idex) *idex = id;
#else
  if (!idex) {
    PetscCall(MPIU_Allreduce(&max, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
  } else {
    struct {
      PetscReal v;
      PetscInt  i;
    } in, out;
    PetscInt rstart;

    PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
    in.v = max;
    in.i = rstart + id;
    PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v)));
    *nrm  = out.v;
    *idex = out.i;
  }
#endif
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideMin - Computes the minimum of subvector of a vector defined
   by a starting point and a stride and optionally its location.

   Collective

   Input Parameters:
+  v - the vector
-  start - starting point of the subvector (defined by a stride)

   Output Parameters:
+  idex - the location where the minimum occurred. (pass `NULL` if not required)
-  nrm - the minimum value in the subvector

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If xa is the array representing the vector x, then this computes the min
   of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)

   This is useful for computing, say the minimum of the pressure variable when
   the pressure is stored (interlaced) with other variables, e.g., density, etc.
   This will only work if the desire subvector is a stride subvector.

.seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
@*/
PetscErrorCode VecStrideMin(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
{
  PetscInt           i, n, bs, id = -1;
  const PetscScalar *x;
  PetscReal          min = PETSC_MAX_REAL;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidRealPointer(nrm, 4);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\nHave you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
  PetscCall(VecGetArrayRead(v, &x));
  for (i = start; i < n; i += bs) {
    if (PetscRealPart(x[i]) < min) {
      min = PetscRealPart(x[i]);
      id  = i;
    }
  }
  PetscCall(VecRestoreArrayRead(v, &x));
#if defined(PETSC_HAVE_MPIUNI)
  *nrm = min;
  if (idex) *idex = id;
#else
  if (!idex) {
    PetscCall(MPIU_Allreduce(&min, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v)));
  } else {
    struct {
      PetscReal v;
      PetscInt  i;
    } in, out;
    PetscInt rstart;

    PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
    in.v = min;
    in.i = rstart + id;
    PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v)));
    *nrm  = out.v;
    *idex = out.i;
  }
#endif
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideSum - Computes the sum of subvector of a vector defined
   by a starting point and a stride.

   Collective

   Input Parameters:
+  v - the vector
-  start - starting point of the subvector (defined by a stride)

   Output Parameter:
.  sum - the sum

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If x is the array representing the vector x then this computes the sum
   of the array (x[start],x[start+stride],x[start+2*stride], ....)

.seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
@*/
PetscErrorCode VecStrideSum(Vec v, PetscInt start, PetscScalar *sum)
{
  PetscInt           i, n, bs;
  const PetscScalar *x;
  PetscScalar        local_sum = 0.0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidScalarPointer(sum, 3);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
  PetscCall(VecGetArrayRead(v, &x));
  for (i = start; i < n; i += bs) local_sum += x[i];
  PetscCall(MPIU_Allreduce(&local_sum, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
  PetscCall(VecRestoreArrayRead(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideScaleAll - Scales the subvectors of a vector defined
   by a starting point and a stride.

   Logically Collective

   Input Parameters:
+  v - the vector
-  scales - values to multiply each subvector entry by

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   The dimension of scales must be the same as the vector block size

.seealso: `Vec`, `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
@*/
PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
{
  PetscInt     i, j, n, bs;
  PetscScalar *x;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidScalarPointer(scales, 2);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCall(VecGetArray(v, &x));
  /* need to provide optimized code for each bs */
  for (i = 0; i < n; i += bs) {
    for (j = 0; j < bs; j++) x[i + j] *= scales[j];
  }
  PetscCall(VecRestoreArray(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideNormAll - Computes the norms of subvectors of a vector defined
   by a starting point and a stride.

   Collective

   Input Parameters:
+  v - the vector
-  ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`

   Output Parameter:
.  nrm - the norms

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If x is the array representing the vector x then this computes the norm
   of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride

   The dimension of nrm must be the same as the vector block size

   This will only work if the desire subvector is a stride subvector

.seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
@*/
PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
{
  PetscInt           i, j, n, bs;
  const PetscScalar *x;
  PetscReal          tnorm[128];
  MPI_Comm           comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidLogicalCollectiveEnum(v, ntype, 2);
  PetscValidRealPointer(nrm, 3);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetArrayRead(v, &x));
  PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

  if (ntype == NORM_2) {
    PetscScalar sum[128];
    for (j = 0; j < bs; j++) sum[j] = 0.0;
    for (i = 0; i < n; i += bs) {
      for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
    }
    for (j = 0; j < bs; j++) tnorm[j] = PetscRealPart(sum[j]);

    PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
    for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
  } else if (ntype == NORM_1) {
    for (j = 0; j < bs; j++) tnorm[j] = 0.0;

    for (i = 0; i < n; i += bs) {
      for (j = 0; j < bs; j++) tnorm[j] += PetscAbsScalar(x[i + j]);
    }

    PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
  } else if (ntype == NORM_INFINITY) {
    PetscReal tmp;
    for (j = 0; j < bs; j++) tnorm[j] = 0.0;

    for (i = 0; i < n; i += bs) {
      for (j = 0; j < bs; j++) {
        if ((tmp = PetscAbsScalar(x[i + j])) > tnorm[j]) tnorm[j] = tmp;
        /* check special case of tmp == NaN */
        if (tmp != tmp) {
          tnorm[j] = tmp;
          break;
        }
      }
    }
    PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_MAX, comm));
  } else SETERRQ(comm, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
  PetscCall(VecRestoreArrayRead(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
   by a starting point and a stride and optionally its location.

   Collective

   Input Parameter:
.  v - the vector

   Output Parameters:
+  index - the location where the maximum occurred (not supported, pass `NULL`,
           if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
-  nrm - the maximum values of each subvector

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   The dimension of nrm must be the same as the vector block size

.seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
@*/
PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
{
  PetscInt           i, j, n, bs;
  const PetscScalar *x;
  PetscReal          max[128], tmp;
  MPI_Comm           comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidRealPointer(nrm, 3);
  PetscCheck(!idex, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetArrayRead(v, &x));
  PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

  if (!n) {
    for (j = 0; j < bs; j++) max[j] = PETSC_MIN_REAL;
  } else {
    for (j = 0; j < bs; j++) max[j] = PetscRealPart(x[j]);

    for (i = bs; i < n; i += bs) {
      for (j = 0; j < bs; j++) {
        if ((tmp = PetscRealPart(x[i + j])) > max[j]) max[j] = tmp;
      }
    }
  }
  PetscCall(MPIU_Allreduce(max, nrm, bs, MPIU_REAL, MPIU_MAX, comm));

  PetscCall(VecRestoreArrayRead(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideMinAll - Computes the minimum of subvector of a vector defined
   by a starting point and a stride and optionally its location.

   Collective

   Input Parameter:
.  v - the vector

   Output Parameters:
+  idex - the location where the minimum occurred (not supported, pass `NULL`,
           if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
-  nrm - the minimums of each subvector

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   The dimension of `nrm` must be the same as the vector block size

.seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
@*/
PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
{
  PetscInt           i, n, bs, j;
  const PetscScalar *x;
  PetscReal          min[128], tmp;
  MPI_Comm           comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidRealPointer(nrm, 3);
  PetscCheck(!idex, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetArrayRead(v, &x));
  PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

  if (!n) {
    for (j = 0; j < bs; j++) min[j] = PETSC_MAX_REAL;
  } else {
    for (j = 0; j < bs; j++) min[j] = PetscRealPart(x[j]);

    for (i = bs; i < n; i += bs) {
      for (j = 0; j < bs; j++) {
        if ((tmp = PetscRealPart(x[i + j])) < min[j]) min[j] = tmp;
      }
    }
  }
  PetscCall(MPIU_Allreduce(min, nrm, bs, MPIU_REAL, MPIU_MIN, comm));

  PetscCall(VecRestoreArrayRead(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideSumAll - Computes the sums of subvectors of a vector defined by a stride.

   Collective

   Input Parameter:
.  v - the vector

   Output Parameter:
.  sums - the sums

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If x is the array representing the vector x then this computes the sum
   of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride

.seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
@*/
PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
{
  PetscInt           i, j, n, bs;
  const PetscScalar *x;
  PetscScalar        local_sums[128];
  MPI_Comm           comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidScalarPointer(sums, 2);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetArrayRead(v, &x));
  PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

  for (j = 0; j < bs; j++) local_sums[j] = 0.0;
  for (i = 0; i < n; i += bs) {
    for (j = 0; j < bs; j++) local_sums[j] += x[i + j];
  }
  PetscCall(MPIU_Allreduce(local_sums, sums, bs, MPIU_SCALAR, MPIU_SUM, comm));

  PetscCall(VecRestoreArrayRead(v, &x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*----------------------------------------------------------------------------------------------*/
/*@
   VecStrideGatherAll - Gathers all the single components from a multi-component vector into
   separate vectors.

   Collective

   Input Parameters:
+  v - the vector
-  addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

   Output Parameter:
.  s - the location where the subvectors are stored

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If x is the array representing the vector x then this gathers
   the arrays (x[start],x[start+stride],x[start+2*stride], ....)
   for start=0,1,2,...bs-1

   The parallel layout of the vector and the subvector must be the same;
   i.e., nlocal of v = stride*(nlocal of s)

   Not optimized; could be easily

.seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
          `VecStrideScatterAll()`
@*/
PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
{
  PetscInt           i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
  PetscScalar      **y;
  const PetscScalar *x;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidPointer(s, 2);
  PetscValidHeaderSpecific(*s, VEC_CLASSID, 2);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetLocalSize(s[0], &n2));
  PetscCall(VecGetArrayRead(v, &x));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");

  PetscCall(PetscMalloc2(bs, &y, bs, &bss));
  nv  = 0;
  nvc = 0;
  for (i = 0; i < bs; i++) {
    PetscCall(VecGetBlockSize(s[i], &bss[i]));
    if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
    PetscCall(VecGetArray(s[i], &y[i]));
    nvc += bss[i];
    nv++;
    PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
    if (nvc == bs) break;
  }

  n = n / bs;

  jj = 0;
  if (addv == INSERT_VALUES) {
    for (j = 0; j < nv; j++) {
      for (k = 0; k < bss[j]; k++) {
        for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
      }
      jj += bss[j];
    }
  } else if (addv == ADD_VALUES) {
    for (j = 0; j < nv; j++) {
      for (k = 0; k < bss[j]; k++) {
        for (i = 0; i < n; i++) y[j][i * bss[j] + k] += x[bs * i + jj + k];
      }
      jj += bss[j];
    }
#if !defined(PETSC_USE_COMPLEX)
  } else if (addv == MAX_VALUES) {
    for (j = 0; j < nv; j++) {
      for (k = 0; k < bss[j]; k++) {
        for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
      }
      jj += bss[j];
    }
#endif
  } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

  PetscCall(VecRestoreArrayRead(v, &x));
  for (i = 0; i < nv; i++) PetscCall(VecRestoreArray(s[i], &y[i]));

  PetscCall(PetscFree2(y, bss));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideScatterAll - Scatters all the single components from separate vectors into
     a multi-component vector.

   Collective

   Input Parameters:
+  s - the location where the subvectors are stored
-  addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

   Output Parameter:
.  v - the multicomponent vector

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   The parallel layout of the vector and the subvector must be the same;
   i.e., nlocal of v = stride*(nlocal of s)

   Not optimized; could be easily

.seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
          `VecStrideScatterAll()`
@*/
PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
{
  PetscInt            i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
  PetscScalar        *x;
  PetscScalar const **y;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 2);
  PetscValidPointer(s, 1);
  PetscValidHeaderSpecific(*s, VEC_CLASSID, 1);
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetLocalSize(s[0], &n2));
  PetscCall(VecGetArray(v, &x));
  PetscCall(VecGetBlockSize(v, &bs));
  PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");

  PetscCall(PetscMalloc2(bs, (PetscScalar ***)&y, bs, &bss));
  nv  = 0;
  nvc = 0;
  for (i = 0; i < bs; i++) {
    PetscCall(VecGetBlockSize(s[i], &bss[i]));
    if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
    PetscCall(VecGetArrayRead(s[i], &y[i]));
    nvc += bss[i];
    nv++;
    PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
    if (nvc == bs) break;
  }

  n = n / bs;

  jj = 0;
  if (addv == INSERT_VALUES) {
    for (j = 0; j < nv; j++) {
      for (k = 0; k < bss[j]; k++) {
        for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
      }
      jj += bss[j];
    }
  } else if (addv == ADD_VALUES) {
    for (j = 0; j < nv; j++) {
      for (k = 0; k < bss[j]; k++) {
        for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
      }
      jj += bss[j];
    }
#if !defined(PETSC_USE_COMPLEX)
  } else if (addv == MAX_VALUES) {
    for (j = 0; j < nv; j++) {
      for (k = 0; k < bss[j]; k++) {
        for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
      }
      jj += bss[j];
    }
#endif
  } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

  PetscCall(VecRestoreArray(v, &x));
  for (i = 0; i < nv; i++) PetscCall(VecRestoreArrayRead(s[i], &y[i]));
  PetscCall(PetscFree2(*(PetscScalar ***)&y, bss));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideGather - Gathers a single component from a multi-component vector into
   another vector.

   Collective

   Input Parameters:
+  v - the vector
.  start - starting point of the subvector (defined by a stride)
-  addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

   Output Parameter:
.  s - the location where the subvector is stored

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   If x is the array representing the vector x then this gathers
   the array (x[start],x[start+stride],x[start+2*stride], ....)

   The parallel layout of the vector and the subvector must be the same;
   i.e., nlocal of v = stride*(nlocal of s)

   Not optimized; could be easily

.seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
          `VecStrideScatterAll()`
@*/
PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidHeaderSpecific(s, VEC_CLASSID, 3);
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
             v->map->bs);
  PetscUseTypeMethod(v, stridegather, start, s, addv);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

   Collective

   Input Parameters:
+  s - the single-component vector
.  start - starting point of the subvector (defined by a stride)
-  addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

   Output Parameter:
.  v - the location where the subvector is scattered (the multi-component vector)

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` on the multi-component vector before this
   routine to set the stride  information, or use a vector created from a multicomponent `DMDA`.

   The parallel layout of the vector and the subvector must be the same;
   i.e., nlocal of v = stride*(nlocal of s)

   Not optimized; could be easily

.seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
          `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
@*/
PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(s, VEC_CLASSID, 1);
  PetscValidHeaderSpecific(v, VEC_CLASSID, 3);
  PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
  PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride\n Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
             v->map->bs);
  PetscCall((*v->ops->stridescatter)(s, start, v, addv));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
   another vector.

   Collective

   Input Parameters:
+  v - the vector
.  nidx - the number of indices
.  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
.  idxs - the indices of the components 0 <= idxs[0] ...idxs[nidx-1] < bs(s), they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
-  addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

   Output Parameter:
.  s - the location where the subvector is stored

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` on both vectors before this routine to set the stride
   information, or use a vector created from a multicomponent `DMDA`.

   The parallel layout of the vector and the subvector must be the same;

   Not optimized; could be easily

.seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
          `VecStrideScatterAll()`
@*/
PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidHeaderSpecific(s, VEC_CLASSID, 5);
  if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
  PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.

   Collective

   Input Parameters:
+  s - the smaller-component vector
.  nidx - the number of indices in idx
.  idxs - the indices of the components in the smaller-component vector, 0 <= idxs[0] ...idxs[nidx-1] < bs(s) they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
.  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
-  addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

   Output Parameter:
.  v - the location where the subvector is into scattered (the multi-component vector)

   Level: advanced

   Notes:
   One must call `VecSetBlockSize()` on the vectors before this
   routine to set the stride  information, or use a vector created from a multicomponent `DMDA`.

   The parallel layout of the vector and the subvector must be the same;

   Not optimized; could be easily

.seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
          `VecStrideScatterAll()`
@*/
PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(s, VEC_CLASSID, 1);
  PetscValidHeaderSpecific(v, VEC_CLASSID, 5);
  if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
  PetscCall((*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
{
  PetscInt           i, n, bs, ns;
  const PetscScalar *x;
  PetscScalar       *y;

  PetscFunctionBegin;
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetLocalSize(s, &ns));
  PetscCall(VecGetArrayRead(v, &x));
  PetscCall(VecGetArray(s, &y));

  bs = v->map->bs;
  PetscCheck(n == ns * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Subvector length * blocksize %" PetscInt_FMT " not correct for gather from original vector %" PetscInt_FMT, ns * bs, n);
  x += start;
  n = n / bs;

  if (addv == INSERT_VALUES) {
    for (i = 0; i < n; i++) y[i] = x[bs * i];
  } else if (addv == ADD_VALUES) {
    for (i = 0; i < n; i++) y[i] += x[bs * i];
#if !defined(PETSC_USE_COMPLEX)
  } else if (addv == MAX_VALUES) {
    for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
#endif
  } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

  PetscCall(VecRestoreArrayRead(v, &x));
  PetscCall(VecRestoreArray(s, &y));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
{
  PetscInt           i, n, bs, ns;
  PetscScalar       *x;
  const PetscScalar *y;

  PetscFunctionBegin;
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetLocalSize(s, &ns));
  PetscCall(VecGetArray(v, &x));
  PetscCall(VecGetArrayRead(s, &y));

  bs = v->map->bs;
  PetscCheck(n == ns * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Subvector length * blocksize %" PetscInt_FMT " not correct for scatter to multicomponent vector %" PetscInt_FMT, ns * bs, n);
  x += start;
  n = n / bs;

  if (addv == INSERT_VALUES) {
    for (i = 0; i < n; i++) x[bs * i] = y[i];
  } else if (addv == ADD_VALUES) {
    for (i = 0; i < n; i++) x[bs * i] += y[i];
#if !defined(PETSC_USE_COMPLEX)
  } else if (addv == MAX_VALUES) {
    for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
#endif
  } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

  PetscCall(VecRestoreArray(v, &x));
  PetscCall(VecRestoreArrayRead(s, &y));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
{
  PetscInt           i, j, n, bs, bss, ns;
  const PetscScalar *x;
  PetscScalar       *y;

  PetscFunctionBegin;
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetLocalSize(s, &ns));
  PetscCall(VecGetArrayRead(v, &x));
  PetscCall(VecGetArray(s, &y));

  bs  = v->map->bs;
  bss = s->map->bs;
  n   = n / bs;

  if (PetscDefined(USE_DEBUG)) {
    PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
    for (j = 0; j < nidx; j++) {
      PetscCheck(idxv[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxv[j]);
      PetscCheck(idxv[j] < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is greater than or equal to vector blocksize %" PetscInt_FMT, j, idxv[j], bs);
    }
    PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not gathering into all locations");
  }

  if (addv == INSERT_VALUES) {
    if (!idxs) {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
      }
    } else {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
      }
    }
  } else if (addv == ADD_VALUES) {
    if (!idxs) {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
      }
    } else {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
      }
    }
#if !defined(PETSC_USE_COMPLEX)
  } else if (addv == MAX_VALUES) {
    if (!idxs) {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
      }
    } else {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
      }
    }
#endif
  } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

  PetscCall(VecRestoreArrayRead(v, &x));
  PetscCall(VecRestoreArray(s, &y));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
{
  PetscInt           j, i, n, bs, ns, bss;
  PetscScalar       *x;
  const PetscScalar *y;

  PetscFunctionBegin;
  PetscCall(VecGetLocalSize(v, &n));
  PetscCall(VecGetLocalSize(s, &ns));
  PetscCall(VecGetArray(v, &x));
  PetscCall(VecGetArrayRead(s, &y));

  bs  = v->map->bs;
  bss = s->map->bs;
  n   = n / bs;

  if (PetscDefined(USE_DEBUG)) {
    PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
    for (j = 0; j < bss; j++) {
      if (idxs) {
        PetscCheck(idxs[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxs[j]);
        PetscCheck(idxs[j] < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is greater than or equal to vector blocksize %" PetscInt_FMT, j, idxs[j], bs);
      }
    }
    PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not scattering from all locations");
  }

  if (addv == INSERT_VALUES) {
    if (!idxs) {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
      }
    } else {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
      }
    }
  } else if (addv == ADD_VALUES) {
    if (!idxs) {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
      }
    } else {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
      }
    }
#if !defined(PETSC_USE_COMPLEX)
  } else if (addv == MAX_VALUES) {
    if (!idxs) {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
      }
    } else {
      for (i = 0; i < n; i++) {
        for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
      }
    }
#endif
  } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown insert type");

  PetscCall(VecRestoreArray(v, &x));
  PetscCall(VecRestoreArrayRead(s, &y));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode VecApplyUnary_Private(Vec v, PetscErrorCode (*unary_op)(Vec), PetscScalar (*UnaryFunc)(PetscScalar))
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecSetErrorIfLocked(v, 1));
  if (unary_op) {
    PetscValidFunction(unary_op, 2);
    PetscCall((*unary_op)(v));
  } else {
    PetscInt     n;
    PetscScalar *x;

    PetscValidFunction(UnaryFunc, 3);
    PetscCall(VecGetLocalSize(v, &n));
    PetscCall(VecGetArray(v, &x));
    for (PetscInt i = 0; i < n; ++i) x[i] = UnaryFunc(x[i]);
    PetscCall(VecRestoreArray(v, &x));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscScalar ScalarReciprocal_Fn(PetscScalar x)
{
  const PetscScalar zero = 0.0;

  return x == zero ? zero : ((PetscScalar)1.0) / x;
}

PetscErrorCode VecReciprocal_Default(Vec v)
{
  PetscFunctionBegin;
  PetscCall(VecApplyUnary_Private(v, NULL, ScalarReciprocal_Fn));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscScalar ScalarExp_Fn(PetscScalar x)
{
  return PetscExpScalar(x);
}

/*@
  VecExp - Replaces each component of a vector by e^x_i

  Not Collective

  Input Parameter:
. v - The vector

  Output Parameter:
. v - The vector of exponents

  Level: beginner

.seealso: `Vec`, `VecLog()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`

@*/
PetscErrorCode VecExp(Vec v)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecApplyUnary_Private(v, v->ops->exp, ScalarExp_Fn));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscScalar ScalarLog_Fn(PetscScalar x)
{
  return PetscLogScalar(x);
}

/*@
  VecLog - Replaces each component of a vector by log(x_i), the natural logarithm

  Not Collective

  Input Parameter:
. v - The vector

  Output Parameter:
. v - The vector of logs

  Level: beginner

.seealso: `Vec`, `VecExp()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`

@*/
PetscErrorCode VecLog(Vec v)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecApplyUnary_Private(v, v->ops->log, ScalarLog_Fn));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscScalar ScalarAbs_Fn(PetscScalar x)
{
  return PetscAbsScalar(x);
}

/*@
   VecAbs - Replaces every element in a vector with its absolute value.

   Logically Collective

   Input Parameters:
.  v - the vector

   Level: intermediate

@*/
PetscErrorCode VecAbs(Vec v)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecApplyUnary_Private(v, v->ops->abs, ScalarAbs_Fn));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscScalar ScalarSqrtAbs_Fn(PetscScalar x)
{
  return PetscSqrtScalar(ScalarAbs_Fn(x));
}

/*@
  VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.

  Not Collective

  Input Parameter:
. v - The vector

  Level: beginner

  Note:
  The actual function is sqrt(|x_i|)

.seealso: `Vec`, `VecLog()`, `VecExp()`, `VecReciprocal()`, `VecAbs()`

@*/
PetscErrorCode VecSqrtAbs(Vec v)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecApplyUnary_Private(v, v->ops->sqrt, ScalarSqrtAbs_Fn));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscScalar ScalarImaginaryPart_Fn(PetscScalar x)
{
  const PetscReal imag = PetscImaginaryPart(x);

#if PetscDefined(USE_COMPLEX)
  return PetscCMPLX(imag, 0.0);
#else
  return imag;
#endif
}

/*@
   VecImaginaryPart - Replaces a complex vector with its imginary part

   Collective

   Input Parameter:
.  v - the vector

   Level: beginner

.seealso: `Vec`, `VecNorm()`, `VecRealPart()`
@*/
PetscErrorCode VecImaginaryPart(Vec v)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecApplyUnary_Private(v, NULL, ScalarImaginaryPart_Fn));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscScalar ScalarRealPart_Fn(PetscScalar x)
{
  const PetscReal real = PetscRealPart(x);

#if PetscDefined(USE_COMPLEX)
  return PetscCMPLX(real, 0.0);
#else
  return real;
#endif
}

/*@
   VecRealPart - Replaces a complex vector with its real part

   Collective

   Input Parameter:
.  v - the vector

   Level: beginner

.seealso: `Vec`, `VecNorm()`, `VecImaginaryPart()`
@*/
PetscErrorCode VecRealPart(Vec v)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscCall(VecApplyUnary_Private(v, NULL, ScalarRealPart_Fn));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector

  Collective

  Input Parameters:
+ s - first vector
- t - second vector

  Output Parameters:
+ dp - s'conj(t)
- nm - t'conj(t)

  Level: advanced

  Note:
    conj(x) is the complex conjugate of x when x is complex

.seealso: `Vec`, `VecDot()`, `VecNorm()`, `VecDotBegin()`, `VecNormBegin()`, `VecDotEnd()`, `VecNormEnd()`

@*/
PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
{
  PetscScalar work[] = {0.0, 0.0};

  PetscFunctionBegin;
  PetscValidHeaderSpecific(s, VEC_CLASSID, 1);
  PetscValidHeaderSpecific(t, VEC_CLASSID, 2);
  PetscValidScalarPointer(dp, 3);
  PetscValidRealPointer(nm, 4);
  PetscValidType(s, 1);
  PetscValidType(t, 2);
  PetscCheckSameTypeAndComm(s, 1, t, 2);
  PetscCheck(s->map->N == t->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector global lengths");
  PetscCheck(s->map->n == t->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector local lengths");

  PetscCall(PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0));
  if (s->ops->dotnorm2) {
    PetscUseTypeMethod(s, dotnorm2, t, work, work + 1);
  } else {
    const PetscScalar *sx, *tx;
    PetscInt           n;

    PetscCall(VecGetLocalSize(s, &n));
    PetscCall(VecGetArrayRead(s, &sx));
    PetscCall(VecGetArrayRead(t, &tx));
    for (PetscInt i = 0; i < n; ++i) {
      const PetscScalar txconj = PetscConj(tx[i]);

      work[0] += sx[i] * txconj;
      work[1] += tx[i] * txconj;
    }
    PetscCall(VecRestoreArrayRead(t, &tx));
    PetscCall(VecRestoreArrayRead(s, &sx));
    PetscCall(MPIU_Allreduce(MPI_IN_PLACE, work, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
    PetscCall(PetscLogFlops(4.0 * n));
  }
  PetscCall(PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0));
  *dp = work[0];
  *nm = PetscRealPart(work[1]);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecSum - Computes the sum of all the components of a vector.

   Collective

   Input Parameter:
.  v - the vector

   Output Parameter:
.  sum - the result

   Level: beginner

.seealso: `Vec`, `VecMean()`, `VecNorm()`
@*/
PetscErrorCode VecSum(Vec v, PetscScalar *sum)
{
  PetscScalar tmp = 0.0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidScalarPointer(sum, 2);
  if (v->ops->sum) {
    PetscUseTypeMethod(v, sum, &tmp);
  } else {
    const PetscScalar *x;
    PetscInt           n;

    PetscCall(VecGetLocalSize(v, &n));
    PetscCall(VecGetArrayRead(v, &x));
    for (PetscInt i = 0; i < n; ++i) tmp += x[i];
    PetscCall(VecRestoreArrayRead(v, &x));
  }
  PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &tmp, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
  *sum = tmp;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecMean - Computes the arithmetic mean of all the components of a vector.

   Collective

   Input Parameter:
.  v - the vector

   Output Parameter:
.  mean - the result

   Level: beginner

.seealso: `Vec`, `VecSum()`, `VecNorm()`
@*/
PetscErrorCode VecMean(Vec v, PetscScalar *mean)
{
  PetscInt n;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidScalarPointer(mean, 2);
  PetscCall(VecGetSize(v, &n));
  PetscCall(VecSum(v, mean));
  *mean /= n;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecShift - Shifts all of the components of a vector by computing
   `x[i] = x[i] + shift`.

   Logically Collective

   Input Parameters:
+  v - the vector
-  shift - the shift

   Level: intermediate

.seealso: `Vec`
@*/
PetscErrorCode VecShift(Vec v, PetscScalar shift)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
  PetscValidLogicalCollectiveScalar(v, shift, 2);
  PetscCall(VecSetErrorIfLocked(v, 1));
  if (shift == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);

  if (v->ops->shift) {
    PetscUseTypeMethod(v, shift, shift);
  } else {
    PetscInt     n;
    PetscScalar *x;

    PetscCall(VecGetLocalSize(v, &n));
    PetscCall(VecGetArray(v, &x));
    for (PetscInt i = 0; i < n; ++i) x[i] += shift;
    PetscCall(VecRestoreArray(v, &x));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  VecPermute - Permutes a vector in place using the given ordering.

  Input Parameters:
+ vec   - The vector
. order - The ordering
- inv   - The flag for inverting the permutation

  Level: beginner

  Note:
  This function does not yet support parallel Index Sets with non-local permutations

.seealso: `Vec`, `MatPermute()`
@*/
PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
{
  const PetscScalar *array;
  PetscScalar       *newArray;
  const PetscInt    *idx;
  PetscInt           i, rstart, rend;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(x, VEC_CLASSID, 1);
  PetscValidHeaderSpecific(row, IS_CLASSID, 2);
  PetscCall(VecSetErrorIfLocked(x, 1));
  PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
  PetscCall(ISGetIndices(row, &idx));
  PetscCall(VecGetArrayRead(x, &array));
  PetscCall(PetscMalloc1(x->map->n, &newArray));
  if (PetscDefined(USE_DEBUG)) {
    for (i = 0; i < x->map->n; i++) PetscCheck(!(idx[i] < rstart) && !(idx[i] >= rend), PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Permutation index %" PetscInt_FMT " is out of bounds: %" PetscInt_FMT, i, idx[i]);
  }
  if (!inv) {
    for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i] - rstart];
  } else {
    for (i = 0; i < x->map->n; i++) newArray[idx[i] - rstart] = array[i];
  }
  PetscCall(VecRestoreArrayRead(x, &array));
  PetscCall(ISRestoreIndices(row, &idx));
  PetscCall(VecReplaceArray(x, newArray));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
   or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
   Does NOT take round-off errors into account.

   Collective

   Input Parameters:
+  vec1 - the first vector
-  vec2 - the second vector

   Output Parameter:
.  flg - `PETSC_TRUE` if the vectors are equal; `PETSC_FALSE` otherwise.

   Level: intermediate

.seealso: `Vec`
@*/
PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
{
  const PetscScalar *v1, *v2;
  PetscInt           n1, n2, N1, N2;
  PetscBool          flg1;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(vec1, VEC_CLASSID, 1);
  PetscValidHeaderSpecific(vec2, VEC_CLASSID, 2);
  PetscValidBoolPointer(flg, 3);
  if (vec1 == vec2) *flg = PETSC_TRUE;
  else {
    PetscCall(VecGetSize(vec1, &N1));
    PetscCall(VecGetSize(vec2, &N2));
    if (N1 != N2) flg1 = PETSC_FALSE;
    else {
      PetscCall(VecGetLocalSize(vec1, &n1));
      PetscCall(VecGetLocalSize(vec2, &n2));
      if (n1 != n2) flg1 = PETSC_FALSE;
      else {
        PetscCall(VecGetArrayRead(vec1, &v1));
        PetscCall(VecGetArrayRead(vec2, &v2));
        PetscCall(PetscArraycmp(v1, v2, n1, &flg1));
        PetscCall(VecRestoreArrayRead(vec1, &v1));
        PetscCall(VecRestoreArrayRead(vec2, &v2));
      }
    }
    /* combine results from all processors */
    PetscCall(MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1)));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
   VecUniqueEntries - Compute the number of unique entries, and those entries

   Collective

   Input Parameter:
.  vec - the vector

   Output Parameters:
+  n - The number of unique entries
-  e - The entries

   Level: intermediate

.seealso: `Vec`
@*/
PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
{
  const PetscScalar *v;
  PetscScalar       *tmp, *vals;
  PetscMPIInt       *N, *displs, l;
  PetscInt           ng, m, i, j, p;
  PetscMPIInt        size;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(vec, VEC_CLASSID, 1);
  PetscValidIntPointer(n, 2);
  PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
  PetscCall(VecGetLocalSize(vec, &m));
  PetscCall(VecGetArrayRead(vec, &v));
  PetscCall(PetscMalloc2(m, &tmp, size, &N));
  for (i = 0, j = 0, l = 0; i < m; ++i) {
    /* Can speed this up with sorting */
    for (j = 0; j < l; ++j) {
      if (v[i] == tmp[j]) break;
    }
    if (j == l) {
      tmp[j] = v[i];
      ++l;
    }
  }
  PetscCall(VecRestoreArrayRead(vec, &v));
  /* Gather serial results */
  PetscCallMPI(MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec)));
  for (p = 0, ng = 0; p < size; ++p) ng += N[p];
  PetscCall(PetscMalloc2(ng, &vals, size + 1, &displs));
  for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
  PetscCallMPI(MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec)));
  /* Find unique entries */
#ifdef PETSC_USE_COMPLEX
  SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
#else
  *n = displs[size];
  PetscCall(PetscSortRemoveDupsReal(n, (PetscReal *)vals));
  if (e) {
    PetscValidPointer(e, 3);
    PetscCall(PetscMalloc1(*n, e));
    for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
  }
  PetscCall(PetscFree2(vals, displs));
  PetscCall(PetscFree2(tmp, N));
  PetscFunctionReturn(PETSC_SUCCESS);
#endif
}