xref: /petsc/src/mat/tests/ex23.c (revision 4f58015e5b963ff70f3b7a63a399e8a46849015e)

static char help[] = "Tests the use of interface functions for MATIS matrices and conversion routines.\n";

#include <petscmat.h>

PetscErrorCode TestMatZeroRows(Mat, Mat, PetscBool, IS, PetscScalar);
PetscErrorCode CheckMat(Mat, Mat, PetscBool, const char *);

int main(int argc, char **args)
{
  Mat                    A, B, A2, B2, T;
  Mat                    Aee, Aeo, Aoe, Aoo;
  Mat                   *mats, *Asub, *Bsub;
  Vec                    x, y;
  MatInfo                info;
  ISLocalToGlobalMapping cmap, rmap;
  IS                     is, is2, reven, rodd, ceven, codd;
  IS                    *rows, *cols;
  IS                     irow[2], icol[2];
  PetscLayout            rlayout, clayout;
  const PetscInt        *rrange, *crange;
  MatType                lmtype;
  PetscScalar            diag = 2.;
  PetscInt               n, m, i, lm, ln;
  PetscInt               rst, ren, cst, cen, nr, nc;
  PetscMPIInt            rank, size, lrank, rrank;
  PetscBool              testT, squaretest, isaij;
  PetscBool              permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE, allow_repeated = PETSC_TRUE;
  PetscBool              diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
  PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
  m = n = 2 * size;
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &symmetric, NULL));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-negmap", &negmap, NULL));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-repmap", &repmap, NULL));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-permmap", &permute, NULL));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-diffmap", &diffmap, NULL));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-allow_repeated", &allow_repeated, NULL));
  PetscCheck(size == 1 || m >= 4, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 4 for parallel runs");
  PetscCheck(size != 1 || m >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 2 for uniprocessor runs");
  PetscCheck(n >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of cols should be larger or equal 2");

  /* create a MATIS matrix */
  PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
  PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
  PetscCall(MatSetType(A, MATIS));
  PetscCall(MatSetFromOptions(A));
  if (!negmap && !repmap) {
    /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines
       Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces
       Equivalent to passing NULL for the mapping */
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &is));
  } else if (negmap && !repmap) { /* non repeated but with negative indices */
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 2, -2, 1, &is));
  } else if (!negmap && repmap) { /* non negative but repeated indices */
    IS isl[2];

    PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &isl[0]));
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, n - 1, -1, &isl[1]));
    PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
    PetscCall(ISDestroy(&isl[0]));
    PetscCall(ISDestroy(&isl[1]));
  } else { /* negative and repeated indices */
    IS isl[2];

    PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, -1, 1, &isl[0]));
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, n - 1, -1, &isl[1]));
    PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
    PetscCall(ISDestroy(&isl[0]));
    PetscCall(ISDestroy(&isl[1]));
  }
  PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmap));
  PetscCall(ISDestroy(&is));

  if (m != n || diffmap) {
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, m, permute ? m - 1 : 0, permute ? -1 : 1, &is));
    PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmap));
    PetscCall(ISDestroy(&is));
  } else {
    PetscCall(PetscObjectReference((PetscObject)cmap));
    rmap = cmap;
  }

  PetscCall(MatISSetAllowRepeated(A, allow_repeated));
  PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
  PetscCall(MatISStoreL2L(A, PETSC_FALSE));
  PetscCall(MatISSetPreallocation(A, 3, NULL, 3, NULL));
  PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool) !(repmap || negmap))); /* I do not want to precompute the pattern */
  PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
  PetscCall(ISLocalToGlobalMappingGetSize(cmap, &ln));
  for (i = 0; i < lm; i++) {
    PetscScalar v[3];
    PetscInt    cols[3];

    cols[0] = (i - 1 + n) % n;
    cols[1] = i % n;
    cols[2] = (i + 1) % n;
    v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
    v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
    v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
    PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
    PetscCall(MatSetValuesLocal(A, 1, &i, 3, cols, v, ADD_VALUES));
  }
  PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

  /* activate tests for square matrices with same maps only */
  PetscCall(MatHasCongruentLayouts(A, &squaretest));
  if (squaretest && rmap != cmap) {
    PetscInt nr, nc;

    PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nr));
    PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nc));
    if (nr != nc) squaretest = PETSC_FALSE;
    else {
      const PetscInt *idxs1, *idxs2;

      PetscCall(ISLocalToGlobalMappingGetIndices(rmap, &idxs1));
      PetscCall(ISLocalToGlobalMappingGetIndices(cmap, &idxs2));
      PetscCall(PetscArraycmp(idxs1, idxs2, nr, &squaretest));
      PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1));
      PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2));
    }
    PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &squaretest, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
  }
  if (negmap && repmap) squaretest = PETSC_FALSE;

  /* test MatISGetLocalMat */
  PetscCall(MatISGetLocalMat(A, &B));
  PetscCall(MatGetType(B, &lmtype));

  /* test MatGetInfo */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetInfo\n"));
  PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
  PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "Process  %2d: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", PetscGlobalRank, (PetscInt)info.nz_used, (PetscInt)info.nz_allocated,
                                               (PetscInt)info.nz_unneeded, (PetscInt)info.assemblies, (PetscInt)info.mallocs));
  PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(MatGetInfo(A, MAT_GLOBAL_MAX, &info));
  PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "GlobalMax  : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated, (PetscInt)info.nz_unneeded,
                                   (PetscInt)info.assemblies, (PetscInt)info.mallocs));
  PetscCall(MatGetInfo(A, MAT_GLOBAL_SUM, &info));
  PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "GlobalSum  : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated, (PetscInt)info.nz_unneeded,
                                   (PetscInt)info.assemblies, (PetscInt)info.mallocs));

  /* test MatIsSymmetric */
  PetscCall(MatIsSymmetric(A, 0.0, &issymmetric));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIsSymmetric: %d\n", issymmetric));

  /* Create a MPIAIJ matrix, same as A */
  PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
  PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
  PetscCall(MatSetType(B, MATAIJ));
  PetscCall(MatSetFromOptions(B));
  PetscCall(MatSetLocalToGlobalMapping(B, rmap, cmap));
  PetscCall(MatMPIAIJSetPreallocation(B, 3, NULL, 3, NULL));
  PetscCall(MatMPIBAIJSetPreallocation(B, 1, 3, NULL, 3, NULL));
