diff --git a/dune/solvers/solvers/quadraticipopt.hh b/dune/solvers/solvers/quadraticipopt.hh index 2307d16b312940a9801869e59aae74d9ffd2728f..6a9ccdf626dfeb03e00e21d616a2856b63563732 100644 --- a/dune/solvers/solvers/quadraticipopt.hh +++ b/dune/solvers/solvers/quadraticipopt.hh @@ -13,6 +13,7 @@ #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> @@ -26,11 +27,14 @@ /** \brief Implementation class used to pipe quadratic problems to the IPOpt interior-point solver + + \tparam JacobianType The type of the jacobian of the linear constraints */ -template <class MatrixType, class VectorType> +template <class MatrixType, class VectorType, class JacobianType=MatrixType> class QuadraticIPOptProblem : public Ipopt::TNLP { enum {blocksize = VectorType::block_type::dimension}; + enum {constraintBlocksize = JacobianType::block_type::rows}; typedef typename VectorType::field_type field_type; @@ -41,15 +45,15 @@ public: std::vector<BoxConstraint<field_type, blocksize> >* obstacles) : hessian_(hessian), x_(x), rhs_(rhs), dirichletNodes_(dirichletNodes), obstacles_(obstacles), - constraintMatrix_(NULL), constraintObstacles_(NULL), jacobianCutoff_(1e-8) + constraintMatrix_(nullptr), constraintObstacles_(nullptr), jacobianCutoff_(1e-8) {} /** Setup problem full problem */ QuadraticIPOptProblem(const MatrixType* hessian, VectorType* x, const VectorType* rhs, const Dune::BitSetVector<blocksize>* dirichletNodes, std::vector<BoxConstraint<field_type, blocksize> >* boundObstacles, - const MatrixType* constraintMatrix, - std::vector<BoxConstraint<field_type, blocksize> >* constraintObstacles) + std::shared_ptr<const JacobianType> constraintMatrix, + std::shared_ptr<const std::vector<BoxConstraint<field_type,constraintBlocksize> > > constraintObstacles) : hessian_(hessian), x_(x), rhs_(rhs), dirichletNodes_(dirichletNodes), obstacles_(boundObstacles), constraintMatrix_(constraintMatrix), constraintObstacles_(constraintObstacles), @@ -135,11 +139,11 @@ public: //! The matrix corresponding to linear constraints // Can stay unset when no linear constraints exist - const MatrixType* constraintMatrix_; + 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::vector<BoxConstraint<field_type, blocksize> >* constraintObstacles_; + std::shared_ptr<const std::vector<BoxConstraint<field_type, constraintBlocksize> > > constraintObstacles_; private: const field_type jacobianCutoff_; @@ -151,8 +155,8 @@ private: }; // returns the size of the problem -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,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) { @@ -200,24 +204,24 @@ get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g, } // We only need the lower left corner of the Hessian (since it is symmetric) - if (!constraintMatrix_) + if (constraintMatrix_==nullptr) return true; assert(constraintMatrix_->N()==constraintObstacles_->size()); // only change the constraint information - m = constraintMatrix_->N()*blocksize; + m = constraintMatrix_->N()*constraintBlocksize; int numJacobianNonzeros = 0; for (size_t i=0; i<constraintMatrix_->N(); i++) { - const typename MatrixType::row_type& row = (*constraintMatrix_)[i]; + const typename JacobianType::row_type& row = (*constraintMatrix_)[i]; // dim for each entry in the constraint matrix - typename MatrixType::row_type::ConstIterator cIt = row.begin(); - typename MatrixType::row_type::ConstIterator cEndIt = row.end(); + 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<blocksize; k++) + for (int k=0; k<constraintBlocksize; k++) for (int l=0; l<blocksize; l++) if (std::abs((*cIt)[k][l]) > jacobianCutoff_ ) numJacobianNonzeros++; @@ -230,8 +234,8 @@ get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g, // returns the variable bounds -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,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) { @@ -279,7 +283,7 @@ get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u, } - if (!constraintMatrix_) + if (constraintMatrix_==nullptr) return true; // initialize with large numbers @@ -289,10 +293,10 @@ get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u, } for (size_t i=0; i<constraintObstacles_->size(); i++) { - for (int k=0; k<blocksize; k++) { + for (int k=0; k<constraintBlocksize; k++) { - g_l[i*blocksize + k] = (*constraintObstacles_)[i].lower(k); - g_u[i*blocksize + k] = (*constraintObstacles_)[i].upper(k); + g_l[i*constraintBlocksize + k] = (*constraintObstacles_)[i].lower(k); + g_u[i*constraintBlocksize + k] = (*constraintObstacles_)[i].upper(k); } } @@ -300,8 +304,8 @@ get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u, } // returns the initial point for the problem -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,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) @@ -322,8 +326,8 @@ get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x, } // returns the value of the objective function -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,VectorType,JacobianType>:: eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value) { // std::cout << "x:" << std::endl; @@ -372,8 +376,8 @@ eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_va } // return the gradient of the objective function grad_{x} f(x) -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,VectorType,JacobianType>:: eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f) { //std::cout << "### eval_grad_f ###" << std::endl; @@ -410,11 +414,11 @@ eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* g } // return the value of the constraints: g(x) -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,VectorType,JacobianType>:: eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g) { - if (!constraintMatrix_) + if (constraintMatrix_==nullptr) return true; for (int i=0; i<m; i++) @@ -429,13 +433,13 @@ eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt // constraintMatrix times x for (; cIt != cEnd; cIt++) - for (int l=0; l<blocksize; l++) + 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 += blocksize; + rowCounter += constraintBlocksize; } return true; @@ -443,13 +447,13 @@ eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt } // return the structure or values of the jacobian -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,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_) + if (constraintMatrix_==nullptr) return true; int idx = 0; @@ -466,7 +470,7 @@ eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, for (; cIt != cEnd; ++cIt) { - for (int l=0; l<blocksize; l++) + for (int l=0; l<constraintBlocksize; l++) for (int k=0; k<blocksize; k++) if ( std::abs((*cIt)[l][k]) > jacobianCutoff_) { @@ -476,7 +480,7 @@ eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, } } - rowCounter += blocksize; + rowCounter += constraintBlocksize; } // If the values are not null IpOpt wants the evaluated constraints gradient @@ -490,13 +494,13 @@ eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, for (; cIt != cEnd; ++cIt) { - for (int l=0; l<blocksize; l++) + 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 += blocksize; + rowCounter += constraintBlocksize; } @@ -507,8 +511,8 @@ eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, } //return the structure or values of the hessian -template <class MatrixType, class VectorType> -bool QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +bool QuadraticIPOptProblem<MatrixType,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, @@ -590,8 +594,8 @@ eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x, return true; } -template <class MatrixType, class VectorType> -void QuadraticIPOptProblem<MatrixType,VectorType>:: +template <class MatrixType, class VectorType, class JacobianType> +void QuadraticIPOptProblem<MatrixType,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, @@ -617,13 +621,16 @@ finalize_solution(Ipopt::SolverReturn status, * 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 MatrixType, class VectorType> +template <class MatrixType, class VectorType, class JacobianType=MatrixType> class QuadraticIPOptSolver : 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; @@ -638,7 +645,7 @@ public: hessian_(NULL), rhs_(NULL), obstacles_(NULL), linearSolverType_("ma27"), - constraintObstacles_(NULL),constraintMatrix_(NULL) + constraintObstacles_(nullptr),constraintMatrix_(nullptr) {} /** \brief Constructor for a linear problem */ @@ -650,7 +657,7 @@ public: : IterativeSolver<VectorType>(1e-8, 100, verbosity), hessian_(&hessian), rhs_(&rhs), obstacles_(NULL), linearSolverType_(linearSolverType), - constraintMatrix_(NULL),constraintObstacles_(NULL) + constraintMatrix_(nullptr),constraintObstacles_(nullptr) { this->x_ = &x; } @@ -663,10 +670,10 @@ public: rhs_ = &rhs; } - void setLinearConstraints(const MatrixType& constraintMatrix, - const std::vector<BoxConstraint<field_type,blocksize> >& obstacles) { - constraintMatrix_ = &constraintMatrix; - constraintObstacles_ = &obstacles; + 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(); @@ -693,19 +700,19 @@ public: //! IpOpt expects constraints to be of the type g_l <= g(x) <= g_u // Can stay unset when no linear constraints exist - std::vector<BoxConstraint<field_type,blocksize> >* constraintObstacles_; + 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 - const MatrixType* constraintMatrix_; + std::shared_ptr<const JacobianType> constraintMatrix_; }; -template <class MatrixType, class VectorType> -void QuadraticIPOptSolver<MatrixType,VectorType>::solve() +template <class MatrixType, class VectorType, class JacobianType> +void QuadraticIPOptSolver<MatrixType,VectorType,JacobianType>::solve() { // Create a new instance of your nlp // (use a SmartPtr, not raw) - Ipopt::SmartPtr<Ipopt::TNLP> mynlp = new QuadraticIPOptProblem<MatrixType,VectorType>(hessian_, x_, rhs_, + Ipopt::SmartPtr<Ipopt::TNLP> mynlp = new QuadraticIPOptProblem<MatrixType,VectorType,JacobianType>(hessian_, x_, rhs_, this->ignoreNodes_, obstacles_, constraintMatrix_,constraintObstacles_);