diff --git a/CMakeLists.txt b/CMakeLists.txt
index 871ec54b56e7d89a9f1bffdc385955f13b43faf1..1a68e818d99b6ebed038779290a866cd02e14384 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -63,6 +63,7 @@ opt(MUMPS "Enable MUMPS sparse direct linear solver" OFF)
 opt(NATIVE_FILE_CHOOSER "Enable native file chooser in GUI" ${DEFAULT})
 opt(NETGEN "Enable Netgen 3D frontal mesh generator" ${DEFAULT})
 opt(NUMPY "Enable conversion between fullMatrix and numpy array (requires SWIG)" OFF)
+opt(PETSC4PY "Enable petsc4py wrappers for petsc matrices" ON)
 opt(OCC "Enable Open CASCADE geometrical models" ${DEFAULT})
 opt(ONELAB "Enable ONELAB solver interface" ${DEFAULT})
 opt(ONELAB2 "Enable experimental ONELAB-Cloud solver interface" OFF)
@@ -1120,6 +1121,18 @@ if(HAVE_PYTHON)
       set_config_option(HAVE_NUMPY "Numpy")
     endif(NUMPY_INC)
   endif(ENABLE_NUMPY)
+  if(HAVE_PETSC)
+    if(ENABLE_PETSC4PY)
+      EXEC_PROGRAM (${PYTHON_EXECUTABLE}
+        ARGS "-c \"import petsc4py; print(petsc4py.get_include())\""
+        OUTPUT_VARIABLE PETSC4PY_INC
+        RETURN_VALUE PETSC4PY_NOT_FOUND)
+      if(PETSC4PY_INC)
+        list(APPEND EXTERNAL_INCLUDES ${PETSC4PY_INC})
+        set_config_option(HAVE_PETSC4PY "PETSc4py")
+      endif(PETSC4PY_INC)
+    endif(ENABLE_PETSC4PY)
+  endif(HAVE_PETSC)
 endif(HAVE_PYTHON)
 
 check_function_exists(vsnprintf HAVE_VSNPRINTF)
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index b6e2c4b18379bd150b094caa85bd0d16565e85da..12baab72de7d6d108633f955e0e3f84079af93c5 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -54,6 +54,7 @@
 #cmakedefine HAVE_OSMESA
 #cmakedefine HAVE_PARSER
 #cmakedefine HAVE_PETSC
+#cmakedefine HAVE_PETSC4PY
 #cmakedefine HAVE_PLUGINS
 #cmakedefine HAVE_POST
 #cmakedefine HAVE_POPPLER
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index a1bc1e4f739ae4ac802df9030891dcd5eb0cc97b..105ae4ad35791fff311518855a17f6a1d071571d 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -114,7 +114,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   void zeroSolution();
   void printMatlab(const char *filename) const;
   virtual int systemSolve();
-  Mat &getMatrix(){ return _a; }
+  Mat getMatrix(){ return _a; }
   //std::vector<scalar> getData();
   //std::vector<int> getRowPointers();
   //std::vector<int> getColumnsIndices();
diff --git a/wrappers/gmshpy/gmshSolver.i b/wrappers/gmshpy/gmshSolver.i
index 6c22c518dd5088610ba635137cdb65074af16b63..2aec16dffa69c2df750daedb4b57bca0167ef3d4 100644
--- a/wrappers/gmshpy/gmshSolver.i
+++ b/wrappers/gmshpy/gmshSolver.i
@@ -36,6 +36,9 @@
 %include "linearSystemFull.h"
 %template(linearSystemFullDouble) linearSystemFull<double> ;
 #if defined(HAVE_PETSC)
+#if defined(HAVE_PETSC4PY)
+%include petsc4py/petsc4py.i
+#endif
 %include "linearSystemPETSc.h"
 %template(linearSystemPETScDouble) linearSystemPETSc<double>;
 %template(linearSystemPETScBlockDouble) linearSystemPETSc<fullMatrix<double> >;