#if defined(PETSC_HAVE_HYPRE)
  PetscCall(MatHYPRESetPreallocation(B, 3, NULL, 3, NULL));
#endif
  PetscCall(MatISSetPreallocation(B, 3, NULL, 3, NULL));
  PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool) !(repmap || negmap))); /* I do not want to precompute the pattern */
  for (i = 0; i < lm; i++) {
    PetscScalar v[3];
    PetscInt    cols[3];

    cols[0] = (i - 1 + n) % n;
    cols[1] = i % n;
    cols[2] = (i + 1) % n;
    v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
    v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
    v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
    PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
    PetscCall(MatSetValuesLocal(B, 1, &i, 3, cols, v, ADD_VALUES));
  }
  PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

  /* test MatView */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView\n"));
  PetscCall(MatView(A, NULL));
  PetscCall(MatView(B, NULL));

  /* test CheckMat */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test CheckMat\n"));
  PetscCall(CheckMat(A, B, PETSC_FALSE, "CheckMat"));

  /* test MatDuplicate and MatAXPY */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDuplicate and MatAXPY\n"));
  PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
  PetscCall(CheckMat(A, A2, PETSC_FALSE, "MatDuplicate and MatAXPY"));

  /* test MatConvert */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_IS_XAIJ\n"));
  PetscCall(MatConvert(A2, MATAIJ, MAT_INITIAL_MATRIX, &B2));
  PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INITIAL_MATRIX"));
  PetscCall(MatConvert(A2, MATAIJ, MAT_REUSE_MATRIX, &B2));
  PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_REUSE_MATRIX"));
  PetscCall(MatConvert(A2, MATAIJ, MAT_INPLACE_MATRIX, &A2));
  PetscCall(CheckMat(B, A2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INPLACE_MATRIX"));
  PetscCall(MatDestroy(&A2));
  PetscCall(MatDestroy(&B2));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_XAIJ_IS\n"));
  PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
  PetscCall(MatConvert(B2, MATIS, MAT_INITIAL_MATRIX, &A2));
  PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INITIAL_MATRIX"));
  PetscCall(MatConvert(B2, MATIS, MAT_REUSE_MATRIX, &A2));
  PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_REUSE_MATRIX"));
  PetscCall(MatConvert(B2, MATIS, MAT_INPLACE_MATRIX, &B2));
  PetscCall(CheckMat(A, B2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INPLACE_MATRIX"));
  PetscCall(MatDestroy(&A2));
  PetscCall(MatDestroy(&B2));
  PetscCall(PetscStrcmp(lmtype, MATSEQAIJ, &isaij));
  if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
    PetscInt               ri, ci, rr[3] = {0, 1, 0}, cr[4] = {1, 2, 0, 1}, rk[3] = {0, 2, 1}, ck[4] = {1, 0, 3, 2};
    ISLocalToGlobalMapping tcmap, trmap;

    for (ri = 0; ri < 2; ri++) {
      PetscInt *r;

      r = (PetscInt *)(ri == 0 ? rr : rk);
      for (ci = 0; ci < 2; ci++) {
        PetscInt *c, rb, cb;

        c = (PetscInt *)(ci == 0 ? cr : ck);
        for (rb = 1; rb < 4; rb++) {
          PetscCall(ISCreateBlock(PETSC_COMM_SELF, rb, 3, r, PETSC_COPY_VALUES, &is));
          PetscCall(ISLocalToGlobalMappingCreateIS(is, &trmap));
          PetscCall(ISDestroy(&is));
          for (cb = 1; cb < 4; cb++) {
            Mat  T, lT, T2;
            char testname[256];

            PetscCall(PetscSNPrintf(testname, sizeof(testname), "MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")", ri, ci, rb, cb));
            PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test %s\n", testname));

            PetscCall(ISCreateBlock(PETSC_COMM_SELF, cb, 4, c, PETSC_COPY_VALUES, &is));
            PetscCall(ISLocalToGlobalMappingCreateIS(is, &tcmap));
            PetscCall(ISDestroy(&is));

            PetscCall(MatCreate(PETSC_COMM_SELF, &T));
            PetscCall(MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, rb * 3, cb * 4));
            PetscCall(MatSetType(T, MATIS));
            PetscCall(MatSetLocalToGlobalMapping(T, trmap, tcmap));
            PetscCall(ISLocalToGlobalMappingDestroy(&tcmap));
            PetscCall(MatISGetLocalMat(T, &lT));
            PetscCall(MatSetType(lT, MATSEQAIJ));
            PetscCall(MatSeqAIJSetPreallocation(lT, cb * 4, NULL));
            PetscCall(MatSetRandom(lT, NULL));
            PetscCall(MatConvert(lT, lmtype, MAT_INPLACE_MATRIX, &lT));
            PetscCall(MatISRestoreLocalMat(T, &lT));
            PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
            PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));

            PetscCall(MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &T2));
            PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INITIAL_MATRIX"));
            PetscCall(MatConvert(T, MATAIJ, MAT_REUSE_MATRIX, &T2));
            PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_REUSE_MATRIX"));
            PetscCall(MatDestroy(&T2));
            PetscCall(MatDuplicate(T, MAT_COPY_VALUES, &T2));
            PetscCall(MatConvert(T2, MATAIJ, MAT_INPLACE_MATRIX, &T2));
            PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INPLACE_MATRIX"));
            PetscCall(MatDestroy(&T));
            PetscCall(MatDestroy(&T2));
          }
          PetscCall(ISLocalToGlobalMappingDestroy(&trmap));
        }
      }
    }
  }

  /* test MatDiagonalScale */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalScale\n"));
  PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
  PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
  PetscCall(MatCreateVecs(A, &x, &y));
  PetscCall(VecSetRandom(x, NULL));
  if (issymmetric) {
    PetscCall(VecCopy(x, y));
  } else {
    PetscCall(VecSetRandom(y, NULL));
    PetscCall(VecScale(y, 8.));
  }
  PetscCall(MatDiagonalScale(A2, y, x));
  PetscCall(MatDiagonalScale(B2, y, x));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalScale"));
  PetscCall(MatDestroy(&A2));
  PetscCall(MatDestroy(&B2));
  PetscCall(VecDestroy(&x));
  PetscCall(VecDestroy(&y));

  /* test MatPtAP (A IS and B AIJ) */
  if (isaij && m == n) {
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatPtAP\n"));
    /* There's a bug in MatCreateSubMatrices_MPIAIJ I cannot figure out */
    if (!allow_repeated || !repmap || size == 1) {
      PetscCall(MatISStoreL2L(A, PETSC_TRUE));
      PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &A2));
      PetscCall(MatPtAP(B, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B2));
      PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_INITIAL_MATRIX"));
      PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &A2));
      PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_REUSE_MATRIX"));
      PetscCall(MatDestroy(&A2));
      PetscCall(MatDestroy(&B2));
    }
  }

  /* test MatGetLocalSubMatrix */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetLocalSubMatrix\n"));
  PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A2));
  PetscCall(ISCreateStride(PETSC_COMM_SELF, lm / 2 + lm % 2, 0, 2, &reven));
  PetscCall(ISComplement(reven, 0, lm, &rodd));
  PetscCall(ISCreateStride(PETSC_COMM_SELF, ln / 2 + ln % 2, 0, 2, &ceven));
  PetscCall(ISComplement(ceven, 0, ln, &codd));
  PetscCall(MatGetLocalSubMatrix(A2, reven, ceven, &Aee));
  PetscCall(MatGetLocalSubMatrix(A2, reven, codd, &Aeo));
  PetscCall(MatGetLocalSubMatrix(A2, rodd, ceven, &Aoe));
  PetscCall(MatGetLocalSubMatrix(A2, rodd, codd, &Aoo));
  for (i = 0; i < lm; i++) {
    PetscInt    j, je, jo, colse[3], colso[3];
    PetscScalar ve[3], vo[3];
    PetscScalar v[3];
    PetscInt    cols[3];
    PetscInt    row = i / 2;

    cols[0] = (i - 1 + n) % n;
    cols[1] = i % n;
    cols[2] = (i + 1) % n;
    v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
    v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
    v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
    PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
    for (j = 0, je = 0, jo = 0; j < 3; j++) {
      if (cols[j] % 2) {
        vo[jo]      = v[j];
        colso[jo++] = cols[j] / 2;
      } else {
        ve[je]      = v[j];
        colse[je++] = cols[j] / 2;
      }
    }
    if (i % 2) {
      PetscCall(MatSetValuesLocal(Aoe, 1, &row, je, colse, ve, ADD_VALUES));
      PetscCall(MatSetValuesLocal(Aoo, 1, &row, jo, colso, vo, ADD_VALUES));
    } else {
      PetscCall(MatSetValuesLocal(Aee, 1, &row, je, colse, ve, ADD_VALUES));
      PetscCall(MatSetValuesLocal(Aeo, 1, &row, jo, colso, vo, ADD_VALUES));
    }
  }
  PetscCall(MatRestoreLocalSubMatrix(A2, reven, ceven, &Aee));
  PetscCall(MatRestoreLocalSubMatrix(A2, reven, codd, &Aeo));
  PetscCall(MatRestoreLocalSubMatrix(A2, rodd, ceven, &Aoe));
  PetscCall(MatRestoreLocalSubMatrix(A2, rodd, codd, &Aoo));
  PetscCall(ISDestroy(&reven));
  PetscCall(ISDestroy(&ceven));
  PetscCall(ISDestroy(&rodd));
  PetscCall(ISDestroy(&codd));
  PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAXPY(A2, -1., A, SAME_NONZERO_PATTERN));
  PetscCall(CheckMat(A2, NULL, PETSC_FALSE, "MatGetLocalSubMatrix"));
  PetscCall(MatDestroy(&A2));

  /* test MatConvert_Nest_IS */
  testT = PETSC_FALSE;
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_trans", &testT, NULL));

  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_Nest_IS\n"));
  nr = 2;
  nc = 2;
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-nr", &nr, NULL));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
  if (testT) {
    PetscCall(MatGetOwnershipRange(A, &cst, &cen));
    PetscCall(MatGetOwnershipRangeColumn(A, &rst, &ren));
  } else {
    PetscCall(MatGetOwnershipRange(A, &rst, &ren));
    PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen));
  }
  PetscCall(PetscMalloc3(nr, &rows, nc, &cols, 2 * nr * nc, &mats));
  for (i = 0; i < nr * nc; i++) {
    if (testT) {
      PetscCall(MatCreateTranspose(A, &mats[i]));
      PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &mats[i + nr * nc]));
    } else {
      PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &mats[i]));
      PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &mats[i + nr * nc]));
    }
  }
  for (i = 0; i < nr; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, ren - rst, i + rst, nr, &rows[i]));
  for (i = 0; i < nc; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, cen - cst, i + cst, nc, &cols[i]));
  PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats, &A2));
  PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats + nr * nc, &B2));
  for (i = 0; i < nr; i++) PetscCall(ISDestroy(&rows[i]));
  for (i = 0; i < nc; i++) PetscCall(ISDestroy(&cols[i]));
  for (i = 0; i < 2 * nr * nc; i++) PetscCall(MatDestroy(&mats[i]));
  PetscCall(PetscFree3(rows, cols, mats));
  PetscCall(MatConvert(B2, MATAIJ, MAT_INITIAL_MATRIX, &T));
  PetscCall(MatDestroy(&B2));
  PetscCall(MatConvert(A2, MATIS, MAT_INITIAL_MATRIX, &B2));
  PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INITIAL_MATRIX"));
  PetscCall(MatConvert(A2, MATIS, MAT_REUSE_MATRIX, &B2));
  PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_REUSE_MATRIX"));
  PetscCall(MatDestroy(&B2));
  PetscCall(MatConvert(A2, MATIS, MAT_INPLACE_MATRIX, &A2));
  PetscCall(CheckMat(A2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INPLACE_MATRIX"));
  PetscCall(MatDestroy(&T));
  PetscCall(MatDestroy(&A2));

  /* test MatCreateSubMatrix */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrix\n"));
  if (rank == 0) {
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 1, 1, &is));
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, &is2));
  } else if (rank == 1) {
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
    if (n > 3) {
      PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 3, 1, &is2));
    } else {
      PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
    }
  } else if (rank == 2 && n > 4) {
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, n - 4, 4, 1, &is2));
  } else {
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
  }
  PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &A2));
  PetscCall(MatCreateSubMatrix(B, is, is, MAT_INITIAL_MATRIX, &B2));
  PetscCall(CheckMat(A2, B2, PETSC_TRUE, "first MatCreateSubMatrix"));

  PetscCall(MatCreateSubMatrix(A, is, is, MAT_REUSE_MATRIX, &A2));
  PetscCall(MatCreateSubMatrix(B, is, is, MAT_REUSE_MATRIX, &B2));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse MatCreateSubMatrix"));
  PetscCall(MatDestroy(&A2));
  PetscCall(MatDestroy(&B2));

  if (!issymmetric) {
    PetscCall(MatCreateSubMatrix(A, is, is2, MAT_INITIAL_MATRIX, &A2));
    PetscCall(MatCreateSubMatrix(B, is, is2, MAT_INITIAL_MATRIX, &B2));
    PetscCall(MatCreateSubMatrix(A, is, is2, MAT_REUSE_MATRIX, &A2));
    PetscCall(MatCreateSubMatrix(B, is, is2, MAT_REUSE_MATRIX, &B2));
    PetscCall(CheckMat(A2, B2, PETSC_FALSE, "second MatCreateSubMatrix"));
  }

  PetscCall(MatDestroy(&A2));
  PetscCall(MatDestroy(&B2));
  PetscCall(ISDestroy(&is));
  PetscCall(ISDestroy(&is2));

  /* test MatCreateSubMatrices */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrices\n"));
  PetscCall(MatGetLayouts(A, &rlayout, &clayout));
  PetscCall(PetscLayoutGetRanges(rlayout, &rrange));
  PetscCall(PetscLayoutGetRanges(clayout, &crange));
  lrank = (size + rank - 1) % size;
  rrank = (rank + 1) % size;
  PetscCall(ISCreateStride(PETSC_COMM_SELF, (rrange[lrank + 1] - rrange[lrank]), rrange[lrank], 1, &irow[0]));
  PetscCall(ISCreateStride(PETSC_COMM_SELF, (crange[rrank + 1] - crange[rrank]), crange[rrank], 1, &icol[0]));
  PetscCall(ISCreateStride(PETSC_COMM_SELF, (rrange[rrank + 1] - rrange[rrank]), rrange[rrank], 1, &irow[1]));
  PetscCall(ISCreateStride(PETSC_COMM_SELF, (crange[lrank + 1] - crange[lrank]), crange[lrank], 1, &icol[1]));
  PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_INITIAL_MATRIX, &Asub));
  PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_INITIAL_MATRIX, &Bsub));
  PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
  PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
  PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_REUSE_MATRIX, &Asub));
  PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_REUSE_MATRIX, &Bsub));
  PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
  PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
  PetscCall(MatDestroySubMatrices(2, &Asub));
  PetscCall(MatDestroySubMatrices(2, &Bsub));
  PetscCall(ISDestroy(&irow[0]));
  PetscCall(ISDestroy(&irow[1]));
  PetscCall(ISDestroy(&icol[0]));
  PetscCall(ISDestroy(&icol[1]));

  /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
  if (size > 1) {
    if (rank == 0) {
      PetscInt st, len;

      st  = (m + 1) / 2;
      len = PetscMin(m / 2, PetscMax(m - (m + 1) / 2 - 1, 0));
      PetscCall(ISCreateStride(PETSC_COMM_WORLD, len, st, 1, &is));
    } else {
      PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
    }
  } else {
    PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
  }

  if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */
    PetscInt *idx0, *idx1, n0, n1;
    IS        Ais[2], Bis[2];

    /* test MatDiagonalSet */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalSet\n"));
    PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
    PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
    PetscCall(MatCreateVecs(A, NULL, &x));
    PetscCall(VecSetRandom(x, NULL));
    PetscCall(MatDiagonalSet(A2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
    PetscCall(MatDiagonalSet(B2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
    PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalSet"));
    PetscCall(VecDestroy(&x));
    PetscCall(MatDestroy(&A2));
    PetscCall(MatDestroy(&B2));

    /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatShift\n"));
    PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
    PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
    PetscCall(MatShift(A2, 2.0));
    PetscCall(MatShift(B2, 2.0));
    PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatShift"));
    PetscCall(MatDestroy(&A2));
    PetscCall(MatDestroy(&B2));

    /* nonzero diag value is supported for square matrices only */
    PetscCall(TestMatZeroRows(A, B, PETSC_TRUE, is, diag));

    /* test MatIncreaseOverlap */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIncreaseOverlap\n"));
    PetscCall(MatGetOwnershipRange(A, &rst, &ren));
    n0 = (ren - rst) / 2;
    n1 = (ren - rst) / 3;
    PetscCall(PetscMalloc1(n0, &idx0));
    PetscCall(PetscMalloc1(n1, &idx1));
    for (i = 0; i < n0; i++) idx0[i] = ren - i - 1;
    for (i = 0; i < n1; i++) idx1[i] = rst + i;
    PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_OWN_POINTER, &Ais[0]));
    PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_OWN_POINTER, &Ais[1]));
    PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_COPY_VALUES, &Bis[0]));
    PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_COPY_VALUES, &Bis[1]));
    PetscCall(MatIncreaseOverlap(A, 2, Ais, 3));
    PetscCall(MatIncreaseOverlap(B, 2, Bis, 3));
    /* Non deterministic output! */
    PetscCall(ISSort(Ais[0]));
    PetscCall(ISSort(Ais[1]));
    PetscCall(ISSort(Bis[0]));
    PetscCall(ISSort(Bis[1]));
    PetscCall(ISView(Ais[0], NULL));
    PetscCall(ISView(Bis[0], NULL));
    PetscCall(ISView(Ais[1], NULL));
    PetscCall(ISView(Bis[1], NULL));
    PetscCall(MatCreateSubMatrices(A, 2, Ais, Ais, MAT_INITIAL_MATRIX, &Asub));
    PetscCall(MatCreateSubMatrices(B, 2, Bis, Ais, MAT_INITIAL_MATRIX, &Bsub));
    PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatIncreaseOverlap[0]"));
    PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatIncreaseOverlap[1]"));
    PetscCall(MatDestroySubMatrices(2, &Asub));
    PetscCall(MatDestroySubMatrices(2, &Bsub));
    PetscCall(ISDestroy(&Ais[0]));
    PetscCall(ISDestroy(&Ais[1]));
    PetscCall(ISDestroy(&Bis[0]));
    PetscCall(ISDestroy(&Bis[1]));
  }
  PetscCall(TestMatZeroRows(A, B, squaretest, is, 0.0));
  PetscCall(ISDestroy(&is));

  /* test MatTranspose */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatTranspose\n"));
  PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
  PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &B2));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "initial matrix MatTranspose"));

  PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &A2));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (not in place) MatTranspose"));
  PetscCall(MatDestroy(&A2));

  PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
  PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (in place) MatTranspose"));
  PetscCall(MatDestroy(&A2));

  PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (different type) MatTranspose"));
  PetscCall(MatDestroy(&A2));
  PetscCall(MatDestroy(&B2));

  /* test MatISFixLocalEmpty */
  if (isaij) {
    PetscInt r[2];

    r[0] = 0;
    r[1] = PetscMin(m, n) - 1;
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISFixLocalEmpty\n"));
    PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));

    PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
    PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
    PetscCall(CheckMat(A2, B, PETSC_FALSE, "MatISFixLocalEmpty (null)"));

    PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
    PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
    PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
    PetscCall(MatZeroRows(B2, 2, r, 0.0, NULL, NULL));
    PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
    PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
    PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
    PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows)"));
    PetscCall(MatDestroy(&A2));

    PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
    PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
    PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
    PetscCall(MatTranspose(B2, MAT_INPLACE_MATRIX, &B2));
    PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
    PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
    PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
    PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
    PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (cols)"));

    PetscCall(MatDestroy(&A2));
    PetscCall(MatDestroy(&B2));

    if (squaretest) {
      PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
      PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
      PetscCall(MatZeroRowsColumns(A2, 2, r, 0.0, NULL, NULL));
      PetscCall(MatZeroRowsColumns(B2, 2, r, 0.0, NULL, NULL));
      PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
      PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
      PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
      PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
      PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
      PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows+cols)"));
      PetscCall(MatDestroy(&A2));
      PetscCall(MatDestroy(&B2));
    }
  }

  /* test MatInvertBlockDiagonal
       special cases for block-diagonal matrices */
  if (m == n) {
    ISLocalToGlobalMapping map;
    Mat                    Abd, Bbd;
    IS                     is, bis;
    const PetscScalar     *isbd, *aijbd;
    PetscScalar           *vals;
    const PetscInt        *sts, *idxs;
    PetscInt              *idxs2, diff, perm, nl, bs, st, en, in;
    PetscBool              ok;

    for (diff = 0; diff < 3; diff++) {
      for (perm = 0; perm < 3; perm++) {
        for (bs = 1; bs < 4; bs++) {
          PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", n, diff, perm, bs));
          PetscCall(PetscMalloc1(bs * bs, &vals));
          PetscCall(MatGetOwnershipRanges(A, &sts));
          switch (diff) {
          case 1: /* inverted layout by processes */
            in = 1;
            st = sts[size - rank - 1];
            en = sts[size - rank];
            nl = en - st;
            break;
          case 2: /* round-robin layout */
            in = size;
            st = rank;
            nl = n / size;
            if (rank < n % size) nl++;
            break;
          default: /* same layout */
            in = 1;
            st = sts[rank];
            en = sts[rank + 1];
            nl = en - st;
            break;
          }
          PetscCall(ISCreateStride(PETSC_COMM_WORLD, nl, st, in, &is));
          PetscCall(ISGetLocalSize(is, &nl));
          PetscCall(ISGetIndices(is, &idxs));
          PetscCall(PetscMalloc1(nl, &idxs2));
          for (i = 0; i < nl; i++) {
            switch (perm) { /* invert some of the indices */
            case 2:
              idxs2[i] = rank % 2 ? idxs[i] : idxs[nl - i - 1];
              break;
            case 1:
              idxs2[i] = rank % 2 ? idxs[nl - i - 1] : idxs[i];
              break;
            default:
              idxs2[i] = idxs[i];
              break;
            }
          }
          PetscCall(ISRestoreIndices(is, &idxs));
          PetscCall(ISCreateBlock(PETSC_COMM_WORLD, bs, nl, idxs2, PETSC_OWN_POINTER, &bis));
          PetscCall(ISLocalToGlobalMappingCreateIS(bis, &map));
          PetscCall(MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, bs * n, bs * n, map, map, &Abd));
          PetscCall(ISLocalToGlobalMappingDestroy(&map));
          PetscCall(MatISSetPreallocation(Abd, bs, NULL, 0, NULL));
          for (i = 0; i < nl; i++) {
            PetscInt b1, b2;

            for (b1 = 0; b1 < bs; b1++) {
              for (b2 = 0; b2 < bs; b2++) vals[b1 * bs + b2] = i * bs * bs + b1 * bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
            }
            PetscCall(MatSetValuesBlockedLocal(Abd, 1, &i, 1, &i, vals, INSERT_VALUES));
          }
          PetscCall(MatAssemblyBegin(Abd, MAT_FINAL_ASSEMBLY));
          PetscCall(MatAssemblyEnd(Abd, MAT_FINAL_ASSEMBLY));
          PetscCall(MatConvert(Abd, MATAIJ, MAT_INITIAL_MATRIX, &Bbd));
          PetscCall(MatInvertBlockDiagonal(Abd, &isbd));
          PetscCall(MatInvertBlockDiagonal(Bbd, &aijbd));
          PetscCall(MatGetLocalSize(Bbd, &nl, NULL));
          ok = PETSC_TRUE;
          for (i = 0; i < nl / bs; i++) {
            PetscInt b1, b2;

            for (b1 = 0; b1 < bs; b1++) {
              for (b2 = 0; b2 < bs; b2++) {
                if (PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2] - aijbd[i * bs * bs + b1 * bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
                if (!ok) {
                  PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] ERROR block %" PetscInt_FMT ", entry %" PetscInt_FMT " %" PetscInt_FMT ": %g %g\n", rank, i, b1, b2, (double)PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2]), (double)PetscAbsScalar(aijbd[i * bs * bs + b1 * bs + b2])));
                  break;
                }
              }
              if (!ok) break;
            }
            if (!ok) break;
          }
          PetscCall(MatDestroy(&Abd));
          PetscCall(MatDestroy(&Bbd));
          PetscCall(PetscFree(vals));
          PetscCall(ISDestroy(&is));
          PetscCall(ISDestroy(&bis));
        }
      }
    }
  }

  /* test MatGetDiagonalBlock */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetDiagonalBlock\n"));
  PetscCall(MatGetDiagonalBlock(A, &A2));
  PetscCall(MatGetDiagonalBlock(B, &B2));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));
  PetscCall(MatScale(A, 2.0));
  PetscCall(MatScale(B, 2.0));
  PetscCall(MatGetDiagonalBlock(A, &A2));
  PetscCall(MatGetDiagonalBlock(B, &B2));
  PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));

  /* test MatISSetAllowRepeated on a MATIS */
  PetscCall(MatISSetAllowRepeated(A, allow_repeated));
  if (allow_repeated) { /* original MATIS maybe with repeated entries, test assembling of local matrices */
    Mat lA, lA2;

    for (PetscInt i = 0; i < 1; i++) { /* TODO: make MatScatter inherit from MATSHELL and support MatProducts */
      PetscBool usemult = PETSC_FALSE;

      PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
      if (i) {
        Mat tA;

        usemult = PETSC_TRUE;
        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries and local shell matrices\n"));
        PetscCall(MatISGetLocalMat(A2, &lA2));
        PetscCall(MatConvert(lA2, MATSHELL, MAT_INITIAL_MATRIX, &tA));
        PetscCall(MatISRestoreLocalMat(A2, &lA2));
        PetscCall(MatISSetLocalMat(A2, tA));
        PetscCall(MatDestroy(&tA));
      } else {
        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries\n"));
      }
      PetscCall(MatISSetAllowRepeated(A2, PETSC_FALSE));
      PetscCall(MatISGetLocalMat(A, &lA));
      PetscCall(MatISGetLocalMat(A2, &lA2));
      if (!repmap) PetscCall(CheckMat(lA, lA2, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
      PetscCall(MatISRestoreLocalMat(A, &lA));
      PetscCall(MatISRestoreLocalMat(A2, &lA2));
      if (repmap) PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with repeated entries"));
      else PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
      PetscCall(MatDestroy(&A2));
    }
  } else { /* original matis with non-repeated entries, this should only recreate the local matrices */
    Mat       lA;
    PetscBool flg;

    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(true) with non repeated entries\n"));
    PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
    PetscCall(MatISSetAllowRepeated(A2, PETSC_TRUE));
    PetscCall(MatISGetLocalMat(A2, &lA));
    PetscCall(MatAssembled(lA, &flg));
    PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local mat should be unassembled");
    PetscCall(MatISRestoreLocalMat(A2, &lA));
    PetscCall(MatDestroy(&A2));
  }

  /* free testing matrices */
  PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
  PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
  PetscCall(MatDestroy(&A));
  PetscCall(MatDestroy(&B));
  PetscCall(PetscFinalize());
  return 0;
}

PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char *func)
{
  Mat       Bcheck;
  PetscReal error;

  PetscFunctionBeginUser;
  if (!usemult && B) {
    PetscBool hasnorm;

    PetscCall(MatHasOperation(B, MATOP_NORM, &hasnorm));
    if (!hasnorm) usemult = PETSC_TRUE;
  }
  if (!usemult) {
    if (B) {
      MatType Btype;

      PetscCall(MatGetType(B, &Btype));
      PetscCall(MatConvert(A, Btype, MAT_INITIAL_MATRIX, &Bcheck));
    } else {
      PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
    }
    if (B) { /* if B is present, subtract it */
      PetscCall(MatAXPY(Bcheck, -1., B, DIFFERENT_NONZERO_PATTERN));
    }
    PetscCall(MatNorm(Bcheck, NORM_INFINITY, &error));
    if (error > PETSC_SQRT_MACHINE_EPSILON) {
      ISLocalToGlobalMapping rl2g, cl2g;

      PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Bcheck"));
      PetscCall(MatView(Bcheck, NULL));
      if (B) {
        PetscCall(PetscObjectSetName((PetscObject)B, "B"));
        PetscCall(MatView(B, NULL));
        PetscCall(MatDestroy(&Bcheck));
        PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
        PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Assembled A"));
        PetscCall(MatView(Bcheck, NULL));
      }
      PetscCall(MatDestroy(&Bcheck));
      PetscCall(PetscObjectSetName((PetscObject)A, "A"));
      PetscCall(MatView(A, NULL));
      PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
      if (rl2g) PetscCall(ISLocalToGlobalMappingView(rl2g, NULL));
      if (cl2g) PetscCall(ISLocalToGlobalMappingView(cl2g, NULL));
      SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "ERROR ON %s: %g", func, (double)error);
    }
    PetscCall(MatDestroy(&Bcheck));
  } else {
    PetscBool ok, okt;

    PetscCall(MatMultEqual(A, B, 3, &ok));
    PetscCall(MatMultTransposeEqual(A, B, 3, &okt));
    PetscCheck(ok && okt, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR ON %s: mult ok ?  %d, multtranspose ok ? %d", func, ok, okt);
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag)
{
  Mat                    B, Bcheck, B2 = NULL, lB;
  Vec                    x = NULL, b = NULL, b2 = NULL;
  ISLocalToGlobalMapping l2gr, l2gc;
  PetscReal              error;
  char                   diagstr[16];
  const PetscInt        *idxs;
  PetscInt               rst, ren, i, n, N, d;
  PetscMPIInt            rank;
  PetscBool              miss, haszerorows;

  PetscFunctionBeginUser;
  if (diag == 0.) {
    PetscCall(PetscStrncpy(diagstr, "zero", sizeof(diagstr)));
  } else {
    PetscCall(PetscStrncpy(diagstr, "nonzero", sizeof(diagstr)));
  }
  PetscCall(ISView(is, NULL));
  PetscCall(MatGetLocalToGlobalMapping(A, &l2gr, &l2gc));
  /* tests MatDuplicate and MatCopy */
  if (diag == 0.) {
    PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
  } else {
    PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
    PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
  }
  PetscCall(MatISGetLocalMat(B, &lB));
  PetscCall(MatHasOperation(lB, MATOP_ZERO_ROWS, &haszerorows));
  if (squaretest && haszerorows) {
    PetscCall(MatCreateVecs(B, &x, &b));
    PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
    PetscCall(VecSetLocalToGlobalMapping(b, l2gr));
    PetscCall(VecSetLocalToGlobalMapping(x, l2gc));
    PetscCall(VecSetRandom(x, NULL));
    PetscCall(VecSetRandom(b, NULL));
    /* mimic b[is] = x[is] */
    PetscCall(VecDuplicate(b, &b2));
    PetscCall(VecSetLocalToGlobalMapping(b2, l2gr));
    PetscCall(VecCopy(b, b2));
    PetscCall(ISGetLocalSize(is, &n));
    PetscCall(ISGetIndices(is, &idxs));
    PetscCall(VecGetSize(x, &N));
    for (i = 0; i < n; i++) {
      if (0 <= idxs[i] && idxs[i] < N) {
        PetscCall(VecSetValue(b2, idxs[i], diag, INSERT_VALUES));
        PetscCall(VecSetValue(x, idxs[i], 1., INSERT_VALUES));
      }
    }
    PetscCall(VecAssemblyBegin(b2));
    PetscCall(VecAssemblyEnd(b2));
    PetscCall(VecAssemblyBegin(x));
    PetscCall(VecAssemblyEnd(x));
    PetscCall(ISRestoreIndices(is, &idxs));
    /*  test ZeroRows on MATIS */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
    PetscCall(MatZeroRowsIS(B, is, diag, x, b));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumns (diag %s)\n", diagstr));
    PetscCall(MatZeroRowsColumnsIS(B2, is, diag, NULL, NULL));
  } else if (haszerorows) {
    /*  test ZeroRows on MATIS */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
    PetscCall(MatZeroRowsIS(B, is, diag, NULL, NULL));
    b = b2 = x = NULL;
  } else {
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Skipping MatZeroRows (diag %s)\n", diagstr));
    b = b2 = x = NULL;
  }

  if (squaretest && haszerorows) {
    PetscCall(VecAXPY(b2, -1., b));
    PetscCall(VecNorm(b2, NORM_INFINITY, &error));
    PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR IN ZEROROWS ON B %g (diag %s)", (double)error, diagstr);
  }

  /* test MatMissingDiagonal */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatMissingDiagonal\n"));
  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
  PetscCall(MatMissingDiagonal(B, &miss, &d));
  PetscCall(MatGetOwnershipRange(B, &rst, &ren));
  PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%" PetscInt_FMT ",%" PetscInt_FMT ") Missing %d, row %" PetscInt_FMT " (diag %s)\n", rank, rst, ren, (int)miss, d, diagstr));
  PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

  PetscCall(VecDestroy(&x));
  PetscCall(VecDestroy(&b));
  PetscCall(VecDestroy(&b2));

  /* check the result of ZeroRows with that from MPIAIJ routines
     assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
  if (haszerorows) {
    PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck));
    PetscCall(MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
    PetscCall(MatZeroRowsIS(Bcheck, is, diag, NULL, NULL));
    PetscCall(CheckMat(B, Bcheck, PETSC_FALSE, "Zerorows"));
    PetscCall(MatDestroy(&Bcheck));
  }
  PetscCall(MatDestroy(&B));

  if (B2) { /* test MatZeroRowsColumns */
    PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &B));
    PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
    PetscCall(MatZeroRowsColumnsIS(B, is, diag, NULL, NULL));
    PetscCall(CheckMat(B2, B, PETSC_FALSE, "MatZeroRowsColumns"));
    PetscCall(MatDestroy(&B));
    PetscCall(MatDestroy(&B2));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*TEST

   test:
      args: -test_trans

   test:
      suffix: 2
      nsize: 4
      args: -mat_is_convert_local_nest -nr 3 -nc 4

   test:
      suffix: 3
      nsize: 5
      args: -m 11 -n 10 -mat_is_convert_local_nest -nr 2 -nc 1

   test:
      suffix: 4
      nsize: 6
      args: -m 9 -n 12 -test_trans -nr 2 -nc 7

   test:
      suffix: 5
      nsize: 6
      args: -m 12 -n 12 -test_trans -nr 3 -nc 1

   test:
      suffix: 6
      args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap

   test:
      suffix: 7
      args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap

   test:
      suffix: 8
      args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap

   test:
      suffix: 9
      nsize: 5
      args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap

   test:
      suffix: 10
      nsize: 5
      args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap

   test:
      suffix: vscat_default
      nsize: 5
      args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
      output_file: output/ex23_11.out

   test:
      suffix: 12
      nsize: 3
      args: -m 12 -n 12 -symmetric -mat_is_localmat_type sbaij -test_trans -nr 2 -nc 3

   testset:
      output_file: output/ex23_13.out
      nsize: 3
      args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
      filter: grep -v "type:"
      test:
        suffix: baij
        args: -mat_is_localmat_type baij
      test:
        requires: viennacl
        suffix: viennacl
        args: -mat_is_localmat_type aijviennacl
      test:
        requires: cuda
        suffix: cusparse
        args: -mat_is_localmat_type aijcusparse

   test:
      suffix: negrep
      nsize: {{1 3}separate output}
      args: -m {{5 7}separate output} -n {{5 7}separate output} -test_trans -nr 2 -nc 3 -negmap {{0 1}separate output} -repmap {{0 1}separate output} -permmap -diffmap {{0 1}separate output} -allow_repeated 0

   test:
      suffix: negrep_allowrep
      nsize: {{1 3}separate output}
      args: -m {{5 7}separate output} -n {{5 7}separate output} -test_trans -nr 2 -nc 3 -negmap {{0 1}separate output} -repmap {{0 1}separate output} -permmap -diffmap {{0 1}separate output} -allow_repeated

TEST*/