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, (char *)0, 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*/