diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 840e119f7e303277a86d67bba562e8034fdf3f62..e39af9541e9bc697c0bb10703734670c3bd2dc56 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -15,8 +15,32 @@
 #include "linearSystemPETSc.hpp"
 
 template class linearSystemPETSc<double>;
+
+
 #ifdef PETSC_USE_COMPLEX
 template class linearSystemPETSc<std::complex<double> >;
+
+// this specialization will cast to a double (take the real part) if "val" is a double wheras Petsc is built in complex mode.
+template <>
+void linearSystemPETSc<double>::getFromRightHandSide(int row, double &val) const
+{
+  PetscScalar *tmp;
+  _try(VecGetArray(_b, &tmp));
+  PetscScalar s = tmp[row];
+  _try(VecRestoreArray(_b, &tmp));
+  val = s.real();
+}
+
+// this specialization will cast to a double (take the real part) if "val" is a double wheras Petsc is built in complex mode.
+template <>
+void linearSystemPETSc<double>::getFromSolution(int row, double &val) const
+{
+  PetscScalar *tmp;
+  _try(VecGetArray(_x, &tmp));
+  PetscScalar s = tmp[row];
+  _try(VecRestoreArray(_x, &tmp));
+  val = s.real();
+}
 #endif
 
 template<>
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index d835b73ecbd5abbc571a08a1f73558c5159687c3..63f0c524baaee40cb02da5f387b22f74339ec338 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/Solver/linearSystemPETSc.hpp
@@ -266,20 +266,16 @@ void linearSystemPETSc<scalar>::addToRightHandSide(int row, const scalar &val)
   PetscScalar s = val;
   _try(VecSetValues(_b, 1, &i, &s, ADD_VALUES));
 }
+#if defined(PETSC_USE_COMPLEX)
+// this specialization will cast to a double (take the real part) if "val" is a double wheras Petsc is built in complex mode.
+template <>
+void linearSystemPETSc<double>::getFromRightHandSide(int row, double &val) const;
+#endif
 
 template <class scalar>
 void linearSystemPETSc<scalar>::getFromRightHandSide(int row, scalar &val) const
 {
-#if defined(PETSC_USE_COMPLEX)
-  PetscScalar *tmp;
-  _try(VecGetArray(_b, &tmp));
-  PetscScalar s = tmp[row];
-  _try(VecRestoreArray(_b, &tmp));
-  // FIXME specialize this routine
-  val = s.real();
-#else
   _try(VecGetValues(_b, 1, &row, &val));
-#endif
 }
 
 template <class scalar>
@@ -310,19 +306,16 @@ void linearSystemPETSc<scalar>::addToSolution(int row, const scalar &val)
   PetscScalar s = val;
   _try(VecSetValues(_x, 1, &i, &s, ADD_VALUES));
 }
+#if defined(PETSC_USE_COMPLEX)
+// this specialization will cast to a double (take the real part) if "val" is a double wheras Petsc is built in complex mode.
+template <>
+void linearSystemPETSc<double>::getFromSolution(int row, double &val) const;
+#endif
 
 template <class scalar>
 void linearSystemPETSc<scalar>::getFromSolution(int row, scalar &val) const
 {
-#if defined(PETSC_USE_COMPLEX)
-  PetscScalar *tmp;
-  _try(VecGetArray(_x, &tmp));
-  PetscScalar s = tmp[row];
-  _try(VecRestoreArray(_x, &tmp));
-  val = s.real();
-#else
   _try(VecGetValues(_x, 1, &row, &val));
-#endif
 }
 
 template <class scalar>