1 static char help[] = "Tests repeated setups and solves of PCFIELDSPLIT.\n\n";
2 #include <petscksp.h>
3
replace_submats(Mat A)4 static PetscErrorCode replace_submats(Mat A)
5 {
6 IS *r, *c;
7 PetscInt i, j, nr, nc;
8
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;
18
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 }
32
main(int argc,char * argv[])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
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));
59
60 PetscCall(MatDestroy(&M));
61
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));
74
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 }
90
91 /*TEST
92
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
97
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
103
104 TEST*/
105