1: #include <petscksp.h>

  3: /* ------------------------------------------------------- */

  5: PetscErrorCode RunTest(void)
  6: {
  7:   PetscInt  N = 100, its = 0;
  8:   PetscBool draw = PETSC_FALSE, test = PETSC_FALSE;
  9:   PetscReal rnorm;
 10:   Mat       A;
 11:   Vec       b, x, r;
 12:   KSP       ksp;
 13:   PC        pc;

 15:   PetscFunctionBegin;

 17:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
 18:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test", &test, NULL));
 19:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw", &draw, NULL));

 21:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 22:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N));
 23:   PetscCall(MatSetType(A, MATPYTHON));
 24:   PetscCall(MatPythonSetType(A, "example100.py:Laplace1D"));
 25:   PetscCall(MatSetUp(A));

 27:   PetscCall(MatCreateVecs(A, &x, &b));
 28:   PetscCall(VecSet(b, 1));

 30:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 31:   PetscCall(KSPSetType(ksp, KSPPYTHON));
 32:   PetscCall(KSPPythonSetType(ksp, "example100.py:ConjGrad"));

 34:   PetscCall(KSPGetPC(ksp, &pc));
 35:   PetscCall(PCSetType(pc, PCPYTHON));
 36:   PetscCall(PCPythonSetType(pc, "example100.py:Jacobi"));

 38:   PetscCall(KSPSetOperators(ksp, A, A));
 39:   PetscCall(KSPSetFromOptions(ksp));
 40:   PetscCall(KSPSolve(ksp, b, x));

 42:   if (test) {
 43:     PetscCall(KSPGetTotalIterations(ksp, &its));
 44:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of KSP iterations = %" PetscInt_FMT "\n", its));
 45:   } else {
 46:     PetscCall(VecDuplicate(b, &r));
 47:     PetscCall(MatMult(A, x, r));
 48:     PetscCall(VecAYPX(r, -1, b));
 49:     PetscCall(VecNorm(r, NORM_2, &rnorm));
 50:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "error norm = %g\n", (double)rnorm));
 51:     PetscCall(VecDestroy(&r));
 52:   }

 54:   if (draw) {
 55:     PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
 56:     PetscCall(PetscSleep(2));
 57:   }

 59:   PetscCall(VecDestroy(&x));
 60:   PetscCall(VecDestroy(&b));
 61:   PetscCall(MatDestroy(&A));
 62:   PetscCall(KSPDestroy(&ksp));

 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: /* ------------------------------------------------------- */

 69: static char help[] = "Python-implemented Mat/KSP/PC.\n\n";

 71: #if !defined(PYTHON_EXE)
 72:   #define PYTHON_EXE 0
 73: #endif
 74: #if !defined(PYTHON_LIB)
 75:   #define PYTHON_LIB 0
 76: #endif

 78: int main(int argc, char *argv[])
 79: {
 80:   PetscFunctionBeginUser;
 81:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
 82:   PetscCall(PetscPythonInitialize(PYTHON_EXE, PYTHON_LIB));
 83:   PetscCall(RunTest());
 84:   PetscCall(PetscPythonPrintError());
 85:   PetscCall(PetscFinalize());
 86:   return 0;
 87: }

 89: /*TEST

 91:     test:
 92:       args: -ksp_monitor_short
 93:       requires: petsc4py
 94:       localrunfiles: example100.py

 96: TEST*/