Actual source code: ex9.c
2: static char help[] = "Tests repeated setups and solves of PCFIELDSPLIT.\n\n";
3: #include <petscksp.h>
5: static PetscErrorCode replace_submats(Mat A)
6: {
7: IS *r, *c;
8: PetscInt i, j, nr, nc;
10: PetscFunctionBeginUser;
11: PetscCall(MatNestGetSubMats(A, &nr, &nc, NULL));
12: PetscCall(PetscMalloc1(nr, &r));
13: PetscCall(PetscMalloc1(nc, &c));
14: PetscCall(MatNestGetISs(A, r, c));
15: for (i = 0; i < nr; i++) {
16: for (j = 0; j < nc; j++) {
17: Mat sA, nA;
18: const char *prefix;
20: PetscCall(MatCreateSubMatrix(A, r[i], c[j], MAT_INITIAL_MATRIX, &sA));
21: PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &nA));
22: PetscCall(MatGetOptionsPrefix(sA, &prefix));
23: PetscCall(MatSetOptionsPrefix(nA, prefix));
24: PetscCall(MatNestSetSubMat(A, i, j, nA));
25: PetscCall(MatDestroy(&nA));
26: PetscCall(MatDestroy(&sA));
27: }
28: }
29: PetscCall(PetscFree(r));
30: PetscCall(PetscFree(c));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: int main(int argc, char *argv[])
35: {
36: KSP ksp;
37: PC pc;
38: Mat M, A, P, sA[2][2], sP[2][2];
39: Vec x, b;
40: IS f[2];
41: PetscInt i, j, rstart, rend;
42: PetscBool missA, missM;
44: PetscFunctionBeginUser;
45: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
46: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &M));
47: PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
48: PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
49: PetscCall(MatShift(M, 1.));
50: PetscCall(MatGetOwnershipRange(M, &rstart, &rend));
51: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)M), 7, rstart, 1, &f[0]));
52: PetscCall(ISComplement(f[0], rstart, rend, &f[1]));
53: for (i = 0; i < 2; i++) {
54: for (j = 0; j < 2; j++) {
55: PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sA[i][j]));
56: PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sP[i][j]));
57: }
58: }
59: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sA[0][0], &A));
60: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sP[0][0], &P));
62: /* Tests MatMissingDiagonal_Nest */
63: PetscCall(MatMissingDiagonal(M, &missM, NULL));
64: PetscCall(MatMissingDiagonal(A, &missA, NULL));
65: if (missM != missA) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Unexpected %s != %s\n", missM ? "true" : "false", missA ? "true" : "false"));
67: PetscCall(MatDestroy(&M));
69: PetscCall(KSPCreate(PetscObjectComm((PetscObject)A), &ksp));
70: PetscCall(KSPSetOperators(ksp, A, P));
71: PetscCall(KSPGetPC(ksp, &pc));
72: PetscCall(PCSetType(pc, PCFIELDSPLIT));
73: PetscCall(KSPSetFromOptions(ksp));
74: PetscCall(MatCreateVecs(A, &x, &b));
75: PetscCall(VecSetRandom(b, NULL));
76: PetscCall(KSPSolve(ksp, b, x));
77: PetscCall(replace_submats(A));
78: PetscCall(replace_submats(P));
79: PetscCall(KSPSolve(ksp, b, x));
81: PetscCall(KSPDestroy(&ksp));
82: PetscCall(VecDestroy(&x));
83: PetscCall(VecDestroy(&b));
84: PetscCall(MatDestroy(&A));
85: PetscCall(MatDestroy(&P));
86: for (i = 0; i < 2; i++) {
87: PetscCall(ISDestroy(&f[i]));
88: for (j = 0; j < 2; j++) {
89: PetscCall(MatDestroy(&sA[i][j]));
90: PetscCall(MatDestroy(&sP[i][j]));
91: }
92: }
93: PetscCall(PetscFinalize());
94: return 0;
95: }
97: /*TEST
99: test:
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 {{additive multiplicative}} -ksp_converged_reason -ksp_error_if_not_converged
104: test:
105: suffix: schur
106: nsize: 1
107: filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
108: 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
110: TEST*/