Actual source code: ex55.c

  1: static const char help[] = "Example demonstrating PCCOMPOSITE where one of the inner PCs uses a different operator\n\
  2: \n";

  4: #include <petscksp.h>

  6: int main(int argc, char **argv)
  7: {
  8:   PetscInt    n = 10, i, col[3];
  9:   Vec         x, b;
 10:   Mat         A, B;
 11:   KSP         ksp;
 12:   PC          pc, subpc;
 13:   PetscScalar value[3];

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 18:   /* Create a diagonal matrix with a given distribution of diagonal elements */
 19:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 20:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 21:   PetscCall(MatSetFromOptions(A));
 22:   PetscCall(MatSetUp(A));
 23:   /*
 24:      Assemble matrix
 25:   */
 26:   value[0] = -1.0;
 27:   value[1] = 2.0;
 28:   value[2] = -1.0;
 29:   for (i = 1; i < n - 1; i++) {
 30:     col[0] = i - 1;
 31:     col[1] = i;
 32:     col[2] = i + 1;
 33:     PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 34:   }
 35:   i      = n - 1;
 36:   col[0] = n - 2;
 37:   col[1] = n - 1;
 38:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 39:   i        = 0;
 40:   col[0]   = 0;
 41:   col[1]   = 1;
 42:   value[0] = 2.0;
 43:   value[1] = -1.0;
 44:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 45:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 46:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 48:   PetscCall(MatCreateVecs(A, &x, &b));

 50:   /* Create a KSP object */
 51:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 52:   PetscCall(KSPSetOperators(ksp, A, A));

 54:   /* Set up a composite preconditioner */
 55:   PetscCall(KSPGetPC(ksp, &pc));
 56:   PetscCall(PCSetType(pc, PCCOMPOSITE)); /* default composite with single Identity PC */
 57:   PetscCall(PCCompositeSetType(pc, PC_COMPOSITE_ADDITIVE));
 58:   PetscCall(PCCompositeAddPCType(pc, PCLU));
 59:   PetscCall(PCCompositeGetPC(pc, 0, &subpc));
 60:   /*  B is set to the diagonal of A; this demonstrates that setting the operator for a subpc changes the preconditioning */
 61:   PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
 62:   PetscCall(MatGetDiagonal(A, b));
 63:   PetscCall(MatDiagonalSet(B, b, ADD_VALUES));
 64:   PetscCall(PCSetOperators(subpc, B, B));
 65:   PetscCall(PCCompositeAddPCType(pc, PCNONE));

 67:   PetscCall(KSPSetFromOptions(ksp));
 68:   PetscCall(KSPSolve(ksp, b, x));

 70:   PetscCall(KSPDestroy(&ksp));
 71:   PetscCall(MatDestroy(&A));
 72:   PetscCall(MatDestroy(&B));
 73:   PetscCall(VecDestroy(&x));
 74:   PetscCall(VecDestroy(&b));
 75:   PetscCall(PetscFinalize());
 76:   return 0;
 77: }

 79: /*TEST

 81:    test:
 82:      args: -ksp_monitor -pc_composite_type multiplicative

 84: TEST*/