static char help[] = "Creates a matrix using 9 pt stencil, and uses it to test MatIncreaseOverlap (needed for additive Schwarz preconditioner). \n\
  -m <size>       : problem size\n\
  -x1, -x2 <size> : no of subdomains in x and y directions\n\n";

#include <petscksp.h>

static PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
{
  PetscFunctionBeginUser;
  Ke[0]  = H / 6.0;
  Ke[1]  = -.125 * H;
  Ke[2]  = H / 12.0;
  Ke[3]  = -.125 * H;
  Ke[4]  = -.125 * H;
  Ke[5]  = H / 6.0;
  Ke[6]  = -.125 * H;
  Ke[7]  = H / 12.0;
  Ke[8]  = H / 12.0;
  Ke[9]  = -.125 * H;
  Ke[10] = H / 6.0;
  Ke[11] = -.125 * H;
  Ke[12] = -.125 * H;
  Ke[13] = H / 12.0;
  Ke[14] = -.125 * H;
  Ke[15] = H / 6.0;
  PetscFunctionReturn(PETSC_SUCCESS);
}

#if 0
// unused
static PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
{
  PetscFunctionBeginUser;
  r[0] = 0.;
  r[1] = 0.;
  r[2] = 0.;
  r[3] = 0.0;
  PetscFunctionReturn(PETSC_SUCCESS);
}
#endif

int main(int argc, char **args)
{
  Mat         C;
  PetscInt    i, m = 2, N, M, idx[4], Nsub1, Nsub2, ol = 1, x1, x2;
  PetscScalar Ke[16];
  PetscReal   h;
  IS         *is1, *is2, *islocal1, *islocal2;
  PetscBool   flg;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &args, NULL, help));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
  N  = (m + 1) * (m + 1); /* dimension of matrix */
  M  = m * m;             /* number of elements */
  h  = 1.0 / m;           /* mesh width */
  x1 = (m + 1) / 2;
  x2 = x1;
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-x1", &x1, NULL));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-x2", &x2, NULL));
  /* create stiffness matrix */
  PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 9, NULL, &C));

  /* forms the element stiffness for the Laplacian */
  PetscCall(FormElementStiffness(h * h, Ke));
  for (i = 0; i < M; i++) {
    /* node numbers for the four corners of element */
    idx[0] = (m + 1) * (i / m) + (i % m);
    idx[1] = idx[0] + 1;
    idx[2] = idx[1] + m + 1;
    idx[3] = idx[2] - 1;
    PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES));
  }
  PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

  for (ol = 0; ol < m + 2; ++ol) {
    PetscCall(PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, 0, &Nsub1, &is1, &islocal1));
    PetscCall(MatIncreaseOverlap(C, Nsub1, is1, ol));
    PetscCall(PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, ol, &Nsub2, &is2, &islocal2));

    PetscCall(PetscPrintf(PETSC_COMM_SELF, "flg == 1 => both index sets are same\n"));
    if (Nsub1 != Nsub2) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: No of index sets don't match\n"));

    for (i = 0; i < Nsub1; ++i) {
      PetscCall(ISEqual(is1[i], is2[i], &flg));
      PetscCall(PetscPrintf(PETSC_COMM_SELF, "i =  %" PetscInt_FMT ",flg = %d \n", i, (int)flg));
    }
    for (i = 0; i < Nsub1; ++i) PetscCall(ISDestroy(&is1[i]));
    for (i = 0; i < Nsub2; ++i) PetscCall(ISDestroy(&is2[i]));
    for (i = 0; i < Nsub1; ++i) PetscCall(ISDestroy(&islocal1[i]));
    for (i = 0; i < Nsub2; ++i) PetscCall(ISDestroy(&islocal2[i]));

    PetscCall(PetscFree(is1));
    PetscCall(PetscFree(is2));
    PetscCall(PetscFree(islocal1));
    PetscCall(PetscFree(islocal2));
  }
  PetscCall(MatDestroy(&C));
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

   test:
      args: -m 7

TEST*/