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;
 16:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
 17:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test", &test, NULL));
 18:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw", &draw, NULL));

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

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

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

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

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

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

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

 58:   PetscCall(VecDestroy(&x));
 59:   PetscCall(VecDestroy(&b));
 60:   PetscCall(MatDestroy(&A));
 61:   PetscCall(KSPDestroy(&ksp));
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

 65: /* ------------------------------------------------------- */

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

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

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

 87: /*TEST

 89:     test:
 90:       args: -ksp_monitor_short
 91:       requires: petsc4py
 92:       localrunfiles: example100.py

 94: TEST*/