Actual source code: ex23.c

  1: const char help[] = "Test whether KSPSolve() can be executed without zeroing output for KSPPREONLY";

  3: #include <petscksp.h>

  5: static PetscErrorCode PCApply_scale(PC pc, Vec x, Vec y)
  6: {
  7:   PetscFunctionBegin;
  8:   PetscCall(VecCopy(x, y));
  9:   PetscCall(VecScale(y, 0.5));
 10:   PetscFunctionReturn(PETSC_SUCCESS);
 11: }

 13: static PetscErrorCode VecSet_Error(Vec x, PetscReal _alpha)
 14: {
 15:   PetscFunctionBegin;
 16:   SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot zero this vector");
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: int main(int argc, char **argv)
 21: {
 22:   Mat      mat;
 23:   PetscInt M = 10;
 24:   MPI_Comm comm;
 25:   Vec      x, y;
 26:   KSP      ksp;
 27:   PC       pc;

 29:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 30:   comm = PETSC_COMM_WORLD;
 31:   PetscCall(MatCreateConstantDiagonal(comm, PETSC_DETERMINE, PETSC_DETERMINE, M, M, 2.0, &mat));
 32:   PetscCall(MatCreateVecs(mat, &x, &y));
 33:   PetscCall(VecSet(x, 1.0));
 34:   PetscCall(VecSet(y, 3.0));
 35:   PetscCall(VecSetOperation(y, VECOP_SET, (void (*)(void))VecSet_Error));

 37:   PetscCall(KSPCreate(comm, &ksp));
 38:   PetscCall(KSPSetOperators(ksp, mat, mat));
 39:   PetscCall(KSPSetType(ksp, KSPPREONLY));
 40:   PetscCall(KSPGetPC(ksp, &pc));
 41:   PetscCall(PCSetType(pc, PCSHELL));

 43:   PetscCall(PCShellSetApply(pc, PCApply_scale));
 44:   PetscCall(KSPSolve(ksp, x, y));

 46:   PetscCall(KSPDestroy(&ksp));
 47:   PetscCall(VecDestroy(&y));
 48:   PetscCall(VecDestroy(&x));
 49:   PetscCall(MatDestroy(&mat));
 50:   PetscCall(PetscFinalize());
 51:   return 0;
 52: }

 54: /*TEST

 56:   test:
 57:     suffix: 0

 59: TEST*/