1: static const char help[] = "Tests PCASMGetSubKSP() ordering and reused submatrices across a change of operator.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat       A;
  8:   PC        pc;
  9:   IS        is;
 10:   Vec       x, b;
 11:   KSP      *subksp;
 12:   PetscInt  n            = 16, rstart, rend, nlocal;
 13:   PetscBool before_setup = PETSC_FALSE;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

 18:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 19:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 20:   PetscCall(MatSetFromOptions(A));
 21:   PetscCall(MatSetUp(A));
 22:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 23:   for (PetscInt row = rstart; row < rend; row++) {
 24:     PetscCall(MatSetValue(A, row, row, 2.0, INSERT_VALUES));
 25:     if (row > 0) PetscCall(MatSetValue(A, row, row - 1, -1.0, INSERT_VALUES));
 26:     if (row < n - 1) PetscCall(MatSetValue(A, row, row + 1, -1.0, INSERT_VALUES));
 27:   }
 28:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 29:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 31:   PetscCall(MatCreateVecs(A, &x, &b));
 32:   PetscCall(VecSet(b, 1.0));

 34:   PetscCall(PCCreate(PETSC_COMM_WORLD, &pc));
 35:   PetscCall(PCSetType(pc, PCASM));
 36:   PetscCall(PCSetOperators(pc, A, A));

 38:   PetscCall(ISCreateStride(PETSC_COMM_SELF, rend - rstart, rstart, 1, &is));
 39:   PetscCall(PCASMSetLocalSubdomains(pc, 1, &is, NULL));

 41:   /* Specifying the subdomains sets n_local_true, but the sub-KSPs are not
 42:      allocated until PCSetUp(). Querying them here must raise an error rather
 43:      than hand back a null array. */
 44:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-before_setup", &before_setup, NULL));
 45:   if (before_setup) PetscCall(PCASMGetSubKSP(pc, &nlocal, NULL, &subksp));

 47:   /* After setup the sub-KSPs exist. */
 48:   PetscCall(PCSetFromOptions(pc));
 49:   PetscCall(PCSetUp(pc));
 50:   PetscCall(PCASMGetSubKSP(pc, &nlocal, NULL, &subksp));
 51:   PetscCheck(nlocal == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected 1 local subdomain, got %" PetscInt_FMT, nlocal);
 52:   PetscCheck(subksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCASMGetSubKSP() returned a null array after PCSetUp()");
 53:   PetscCall(PCApply(pc, b, x));

 55:   /* A subsolver told to factor in place, e.g. -sub_pc_type ilu
 56:      -sub_pc_factor_in_place, has overwritten its submatrix with the factors
 57:      and left it flagged as factored. Changing the operator refills those
 58:      submatrices with MAT_REUSE_MATRIX, which a matrix still flagged as
 59:      factored rejects. */
 60:   PetscCall(MatScale(A, 2.0));
 61:   PetscCall(PCSetOperators(pc, A, A));
 62:   PetscCall(PCSetUp(pc));
 63:   PetscCall(PCApply(pc, b, x));

 65:   PetscCall(VecDestroy(&x));
 66:   PetscCall(VecDestroy(&b));
 67:   PetscCall(ISDestroy(&is));
 68:   PetscCall(PCDestroy(&pc));
 69:   PetscCall(MatDestroy(&A));
 70:   PetscCall(PetscFinalize());
 71:   return 0;
 72: }

 74: /*TEST

 76:    test:
 77:       suffix: 1
 78:       nsize: {{1 2}}
 79:       args: -sub_pc_type ilu -sub_pc_factor_in_place
 80:       output_file: output/empty.out

 82:    test:
 83:       suffix: 2
 84:       requires: !defined(PETSCTEST_VALGRIND) !defined(PETSC_HAVE_SANITIZER)
 85:       args: -before_setup -petsc_ci_portable_error_output -error_output_stdout
 86:       filter: grep -E "(PETSC ERROR)" | grep -E "(wrong order|Need to call|PCASMGetSubKSP\(\)|main\(\))"

 88: TEST*/