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;
23: PetscInitialize(&argc, &argv, (char *)0, PETSC_NULL);
25: MatCreateAIJ(PETSC_COMM_WORLD, 5, 5, PETSC_DETERMINE, PETSC_DETERMINE, 3, NULL, 0, NULL, &subA[0]);
26: MatGetOwnershipRange(subA[0], &mstart, &mend);
27: for (row = mstart; row < mend; ++row) {
28: val = 1.0;
29: col = row;
30: MatSetValues(subA[0], 1, &row, 1, &col, &val, INSERT_VALUES);
31: }
32: MatAssemblyBegin(subA[0], MAT_FINAL_ASSEMBLY);
33: MatAssemblyEnd(subA[0], MAT_FINAL_ASSEMBLY);
35: MatCreateAIJ(PETSC_COMM_WORLD, 5, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[1]);
36: MatGetOwnershipRange(subA[1], &mstart, &mend);
37: for (row = mstart; row < mend; ++row) {
38: col = 1;
39: val = 0.0;
40: MatSetValues(subA[1], 1, &row, 1, &col, &val, INSERT_VALUES);
41: }
42: MatAssemblyBegin(subA[1], MAT_FINAL_ASSEMBLY);
43: MatAssemblyEnd(subA[1], MAT_FINAL_ASSEMBLY);
45: MatCreateAIJ(PETSC_COMM_WORLD, 5, 4, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[2]);
46: MatGetOwnershipRange(subA[2], &mstart, &mend);
47: for (row = mstart; row < mend; ++row) {
48: col = 1;
49: val = 0.0;
50: MatSetValues(subA[2], 1, &row, 1, &col, &val, INSERT_VALUES);
51: }
52: MatAssemblyBegin(subA[2], MAT_FINAL_ASSEMBLY);
53: MatAssemblyEnd(subA[2], MAT_FINAL_ASSEMBLY);
55: MatCreateAIJ(PETSC_COMM_WORLD, 3, 5, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[3]);
56: MatGetOwnershipRange(subA[3], &mstart, &mend);
57: for (row = mstart; row < mend; ++row) {
58: col = row;
59: val = 0.0;
60: MatSetValues(subA[3], 1, &row, 1, &col, &val, INSERT_VALUES);
61: }
62: MatAssemblyBegin(subA[3], MAT_FINAL_ASSEMBLY);
63: MatAssemblyEnd(subA[3], MAT_FINAL_ASSEMBLY);
65: MatCreateAIJ(PETSC_COMM_WORLD, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[4]);
66: MatGetOwnershipRange(subA[4], &mstart, &mend);
67: for (row = mstart; row < mend; ++row) {
68: col = row;
69: val = 4.0;
70: MatSetValues(subA[4], 1, &row, 1, &col, &val, INSERT_VALUES);
71: }
72: MatAssemblyBegin(subA[4], MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(subA[4], MAT_FINAL_ASSEMBLY);
75: MatCreateAIJ(PETSC_COMM_WORLD, 3, 4, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[5]);
76: MatGetOwnershipRange(subA[5], &mstart, &mend);
77: for (row = mstart; row < mend; ++row) {
78: col = row;
79: val = 0.0;
80: MatSetValues(subA[5], 1, &row, 1, &col, &val, INSERT_VALUES);
81: }
82: MatAssemblyBegin(subA[5], MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(subA[5], MAT_FINAL_ASSEMBLY);
85: MatCreateAIJ(PETSC_COMM_WORLD, 4, 5, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[6]);
86: MatGetOwnershipRange(subA[6], &mstart, &mend);
87: for (row = mstart; row < mend; ++row) {
88: col = 2;
89: val = 0.0;
90: MatSetValues(subA[6], 1, &row, 1, &col, &val, INSERT_VALUES);
91: }
92: MatAssemblyBegin(subA[6], MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(subA[6], MAT_FINAL_ASSEMBLY);
95: MatCreateAIJ(PETSC_COMM_WORLD, 4, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[7]);
96: MatGetOwnershipRange(subA[7], &mstart, &mend);
97: for (row = mstart; row < mend; ++row) {
98: col = 1;
99: val = 0.0;
100: MatSetValues(subA[7], 1, &row, 1, &col, &val, INSERT_VALUES);
101: }
102: MatAssemblyBegin(subA[7], MAT_FINAL_ASSEMBLY);
103: MatAssemblyEnd(subA[7], MAT_FINAL_ASSEMBLY);
105: MatCreateAIJ(PETSC_COMM_WORLD, 4, 4, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &subA[8]);
106: MatGetOwnershipRange(subA[8], &mstart, &mend);
107: for (row = mstart; row < mend; ++row) {
108: col = row;
109: val = 8.0;
110: MatSetValues(subA[8], 1, &row, 1, &col, &val, INSERT_VALUES);
111: }
112: MatAssemblyBegin(subA[8], MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(subA[8], MAT_FINAL_ASSEMBLY);
115: MatCreateNest(PETSC_COMM_WORLD, 3, NULL, 3, NULL, subA, &A);
116: MatNestGetISs(A, isg, NULL);
117: VecCreate(PETSC_COMM_WORLD, &subb[0]);
118: VecSetSizes(subb[0], 5, PETSC_DECIDE);
119: VecSetFromOptions(subb[0]);
120: VecSet(subb[0], 1.0);
122: VecCreate(PETSC_COMM_WORLD, &subb[1]);
123: VecSetSizes(subb[1], 3, PETSC_DECIDE);
124: VecSetFromOptions(subb[1]);
125: VecSet(subb[1], 2.0);
127: VecCreate(PETSC_COMM_WORLD, &subb[2]);
128: VecSetSizes(subb[2], 4, PETSC_DECIDE);
129: VecSetFromOptions(subb[2]);
130: VecSet(subb[2], 3.0);
132: VecCreateNest(PETSC_COMM_WORLD, 3, NULL, subb, &b);
133: VecDuplicate(b, &r);
134: VecCopy(b, r);
136: MatMult(A, b, r);
137: VecSet(b, 0.0);
139: KSPCreate(PETSC_COMM_WORLD, &ksp);
140: KSPSetOperators(ksp, A, A);
141: KSPSetFromOptions(ksp);
142: KSPGetPC(ksp, &pc);
143: PCFieldSplitSetIS(pc, "a", isg[0]);
144: PCFieldSplitSetIS(pc, "b", isg[1]);
145: PCFieldSplitSetIS(pc, "c", isg[2]);
147: KSPSolve(ksp, r, b);
148: KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);
150: MatDestroy(&subA[0]);
151: MatDestroy(&subA[1]);
152: MatDestroy(&subA[2]);
153: MatDestroy(&subA[3]);
154: MatDestroy(&subA[4]);
155: MatDestroy(&subA[5]);
156: MatDestroy(&subA[6]);
157: MatDestroy(&subA[7]);
158: MatDestroy(&subA[8]);
159: MatDestroy(&A);
160: VecDestroy(&subb[0]);
161: VecDestroy(&subb[1]);
162: VecDestroy(&subb[2]);
163: VecDestroy(&b);
164: VecDestroy(&r);
165: KSPDestroy(&ksp);
167: PetscFinalize();
168: return 0;
169: }
171: /*TEST
173: test:
174: args: -pc_type fieldsplit -ksp_monitor
176: TEST*/