diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
index d4e0e3300..f7d3b2281 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
@@ -107,9 +107,10 @@ namespace Opm
 
 
     /// Construct a system solver.
-    /// \param[in] linsolver   linear solver to use
-    NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& /*param*/)
+    NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
     {
+        use_amg_ = param.getDefault("cpr_use_amg", false);
+        use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true);
     }
 
 
@@ -185,7 +186,7 @@ namespace Opm
         // typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
         typedef Opm::CPRPreconditioner<Mat,Vector,Vector> Preconditioner;
         const double relax = 1.0;
-        Preconditioner precond(istlA, istlAe, relax);
+        Preconditioner precond(istlA, istlAe, relax, use_amg_, use_bicgstab_);
 
         // Construct linear solver.
         const double tolerance = 1e-3;
diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
index c8101ccef..fb2bd2b73 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
@@ -42,7 +42,9 @@ namespace Opm
         /// \param[in] param   parameters controlling the behaviour of
         ///                    the preconditioning and choice of
         ///                    linear solvers.
-        ///                    Note: parameters currently unused.
+        ///                    Parameters:
+        ///                        cpr_use_amg      (default false) if true, use AMG preconditioner for elliptic part
+        ///                        cpr_use_bicgstab (default true)  if true, use BiCGStab (else use CG) for elliptic part
         NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param);
 
         /// Solve the system of linear equations Ax = b, with A being the
@@ -51,6 +53,10 @@ namespace Opm
         /// \param[in] residual   residual object containing A and b.
         /// \return               the solution x
         virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
+
+    private:
+        bool use_amg_;
+        bool use_bicgstab_;
     };
 
 } // namespace Opm