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;
11: MatNestGetSubMats(A, &nr, &nc, NULL);
12: PetscMalloc1(nr, &r);
13: PetscMalloc1(nc, &c);
14: 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: MatCreateSubMatrix(A, r[i], c[j], MAT_INITIAL_MATRIX, &sA);
21: MatDuplicate(sA, MAT_COPY_VALUES, &nA);
22: MatGetOptionsPrefix(sA, &prefix);
23: MatSetOptionsPrefix(nA, prefix);
24: MatNestSetSubMat(A, i, j, nA);
25: MatDestroy(&nA);
26: MatDestroy(&sA);
27: }
28: }
29: PetscFree(r);
30: PetscFree(c);
31: return 0;
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;
45: PetscInitialize(&argc, &argv, (char *)0, help);
46: MatCreateAIJ(PETSC_COMM_WORLD, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &M);
47: MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
48: MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
49: MatShift(M, 1.);
50: MatGetOwnershipRange(M, &rstart, &rend);
51: ISCreateStride(PetscObjectComm((PetscObject)M), 7, rstart, 1, &f[0]);
52: ISComplement(f[0], rstart, rend, &f[1]);
53: for (i = 0; i < 2; i++) {
54: for (j = 0; j < 2; j++) {
55: MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sA[i][j]);
56: MatCreateSubMatrix(M, f[i], f[j], MAT_INITIAL_MATRIX, &sP[i][j]);
57: }
58: }
59: MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sA[0][0], &A);
60: MatCreateNest(PetscObjectComm((PetscObject)M), 2, f, 2, f, &sP[0][0], &P);
62: /* Tests MatMissingDiagonal_Nest */
63: MatMissingDiagonal(M, &missM, NULL);
64: MatMissingDiagonal(A, &missA, NULL);
65: if (missM != missA) PetscPrintf(PETSC_COMM_WORLD, "Unexpected %s != %s\n", missM ? "true" : "false", missA ? "true" : "false");
67: MatDestroy(&M);
69: KSPCreate(PetscObjectComm((PetscObject)A), &ksp);
70: KSPSetOperators(ksp, A, P);
71: KSPGetPC(ksp, &pc);
72: PCSetType(pc, PCFIELDSPLIT);
73: KSPSetFromOptions(ksp);
74: MatCreateVecs(A, &x, &b);
75: VecSetRandom(b, NULL);
76: KSPSolve(ksp, b, x);
77: replace_submats(A);
78: replace_submats(P);
79: KSPSolve(ksp, b, x);
81: KSPDestroy(&ksp);
82: VecDestroy(&x);
83: VecDestroy(&b);
84: MatDestroy(&A);
85: MatDestroy(&P);
86: for (i = 0; i < 2; i++) {
87: ISDestroy(&f[i]);
88: for (j = 0; j < 2; j++) {
89: MatDestroy(&sA[i][j]);
90: MatDestroy(&sP[i][j]);
91: }
92: }
93: 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*/