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