diff --git a/dune/solvers/solvers/quadraticipopt.hh b/dune/solvers/solvers/quadraticipopt.hh
index 621e2928d1b045f845214e560fdb49f9bf4fddb1..76bb1d7056dd1a10dfbc096c93707db2f7ec88ec 100644
--- a/dune/solvers/solvers/quadraticipopt.hh
+++ b/dune/solvers/solvers/quadraticipopt.hh
@@ -129,7 +129,7 @@ private:
template <class MatrixType, class VectorType>
bool QuadraticIPOptProblem<MatrixType,VectorType>::
get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
- Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style)
+ Ipopt::Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style)
{
// Number of variables
n = x_->dim();
@@ -186,7 +186,7 @@ get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_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());
- assert(m == 0);
+ //assert(m == 0);
if (obstacles_) {
@@ -461,9 +461,234 @@ finalize_solution(Ipopt::SolverReturn status,
}
+/** \brief Implementation class used to pipe quadratic problems to
+ the IPOpt interior-point solver
+*/
+template <class MatrixType, class VectorType>
+class QuadraticIPOptObstacleProblem : public QuadraticIPOptProblem<MatrixType, VectorType>
+{
+ enum {blocksize = VectorType::block_type::dimension};
+
+ typedef typename VectorType::field_type field_type;
+ typedef QuadraticIPOptProblem<MatrixType, VectorType> Base;
+
+public:
+ /** default constructor */
+ QuadraticIPOptObstacleProblem(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)
+ : Base(hessian,x,rhs,dirichletNodes,boundObstacles),
+ constraintMatrix_(constraintMatrix), constraintObstacles_(constraintObstacles),
+ jacobianCutoff_(1e-8)
+ {}
+
+ /** default destructor */
+ virtual ~QuadraticIPOptObstacleProblem() {};
+
+ /**@name Overloaded from QuadraticIPOptProblem */
+ //@{
+ /** 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, Ipopt::TNLP::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 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);
+ //@}
+
+
+ // /////////////////////////////////
+ // Data
+ // /////////////////////////////////
+
+ const MatrixType* constraintMatrix_;
+
+ std::vector<BoxConstraint<field_type, blocksize> >* constraintObstacles_;
+
+private:
+ const field_type jacobianCutoff_;
+ /**@name Methods to block default compiler methods.
+ */
+ //@{
+ // QuadraticIPOptObstacleProblem();
+ QuadraticIPOptObstacleProblem(const QuadraticIPOptObstacleProblem&);
+ QuadraticIPOptObstacleProblem& operator=(const QuadraticIPOptObstacleProblem&);
+ //@}
+};
+
+// returns the size of the problem
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::
+get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+ Ipopt::Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style)
+{
+ Base::get_nlp_info(n,m,nnz_jac_g,nnz_h_lag,index_style);
+
+ // We only need the lower left corner of the Hessian (since it is symmetric)
+ if (!constraintMatrix_)
+ DUNE_THROW(SolverError, "No constraint matrix has been supplied!");
+ assert(constraintMatrix_->N()==constraintObstacles_->size());
+
+ // only change the constraint information
+ m = constraintMatrix_->N()*blocksize;
+
+ int numJacobianNonzeros = 0;
+ for (size_t i=0; i<constraintMatrix_->N(); i++) {
+
+ const auto& row = (*constraintMatrix_)[i];
+
+ // dim for entry in the constraint matrix
+ auto cEnd = row.end();
+ for (auto cIt = row.begin(); cIt!=cEnd; cIt++)
+ for (int k=0; k<blocksize; 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 MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::
+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)
+{
+ Base::get_bounds_info(n,x_l, x_u, m, g_l ,g_u);
+
+ // 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<blocksize; k++) {
+
+ g_l[i*blocksize + k] = (*constraintObstacles_)[i].lower(k);
+ g_u[i*blocksize + k] = (*constraintObstacles_)[i].upper(k);
+ }
+ }
+
+ return true;
+}
+
+// Evaluate the constraints
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
+{
+
+ 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<blocksize; 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;
+ }
+
+ return true;
+}
+
+
+// Evaluate gradient of the constraints
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::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)
+{
+
+ 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<blocksize; 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 += blocksize;
+ }
+
+ // 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<blocksize; l++)
+ for (int k=0; k<blocksize; k++)
+ if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+ values[idx++] = (*cIt)[l][k];
+ }
+
+ rowCounter += blocksize;
+
+ }
+
+ }
+
+ return true;
+}
/** \brief Wraps the IPOpt interior-point solver for quadratic problems with
- bound constraints.
+ * 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]
*/
template <class MatrixType, class VectorType>
class QuadraticIPOptSolver
@@ -484,7 +709,8 @@ public:
QuadraticIPOptSolver () : IterativeSolver<VectorType>(1e-8, 100, NumProc::FULL),
hessian_(NULL), rhs_(NULL),
obstacles_(NULL),
- linearSolverType_("ma27")
+ linearSolverType_("ma27"),
+ constraintMatrix_(NULL),constraintObstacles_(NULL)
{}
/** \brief Constructor for a linear problem */
@@ -495,7 +721,8 @@ public:
std::string linearSolverType = "ma27")
: IterativeSolver<VectorType>(1e-8, 100, verbosity),
hessian_(&hessian), rhs_(&rhs), obstacles_(NULL),
- linearSolverType_(linearSolverType)
+ linearSolverType_(linearSolverType),
+ constraintMatrix_(NULL),constraintObstacles_(NULL)
{
this->x_ = &x;
}
@@ -508,22 +735,41 @@ public:
rhs_ = &rhs;
}
+ void setLinearConstraints(const MatrixType& constraintMatrix,
+ const std::vector<BoxConstraint<field_type,blocksize> >& obstacles) {
+ constraintMatrix_ = &constraintMatrix;
+ constraintObstacles_ = &obstacles;
+ }
+
virtual void solve();
///////////////////////////////////////////////////////
// Data
///////////////////////////////////////////////////////
+ //! The quadratic term in the quadratic energy, is assumed to be symmetric
const MatrixType* hessian_;
+ //! Vector to store the solution
VectorType* x_;
+ //! The linear term in the quadratic 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 ma27
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::vector<BoxConstraint<field_type,blocksize> >* constraintObstacles_;
+
+ //! The matrix corresponding to linear constraints
+ // Can stay unset when no linear constraints exist
+ const MatrixType* constraintMatrix_;
};
template <class MatrixType, class VectorType>
@@ -531,8 +777,13 @@ void QuadraticIPOptSolver<MatrixType,VectorType>::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;
+ if (constraintMatrix_) {
+ mynlp= new QuadraticIPOptObstacleProblem<MatrixType,VectorType>(hessian_, x_, rhs_,
+ this->ignoreNodes_, obstacles_,
+ constraintMatrix_,constraintObstacles_);
+ } else
+ mynlp= new QuadraticIPOptProblem<MatrixType,VectorType>(hessian_, x_, rhs_,
this->ignoreNodes_, obstacles_);
// Create a new instance of IpoptApplication
@@ -546,6 +797,13 @@ void QuadraticIPOptSolver<MatrixType,VectorType>::solve()
app->Options()->SetStringValue("linear_solver",linearSolverType_);
app->Options()->SetStringValue("output_file", "ipopt.out");
app->Options()->SetStringValue("hessian_constant", "yes");
+
+ // We only support linear constraints
+ if (constraintMatrix_) {
+ 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);