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*/