diff --git a/dune/solvers/solvers/linearipopt.hh b/dune/solvers/solvers/linearipopt.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2d9689b73ad14601d6a5d97e459593569c8e3589
--- /dev/null
+++ b/dune/solvers/solvers/linearipopt.hh
@@ -0,0 +1,604 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_LINEAR_IPOPT_SOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_LINEAR_IPOPT_SOLVER_HH
+
+/** \file
+ \brief A wrapper for the IPOpt solver for linear optimization problem
+ with box constraints and linear constraints
+*/
+#include <algorithm>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+
+#include "IpTNLP.hpp"
+#include "IpIpoptApplication.hpp"
+
+/** \brief Implementation class used to pipe linear problems to
+ the IPOpt interior-point solver
+
+ \tparam JacobianType The type of the jacobian of the linear constraints
+*/
+template <class VectorType, class JacobianType>
+class LinearIPOptProblem : public Ipopt::TNLP
+{
+ enum {blocksize = VectorType::block_type::dimension};
+ enum {constraintBlocksize = JacobianType::block_type::rows};
+
+ typedef typename VectorType::field_type field_type;
+
+public:
+ /** Setup problem if no linear constraints are present */
+ LinearIPOptProblem( VectorType* x, const VectorType* rhs,
+ const Dune::BitSetVector<blocksize>* dirichletNodes,
+ std::vector<BoxConstraint<field_type, blocksize> >* obstacles)
+ : x_(x), rhs_(rhs),
+ dirichletNodes_(dirichletNodes), obstacles_(obstacles),
+ constraintMatrix_(nullptr), constraintObstacles_(nullptr), jacobianCutoff_(1e-8)
+ {}
+
+ /** Setup problem full problem */
+ LinearIPOptProblem(VectorType* x, const VectorType* rhs,
+ const Dune::BitSetVector<blocksize>* dirichletNodes,
+ std::vector<BoxConstraint<field_type, blocksize> >* boundObstacles,
+ std::shared_ptr<const JacobianType> constraintMatrix,
+ std::shared_ptr<const std::vector<BoxConstraint<field_type,constraintBlocksize> > > constraintObstacles)
+ : x_(x), rhs_(rhs),
+ dirichletNodes_(dirichletNodes), obstacles_(boundObstacles),
+ constraintMatrix_(constraintMatrix), constraintObstacles_(constraintObstacles),
+ jacobianCutoff_(1e-8)
+ {}
+
+ /** default destructor */
+ virtual ~LinearIPOptProblem() {}
+
+ /**@name Overloaded from TNLP */
+ //@{
+ /** Method to return some info about the nlp */
+ virtual bool get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+ Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style);
+
+ /** Method to return the bounds for my problem */
+ virtual bool get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+ Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u);
+
+ /** Method to return the starting point for the algorithm */
+ virtual bool get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+ bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+ Ipopt::Index m, bool init_lambda,
+ Ipopt::Number* lambda);
+
+ /** Method to return the objective value */
+ virtual bool eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value);
+
+ /** Method to return the gradient of the objective */
+ virtual bool eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f);
+
+ /** Method to return the constraint residuals */
+ virtual bool eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g);
+
+ /** Method to return:
+ * 1) The structure of the jacobian (if "values" is NULL)
+ * 2) The values of the jacobian (if "values" is not NULL)
+ */
+ virtual bool eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+ Ipopt::Number* values);
+
+ /** Method to return:
+ * 1) The structure of the hessian of the lagrangian (if "values" is NULL)
+ * 2) The values of the hessian of the lagrangian (if "values" is not NULL)
+ */
+ virtual bool eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+ bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+ Ipopt::Index* jCol, Ipopt::Number* values);
+
+ //@}
+
+ /** @name Solution Methods */
+ //@{
+ /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
+ virtual void finalize_solution(Ipopt::SolverReturn status,
+ Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+ Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+ Ipopt::Number obj_value,
+ const Ipopt::IpoptData* ip_data,
+ Ipopt::IpoptCalculatedQuantities* ip_cq);
+ //@}
+
+ // /////////////////////////////////
+ // Data
+ // /////////////////////////////////
+
+ //! Vector to store the solution
+ VectorType* x_;
+
+ //! The linear term of the linear energy
+ const VectorType* rhs_;
+
+ const Dune::BitSetVector<blocksize>* dirichletNodes_;
+
+ //! Vector containing the bound constraints of the variables
+ // Can stay unset when no bound constraints exist
+ std::vector<BoxConstraint<field_type, blocksize> >* obstacles_;
+
+ //! The matrix corresponding to linear constraints
+ // Can stay unset when no linear constraints exist
+ std::shared_ptr<const JacobianType> constraintMatrix_;
+
+ //! IpOpt expects constraints to be of the type g_l <= g(x) <= g_u
+ // Can stay unset when no linear constraints exist
+ std::shared_ptr<const std::vector<BoxConstraint<field_type, constraintBlocksize> > > constraintObstacles_;
+
+private:
+ const field_type jacobianCutoff_;
+ /**@name Methods to block default compiler methods.*/
+ //@{
+ LinearIPOptProblem(const LinearIPOptProblem&);
+ LinearIPOptProblem& operator=(const LinearIPOptProblem&);
+ //@}
+};
+
+// returns the size of the problem
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType, JacobianType>::
+get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+ Ipopt::Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style)
+{
+ // use the C style indexing (0-based)
+ index_style = Ipopt::TNLP::C_STYLE;
+
+ // Number of variables
+ n = x_->dim();
+
+ // Number of constraints
+ // Dirichlet conditions are actually bound inequality constraints
+ // bound constraints are handled differently by IpOpt
+ m = 0;
+
+ // hence the constraint jacobian is empty
+ nnz_jac_g = 0;
+
+ // energy is linear, thus the hessian is zero
+ nnz_h_lag = 0;
+
+ // We only need the lower left corner of the Hessian (since it is symmetric)
+ if (constraintMatrix_==nullptr)
+ return true;
+
+ assert(constraintMatrix_->N()==constraintObstacles_->size());
+
+ // only change the constraint information
+ m = constraintMatrix_->N()*constraintBlocksize;
+
+ int numJacobianNonzeros = 0;
+ for (size_t i=0; i<constraintMatrix_->N(); i++) {
+
+ const typename JacobianType::row_type& row = (*constraintMatrix_)[i];
+
+ // dim for each entry in the constraint matrix
+ typename JacobianType::row_type::ConstIterator cIt = row.begin();
+ typename JacobianType::row_type::ConstIterator cEndIt = row.end();
+ for (; cIt!=cEndIt; cIt++)
+ for (int k=0; k<constraintBlocksize; k++)
+ for (int l=0; l<blocksize; l++)
+ if (std::abs((*cIt)[k][l]) > jacobianCutoff_ )
+ numJacobianNonzeros++;
+ }
+
+ nnz_jac_g = numJacobianNonzeros;
+
+ return true;
+}
+
+
+// returns the variable bounds
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType, JacobianType>::
+get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+ Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u)
+{
+ // here, the n and m we gave IPOPT in get_nlp_info are passed back to us.
+ // If desired, we could assert to make sure they are what we think they are.
+ assert(n == (Ipopt::Index)x_->dim());
+
+ if (obstacles_) {
+
+ for (size_t i=0; i<x_->size(); i++) {
+
+ for (size_t j=0; j<blocksize; j++) {
+
+ // Simply copy the values. If they exceed a certain limit they are
+ // considered 'off'.
+ x_l[i*blocksize+j] = (*obstacles_)[i].lower(j);
+ x_u[i*blocksize+j] = (*obstacles_)[i].upper(j);
+
+ }
+
+ }
+
+ } else {
+
+ // unset all bounds
+ for (size_t i=0; i<x_->dim(); i++) {
+ x_l[i] = -std::numeric_limits<double>::max();
+ x_u[i] = std::numeric_limits<double>::max();
+ }
+
+ }
+
+ if (dirichletNodes_) {
+
+ for (size_t i=0; i<x_->size(); i++) {
+
+ for (size_t j=0; j<blocksize; j++) {
+
+ if ( (*dirichletNodes_)[i][j])
+ x_l[i*blocksize+j] = x_u[i*blocksize+j] = (*this->x_)[i][j];
+
+ }
+
+ }
+
+ }
+
+ if (constraintMatrix_==nullptr)
+ return true;
+
+ // initialize with large numbers
+ for (int i=0; i<m; i++) {
+ g_l[i] = -std::numeric_limits<field_type>::max();
+ g_u[i] = std::numeric_limits<field_type>::max();
+ }
+
+ for (size_t i=0; i<constraintObstacles_->size(); i++) {
+ for (int k=0; k<constraintBlocksize; k++) {
+
+ g_l[i*constraintBlocksize + k] = (*constraintObstacles_)[i].lower(k);
+ g_u[i*constraintBlocksize + k] = (*constraintObstacles_)[i].upper(k);
+ }
+ }
+
+ return true;
+}
+
+// returns the initial point for the problem
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+ bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+ Ipopt::Index m, bool init_lambda, Ipopt::Number* lambda)
+{
+ // Here, we assume we only have starting values for x, if you code
+ // your own NLP, you can provide starting values for the dual variables
+ // if you wish
+ assert(init_x == true);
+ assert(init_z == false);
+ assert(init_lambda == false);
+
+ // initialize to the given starting point
+ for (size_t i=0; i<x_->size(); i++)
+ for (int j=0; j<blocksize; j++)
+ x[i*blocksize+j] = (*x_)[i][j];
+
+ return true;
+}
+
+// returns the value of the objective function
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value)
+{
+ // Init return value
+ obj_value = 0;
+
+ ////////////////////////////////////
+ // Compute rhs*x
+ ////////////////////////////////////
+
+ for (size_t i=0; i<rhs_->size(); i++)
+ for (int j=0; j<blocksize; j++)
+ obj_value += x[i*blocksize+j] * (*rhs_)[i][j];
+
+ return true;
+}
+
+// return the gradient of the objective function grad_{x} f(x)
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f)
+{
+ // g = g - b
+ for (size_t i=0; i<rhs_->size(); i++)
+ for (int j=0; j<blocksize; j++)
+ grad_f[i*blocksize+j] = (*rhs_)[i][j];
+
+ return true;
+}
+
+// return the value of the constraints: g(x)
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
+{
+ if (constraintMatrix_==nullptr)
+ return true;
+
+ for (int i=0; i<m; i++)
+ g[i] = 0;
+
+ int rowCounter = 0;
+ for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+ const auto& row = (*constraintMatrix_)[rowIdx];
+ auto cIt = row.begin();
+ auto cEnd = row.end();
+
+ // constraintMatrix times x
+ for (; cIt != cEnd; cIt++)
+ for (int l=0; l<constraintBlocksize; l++)
+ for (int k=0; k<blocksize; k++)
+ if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+ g[rowCounter+l] += (*cIt)[l][k] * x[cIt.index()*blocksize+k];
+
+ // increment constraint counter
+ rowCounter += constraintBlocksize;
+ }
+
+ return true;
+
+}
+
+// return the structure or values of the jacobian
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+ Ipopt::Number* values)
+{
+ if (constraintMatrix_==nullptr)
+ return true;
+
+ int idx = 0;
+ int rowCounter = 0;
+
+ // If the values are null then IpOpt wants the sparsity pattern
+ if (values==NULL) {
+
+ for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+ const auto& row = (*constraintMatrix_)[rowIdx];
+ auto cIt = row.begin();
+ auto cEnd = row.end();
+
+ for (; cIt != cEnd; ++cIt) {
+
+ for (int l=0; l<constraintBlocksize; l++)
+ for (int k=0; k<blocksize; k++)
+ if ( std::abs((*cIt)[l][k]) > jacobianCutoff_) {
+
+ iRow[idx] = rowCounter + l;
+ jCol[idx] = cIt.index()*blocksize+k;
+ idx++;
+ }
+ }
+
+ rowCounter += constraintBlocksize;
+ }
+
+ // If the values are not null IpOpt wants the evaluated constraints gradient
+ } else {
+
+ for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+ const auto& row = (*constraintMatrix_)[rowIdx];
+ auto cIt = row.begin();
+ auto cEnd = row.end();
+
+ for (; cIt != cEnd; ++cIt) {
+
+ for (int l=0; l<constraintBlocksize; l++)
+ for (int k=0; k<blocksize; k++)
+ if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+ values[idx++] = (*cIt)[l][k];
+ }
+
+ rowCounter += constraintBlocksize;
+
+ }
+
+ }
+
+ return true;
+
+}
+
+//return the structure or values of the hessian
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+ bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+ Ipopt::Index* jCol, Ipopt::Number* values)
+{
+ return true;
+}
+
+template <class VectorType, class JacobianType>
+void LinearIPOptProblem<VectorType,JacobianType>::
+finalize_solution(Ipopt::SolverReturn status,
+ Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+ Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+ Ipopt::Number obj_value,
+ const Ipopt::IpoptData* ip_data,
+ Ipopt::IpoptCalculatedQuantities* ip_cq)
+{
+ // here is where we would store the solution to variables, or write to a file, etc
+ // so we could use the solution.
+
+ for (size_t i=0; i<x_->size(); i++)
+ for (int j=0; j<blocksize; j++)
+ (*x_)[i][j] = x[i*blocksize+j];
+
+}
+
+/** \brief Wraps the IPOpt interior-point solver for quadratic problems with
+ * linear constraints and bound constraints.
+ *
+ * The problems that can be solved are of the form
+ * \f[
+ * \min_x \frac{1}{2}x^T A x - b^T x
+ * x_l \leq x \leq x_u
+ * g_l \leq G x \leq g_u
+ * \f]
+ *
+ * \tparam JacobianType The type of the jacobian of the linear constraints
+*/
+template <class VectorType, class JacobianType>
+class LinearIPOptSolver
+ : public IterativeSolver<VectorType>, public CanIgnore<Dune::BitSetVector<VectorType::block_type::dimension> >
+{
+
+ enum {blocksize = VectorType::block_type::dimension};
+ enum {constraintBlocksize = JacobianType::block_type::rows};
+
+ typedef typename VectorType::field_type field_type;
+
+ typedef Dune::FieldMatrix<field_type, blocksize, blocksize> MatrixBlock;
+ typedef Dune::FieldVector<field_type, blocksize> VectorBlock;
+
+
+public:
+
+ /** \brief Default constructor */
+ LinearIPOptSolver () : IterativeSolver<VectorType>(1e-8, 100, NumProc::FULL),
+ rhs_(NULL),
+ obstacles_(NULL),
+ linearSolverType_(""),
+ constraintObstacles_(nullptr),
+ constraintMatrix_(nullptr)
+ {}
+
+ /** \brief Constructor for a linear problem */
+ LinearIPOptSolver (VectorType& x,
+ const VectorType& rhs,
+ NumProc::VerbosityMode verbosity=NumProc::FULL,
+ std::string linearSolverType = "")
+ : IterativeSolver<VectorType>(1e-8, 100, verbosity),
+ rhs_(&rhs), obstacles_(NULL),
+ linearSolverType_(linearSolverType),
+ constraintObstacles_(nullptr),
+ constraintMatrix_(nullptr)
+ {
+ this->x_ = &x;
+ }
+
+ void setProblem(VectorType& x,
+ const VectorType& rhs) {
+ x_ = &x;
+ rhs_ = &rhs;
+ }
+
+ void setLinearConstraints(const JacobianType& constraintMatrix,
+ const std::vector<BoxConstraint<field_type,constraintBlocksize> >& obstacles) {
+ constraintMatrix_ = Dune::stackobject_to_shared_ptr(constraintMatrix);
+ constraintObstacles_ = Dune::stackobject_to_shared_ptr(obstacles);
+ }
+
+ virtual void solve();
+
+ ///////////////////////////////////////////////////////
+ // Data
+ ///////////////////////////////////////////////////////
+
+ //! Vector to store the solution
+ VectorType* x_;
+
+ //! The linear term in the linear energy
+ const VectorType* rhs_;
+
+ //! Vector containing the bound constraints of the variables
+ // Can stay unset when no bound constraints exist
+ std::vector<BoxConstraint<field_type,blocksize> >* obstacles_;
+
+ //! The type of the linear solver to be used within IpOpt, default is ""
+ std::string linearSolverType_;
+
+ //! IpOpt expects constraints to be of the type g_l <= g(x) <= g_u
+ // Can stay unset when no linear constraints exist
+ std::shared_ptr<const std::vector<BoxConstraint<field_type,constraintBlocksize> > > constraintObstacles_;
+
+ //! The matrix corresponding to linear constraints
+ // Can stay unset when no linear constraints exist
+ std::shared_ptr<const JacobianType> constraintMatrix_;
+};
+
+template <class VectorType, class JacobianType>
+void LinearIPOptSolver<VectorType,JacobianType>::solve()
+{
+ // Create a new instance of your nlp
+ // (use a SmartPtr, not raw)
+ Ipopt::SmartPtr<Ipopt::TNLP> mynlp = new LinearIPOptProblem<VectorType,JacobianType>(x_, rhs_,
+ this->ignoreNodes_, obstacles_,
+ constraintMatrix_,constraintObstacles_);
+
+ // Create a new instance of IpoptApplication
+ // (use a SmartPtr, not raw)
+ Ipopt::SmartPtr<Ipopt::IpoptApplication> app = new Ipopt::IpoptApplication();
+
+ // Change some options
+ app->Options()->SetNumericValue("tol", this->tolerance_);
+ app->Options()->SetIntegerValue("max_iter", this->maxIterations_);
+ app->Options()->SetStringValue("mu_strategy", "adaptive");
+ if (linearSolverType_ != "")
+ app->Options()->SetStringValue("linear_solver",linearSolverType_);
+ app->Options()->SetStringValue("output_file", "ipopt.out");
+ app->Options()->SetStringValue("hessian_constant", "yes");
+
+ // We only support linear constraints
+ app->Options()->SetStringValue("jac_d_constant","yes");
+ app->Options()->SetStringValue("jac_c_constant","yes");
+
+ switch (this->verbosity_) {
+ case NumProc::QUIET:
+ app->Options()->SetIntegerValue("print_level", 0);
+ break;
+ case NumProc::REDUCED:
+ app->Options()->SetIntegerValue("print_level", 5);
+ break;
+ case NumProc::FULL:
+ app->Options()->SetIntegerValue("print_level", 12);
+ }
+
+ // Intialize the IpoptApplication and process the options
+ Ipopt::ApplicationReturnStatus status;
+ status = app->Initialize();
+ if (status != Ipopt::Solve_Succeeded)
+ DUNE_THROW(Dune::Exception, "IPOpt: Error during initialization!");
+
+ // Ask Ipopt to solve the problem
+ status = app->OptimizeTNLP(mynlp);
+
+ // Print helpful text for at least some cases
+ if (status == Ipopt::Invalid_Option)
+ DUNE_THROW(Dune::Exception, "IPOpt returned 'Invalid_Option' error!");
+
+ if (status == Ipopt::Solved_To_Acceptable_Level)
+ std::cout<<"WARNING: Desired tolerance could not be reached, but still accetable tolerance is reached.\n";
+ else if (status != Ipopt::Solve_Succeeded)
+ DUNE_THROW(Dune::Exception, "IPOpt: Error during optimization!");
+}
+
+#endif