Actual source code: ex57.c
1: /*
2: Tests PCFIELDSPLIT and hence VecGetRestoreArray_Nest() usage in VecScatter
4: Example contributed by: Mike Wick <michael.wick.1980@gmail.com>
5: */
6: #include <petscksp.h>
8: int main(int argc, char **argv)
9: {
10: Mat A;
11: Mat subA[9];
12: IS isg[3];
13: PetscInt row, col, mstart, mend;
14: PetscScalar val;
15: Vec subb[3];
16: Vec b;
17: Vec r;
18: KSP ksp;
19: PC pc;
21: PetscFunctionBeginUser;
22: PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
24: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 5, 5, PETSC_DETERMINE, PETSC_DETERMINE, &subA[0]));
25: PetscCall(MatGetOwnershipRange(subA[0], &mstart, &mend));
26: for (row = mstart; row < mend; ++row) {
27: val = 1.0;
28: col = row;
29: PetscCall(MatSetValues(subA[0], 1, &row, 1, &col, &val, INSERT_VALUES));
30: }
31: PetscCall(MatAssemblyBegin(subA[0], MAT_FINAL_ASSEMBLY));
32: PetscCall(MatAssemblyEnd(subA[0], MAT_FINAL_ASSEMBLY));
34: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 5, 3, PETSC_DETERMINE, PETSC_DETERMINE, &subA[1]));
35: PetscCall(MatGetOwnershipRange(subA[1], &mstart, &mend));
36: for (row = mstart; row < mend; ++row) {
37: col = 1;
38: val = 0.0;
39: PetscCall(MatSetValues(subA[1], 1, &row, 1, &col, &val, INSERT_VALUES));
40: }
41: PetscCall(MatAssemblyBegin(subA[1], MAT_FINAL_ASSEMBLY));
42: PetscCall(MatAssemblyEnd(subA[1], MAT_FINAL_ASSEMBLY));
44: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 5, 4, PETSC_DETERMINE, PETSC_DETERMINE, &subA[2]));
45: PetscCall(MatGetOwnershipRange(subA[2], &mstart, &mend));
46: for (row = mstart; row < mend; ++row) {
47: col = 1;
48: val = 0.0;
49: PetscCall(MatSetValues(subA[2], 1, &row, 1, &col, &val, INSERT_VALUES));
50: }
51: PetscCall(MatAssemblyBegin(subA[2], MAT_FINAL_ASSEMBLY));
52: PetscCall(MatAssemblyEnd(subA[2], MAT_FINAL_ASSEMBLY));
54: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 3, 5, PETSC_DETERMINE, PETSC_DETERMINE, &subA[3]));
55: PetscCall(MatGetOwnershipRange(subA[3], &mstart, &mend));
56: for (row = mstart; row < mend; ++row) {
57: col = row;
58: val = 0.0;
59: PetscCall(MatSetValues(subA[3], 1, &row, 1, &col, &val, INSERT_VALUES));
60: }
61: PetscCall(MatAssemblyBegin(subA[3], MAT_FINAL_ASSEMBLY));
62: PetscCall(MatAssemblyEnd(subA[3], MAT_FINAL_ASSEMBLY));
64: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, &subA[4]));
65: PetscCall(MatGetOwnershipRange(subA[4], &mstart, &mend));
66: for (row = mstart; row < mend; ++row) {
67: col = row;
68: val = 4.0;
69: PetscCall(MatSetValues(subA[4], 1, &row, 1, &col, &val, INSERT_VALUES));
70: }
71: PetscCall(MatAssemblyBegin(subA[4], MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(subA[4], MAT_FINAL_ASSEMBLY));
74: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 3, 4, PETSC_DETERMINE, PETSC_DETERMINE, &subA[5]));
75: PetscCall(MatGetOwnershipRange(subA[5], &mstart, &mend));
76: for (row = mstart; row < mend; ++row) {
77: col = row;
78: val = 0.0;
79: PetscCall(MatSetValues(subA[5], 1, &row, 1, &col, &val, INSERT_VALUES));
80: }
81: PetscCall(MatAssemblyBegin(subA[5], MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(subA[5], MAT_FINAL_ASSEMBLY));
84: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 4, 5, PETSC_DETERMINE, PETSC_DETERMINE, &subA[6]));
85: PetscCall(MatGetOwnershipRange(subA[6], &mstart, &mend));
86: for (row = mstart; row < mend; ++row) {
87: col = 2;
88: val = 0.0;
89: PetscCall(MatSetValues(subA[6], 1, &row, 1, &col, &val, INSERT_VALUES));
90: }
91: PetscCall(MatAssemblyBegin(subA[6], MAT_FINAL_ASSEMBLY));
92: PetscCall(MatAssemblyEnd(subA[6], MAT_FINAL_ASSEMBLY));
94: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 4, 3, PETSC_DETERMINE, PETSC_DETERMINE, &subA[7]));
95: PetscCall(MatGetOwnershipRange(subA[7], &mstart, &mend));
96: for (row = mstart; row < mend; ++row) {
97: col = 1;
98: val = 0.0;
99: PetscCall(MatSetValues(subA[7], 1, &row, 1, &col, &val, INSERT_VALUES));
100: }
101: PetscCall(MatAssemblyBegin(subA[7], MAT_FINAL_ASSEMBLY));
102: PetscCall(MatAssemblyEnd(subA[7], MAT_FINAL_ASSEMBLY));
104: PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 4, 4, PETSC_DETERMINE, PETSC_DETERMINE, &subA[8]));
105: PetscCall(MatGetOwnershipRange(subA[8], &mstart, &mend));
106: for (row = mstart; row < mend; ++row) {
107: col = row;
108: val = 8.0;
109: PetscCall(MatSetValues(subA[8], 1, &row, 1, &col, &val, INSERT_VALUES));
110: }
111: PetscCall(MatAssemblyBegin(subA[8], MAT_FINAL_ASSEMBLY));
112: PetscCall(MatAssemblyEnd(subA[8], MAT_FINAL_ASSEMBLY));
114: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 3, NULL, 3, NULL, subA, &A));
115: PetscCall(MatNestGetISs(A, isg, NULL));
116: PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 5, PETSC_DECIDE, &subb[0]));
117: PetscCall(VecSet(subb[0], 1.0));
119: PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 3, PETSC_DECIDE, &subb[1]));
120: PetscCall(VecSet(subb[1], 2.0));
122: PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 4, PETSC_DECIDE, &subb[2]));
123: PetscCall(VecSet(subb[2], 3.0));
125: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 3, NULL, subb, &b));
126: PetscCall(VecDuplicate(b, &r));
127: PetscCall(VecCopy(b, r));
129: PetscCall(MatMult(A, b, r));
130: PetscCall(VecSet(b, 0.0));
132: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
133: PetscCall(KSPSetOperators(ksp, A, A));
134: PetscCall(KSPSetFromOptions(ksp));
135: PetscCall(KSPGetPC(ksp, &pc));
136: PetscCall(PCFieldSplitSetIS(pc, "a", isg[0]));
137: PetscCall(PCFieldSplitSetIS(pc, "b", isg[1]));
138: PetscCall(PCFieldSplitSetIS(pc, "c", isg[2]));
140: PetscCall(KSPSolve(ksp, r, b));
141: PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));
143: PetscCall(MatDestroy(&subA[0]));
144: PetscCall(MatDestroy(&subA[1]));
145: PetscCall(MatDestroy(&subA[2]));
146: PetscCall(MatDestroy(&subA[3]));
147: PetscCall(MatDestroy(&subA[4]));
148: PetscCall(MatDestroy(&subA[5]));
149: PetscCall(MatDestroy(&subA[6]));
150: PetscCall(MatDestroy(&subA[7]));
151: PetscCall(MatDestroy(&subA[8]));
152: PetscCall(MatDestroy(&A));
153: PetscCall(VecDestroy(&subb[0]));
154: PetscCall(VecDestroy(&subb[1]));
155: PetscCall(VecDestroy(&subb[2]));
156: PetscCall(VecDestroy(&b));
157: PetscCall(VecDestroy(&r));
158: PetscCall(KSPDestroy(&ksp));
160: PetscCall(PetscFinalize());
161: return 0;
162: }
164: /*TEST
166: test:
167: args: -pc_type fieldsplit -ksp_monitor
169: TEST*/