1: static char help[] = "Tests MatISFixLocalEmpty() with a blocked local to global mapping.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat                    A, B;
  8:   ISLocalToGlobalMapping map, rl2g, cl2g;
  9:   PetscScalar            v;
 10:   PetscInt              *idxs, N, i, j, rbs, cbs, bs = 2, nl = 4;
 11:   PetscMPIInt            rank, size;
 12:   PetscBool              partial = PETSC_FALSE, drop, flg;

 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 16:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 17:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 18:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 19:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, NULL));
 20:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-partial", &partial, NULL));

 22:   /* subdomains sharing a block of degrees of freedom, as in a finite element decomposition */
 23:   PetscCall(PetscMalloc1(nl, &idxs));
 24:   for (i = 0; i < nl; i++) idxs[i] = rank * (nl - 1) + i;
 25:   PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, bs, nl, idxs, PETSC_OWN_POINTER, &map));

 27:   N = (size * (nl - 1) + 1) * bs;
 28:   PetscCall(MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, N, N, map, map, &A));
 29:   PetscCall(MatISSetPreallocation(A, 3, NULL, 3, NULL));

 31:   /* leave the last local block without any entry, so that MatISFixLocalEmpty() drops it;
 32:      with -partial the last subdomain instead leaves a single degree of freedom out, which
 33:      breaks the block and must make every process give up the block size
 34:   */
 35:   for (i = 0; i < nl * bs; i++) {
 36:     drop = (PetscBool)(i / bs == nl - 1);
 37:     if (partial && rank == size - 1) drop = (PetscBool)(i == nl * bs - 1);
 38:     if (drop) continue;
 39:     v = 2.0;
 40:     PetscCall(MatSetValuesLocal(A, 1, &i, 1, &i, &v, ADD_VALUES));
 41:     for (j = i - 1; j <= i + 1; j += 2) {
 42:       if (j < 0 || j >= (partial && rank == size - 1 ? nl * bs - 1 : (nl - 1) * bs)) continue;
 43:       v = -1.0;
 44:       PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &v, ADD_VALUES));
 45:     }
 46:   }
 47:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 48:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 49:   PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));

 51:   PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
 52:   PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs));
 53:   PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs));
 54:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Block sizes before MatISFixLocalEmpty %" PetscInt_FMT " %" PetscInt_FMT "\n", rbs, cbs));

 56:   PetscCall(MatISFixLocalEmpty(A, PETSC_TRUE));
 57:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 58:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 60:   PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
 61:   PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs));
 62:   PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs));
 63:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Block sizes after MatISFixLocalEmpty %" PetscInt_FMT " %" PetscInt_FMT "\n", rbs, cbs));

 65:   PetscCall(MatMultEqual(A, B, 5, &flg));
 66:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatISFixLocalEmpty changed the operator");

 68:   PetscCall(MatDestroy(&A));
 69:   PetscCall(MatDestroy(&B));
 70:   PetscCall(ISLocalToGlobalMappingDestroy(&map));
 71:   PetscCall(PetscFinalize());
 72:   return 0;
 73: }

 75: /*TEST

 77:    test:
 78:       suffix: 1
 79:       nsize: 1

 81:    test:
 82:       suffix: 2
 83:       nsize: 3

 85:    test:
 86:       suffix: 1_partial
 87:       nsize: 1
 88:       args: -partial

 90:    test:
 91:       suffix: 2_partial
 92:       nsize: 3
 93:       args: -partial

 95: TEST*/