Actual source code: ex9.c
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;
41: PetscBool missA, missM;
43: PetscFunctionBeginUser;
44: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
45: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 10, 10, PETSC_DECIDE, PETSC_DECIDE, &M));
46: PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
47: PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
48: PetscCall(MatShift(M, 1.));
49: PetscCall(MatGetOwnershipRange(M, &rstart, &rend));
50: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)M), 7, rstart, 1, &f[0]));
51: PetscCall(ISComplement(f[0], rstart, rend, &f[1]));
52: for (i = 0; i < 2; i++) {
53: for (j = 0; j < 2; j++) {
54: PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sA[i][j]));
55: PetscCall(MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sP[i][j]));
56: }
57: }
58: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sA[0][0], &A));
59: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sP[0][0], &P));
61: /* Tests MatMissingDiagonal_Nest */
62: PetscCall(MatMissingDiagonal(M, &missM, NULL));
63: PetscCall(MatMissingDiagonal(A, &missA, NULL));
64: if (missM != missA) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Unexpected %s != %s\n", missM ? "true" : "false", missA ? "true" : "false"));
66: PetscCall(MatDestroy(&M));
68: PetscCall(KSPCreate(PetscObjectComm((PetscObject)A), &ksp));
69: PetscCall(KSPSetOperators(ksp, A, P));
70: PetscCall(KSPGetPC(ksp, &pc));
71: PetscCall(PCSetType(pc, PCFIELDSPLIT));
72: PetscCall(KSPSetFromOptions(ksp));
73: PetscCall(MatCreateVecs(A, &x, &b));
74: PetscCall(VecSetRandom(b, NULL));
75: PetscCall(KSPSolve(ksp, b, x));
76: PetscCall(replace_submats(A));
77: PetscCall(replace_submats(P));
78: PetscCall(KSPSolve(ksp, b, x));
79: PetscCall(KSPSolveTranspose(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*/