1: static char help[] = "Tests repeated setups and solves of PCFIELDSPLIT.\n\n";
  2: #include <petscksp.h>

  4: static PetscErrorCode replace_submats(Mat A)
  5: {
  6:   IS      *r, *c;
  7:   PetscInt i, j, nr, nc;

  9:   PetscFunctionBeginUser;
 10:   PetscCall(MatNestGetSubMats(A, &nr, &nc, NULL));
 11:   PetscCall(PetscMalloc1(nr, &r));
 12:   PetscCall(PetscMalloc1(nc, &c));
 13:   PetscCall(MatNestGetISs(A, r, c));
 14:   for (i = 0; i < nr; i++) {
 15:     for (j = 0; j < nc; j++) {
 16:       Mat         sA, nA;
 17:       const char *prefix;

 19:       PetscCall(MatCreateSubMatrix(A, r[i], c[j], MAT_INITIAL_MATRIX, &sA));
 20:       PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &nA));
 21:       PetscCall(MatGetOptionsPrefix(sA, &prefix));
 22:       PetscCall(MatSetOptionsPrefix(nA, prefix));
 23:       PetscCall(MatNestSetSubMat(A, i, j, nA));
 24:       PetscCall(MatDestroy(&nA));
 25:       PetscCall(MatDestroy(&sA));
 26:     }
 27:   }
 28:   PetscCall(PetscFree(r));
 29:   PetscCall(PetscFree(c));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: int main(int argc, char *argv[])
 34: {
 35:   KSP      ksp;
 36:   PC       pc;
 37:   Mat      M, A, P, sA[2][2], sP[2][2];
 38:   Vec      x, b;
 39:   IS       f[2];
 40:   PetscInt i, j, rstart, rend;

 42:   PetscFunctionBeginUser;
 43:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 44:   PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 10, 10, PETSC_DECIDE, PETSC_DECIDE, &M));
 45:   PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
 46:   PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
 47:   PetscCall(MatShift(M, 1.));
 48:   PetscCall(MatGetOwnershipRange(M, &rstart, &rend));
 49:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)M), 7, rstart, 1, &f[0]));
 50:   PetscCall(ISComplement(f[0], rstart, rend, &f[1]));
 51:   for (i = 0; i < 2; i++) {
 52:     for (j = 0; j < 2; j++) {
 53:       PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sA[i][j]));
 54:       PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sP[i][j]));
 55:     }
 56:   }
 57:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sA[0][0], &A));
 58:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sP[0][0], &P));

 60:   PetscCall(MatDestroy(&M));

 62:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)A), &ksp));
 63:   PetscCall(KSPSetOperators(ksp, A, P));
 64:   PetscCall(KSPGetPC(ksp, &pc));
 65:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
 66:   PetscCall(KSPSetFromOptions(ksp));
 67:   PetscCall(MatCreateVecs(A, &x, &b));
 68:   PetscCall(VecSetRandom(b, NULL));
 69:   PetscCall(KSPSolve(ksp, b, x));
 70:   PetscCall(replace_submats(A));
 71:   PetscCall(replace_submats(P));
 72:   PetscCall(KSPSolve(ksp, b, x));
 73:   PetscCall(KSPSolveTranspose(ksp, b, x));

 75:   PetscCall(KSPDestroy(&ksp));
 76:   PetscCall(VecDestroy(&x));
 77:   PetscCall(VecDestroy(&b));
 78:   PetscCall(MatDestroy(&A));
 79:   PetscCall(MatDestroy(&P));
 80:   for (i = 0; i < 2; i++) {
 81:     PetscCall(ISDestroy(&f[i]));
 82:     for (j = 0; j < 2; j++) {
 83:       PetscCall(MatDestroy(&sA[i][j]));
 84:       PetscCall(MatDestroy(&sP[i][j]));
 85:     }
 86:   }
 87:   PetscCall(PetscFinalize());
 88:   return 0;
 89: }

 91: /*TEST

 93:    test:
 94:      nsize: 1
 95:      filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
 96:      args: -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_type {{additive multiplicative}} -ksp_converged_reason -ksp_error_if_not_converged

 98:    test:
 99:      suffix: schur
100:      nsize: 1
101:      filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
102:      args: -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_diag_use_amat {{0 1}} -pc_fieldsplit_type schur -pc_fieldsplit_schur_scale 1.0 -pc_fieldsplit_schur_fact_type {{diag lower upper full}} -ksp_converged_reason -ksp_error_if_not_converged

104: TEST*/