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