From 51ef2200656abe1eb024d84d9c831f331476d2c6 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 10 May 2009 21:00:47 +0000
Subject: [PATCH] IPOpt solver moved from ag-common to dune-solvers
[[Imported from SVN: r2429]]
---
Makefile.am | 7 +-
configure.ac | 9 +-
dune-solvers/Makefile.am | 7 +
dune-solvers/iterationsteps/Makefile.am | 9 +
dune-solvers/norms/Makefile.am | 7 +
dune-solvers/solvers/Makefile.am | 7 +
dune-solvers/solvers/quadraticipopt.hh | 605 +++++++++++++++++++++
dune-solvers/transferoperators/Makefile.am | 10 +
m4/ipopt.m4 | 97 ++++
9 files changed, 750 insertions(+), 8 deletions(-)
create mode 100644 dune-solvers/Makefile.am
create mode 100644 dune-solvers/iterationsteps/Makefile.am
create mode 100644 dune-solvers/norms/Makefile.am
create mode 100644 dune-solvers/solvers/Makefile.am
create mode 100644 dune-solvers/solvers/quadraticipopt.hh
create mode 100644 dune-solvers/transferoperators/Makefile.am
create mode 100644 m4/ipopt.m4
diff --git a/Makefile.am b/Makefile.am
index 68eb1740..a6aa5e04 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -3,13 +3,10 @@
# we need the module file to be able to build via dunecontrol
EXTRA_DIST=dune.module
-DIST_SUBDIRS = doc src m4
-
-SUBDIRS = src m4
+SUBDIRS = src m4 dune-solvers
if BUILD_DOCS
-# TODO: set up documentation tree automatically
-#SUBDIRS += doc
+SUBDIRS += doc
endif
# don't follow the full GNU-standard
diff --git a/configure.ac b/configure.ac
index 95b8d267..d4a970ff 100644
--- a/configure.ac
+++ b/configure.ac
@@ -7,9 +7,7 @@ AC_CONFIG_SRCDIR([src/dune_solvers.cc])
AM_CONFIG_HEADER([config.h])
-# we need no more than the standard DE-stuff
-# this module depends on dune-common dune-grid dune-istl
-# this implies checking for [dune-common], [dune-grid], [dune-istl]
+# Run all checks
DUNE_CHECK_ALL
# Create a symlink dune --> . in the main module directory
@@ -23,6 +21,11 @@ LIBS="$DUNE_LIBS"
AC_CONFIG_FILES([
Makefile
src/Makefile
+ dune-solvers/Makefile
+ dune-solvers/iterationsteps/Makefile
+ dune-solvers/norms/Makefile
+ dune-solvers/solvers/Makefile
+ dune-solvers/transferoperators/Makefile
doc/Makefile
doc/doxygen/Makefile
m4/Makefile
diff --git a/dune-solvers/Makefile.am b/dune-solvers/Makefile.am
new file mode 100644
index 00000000..89ba8da5
--- /dev/null
+++ b/dune-solvers/Makefile.am
@@ -0,0 +1,7 @@
+
+SUBDIRS = iterationsteps norms solvers transferoperators
+
+dune_solversdir = $(includedir)/dune/dune-solvers
+dune_solvers_HEADERS = boxconstraint.hh numproc.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune-solvers/iterationsteps/Makefile.am b/dune-solvers/iterationsteps/Makefile.am
new file mode 100644
index 00000000..9422dcf3
--- /dev/null
+++ b/dune-solvers/iterationsteps/Makefile.am
@@ -0,0 +1,9 @@
+SUBDIRS =
+
+iterationstepsdir = $(includedir)/dune/dune-solvers/iterationsteps
+iterationsteps_HEADERS = amgstep.hh blockgsstep.hh blockgsstep.cc iterationstep.hh \
+ lineariterationstep.hh mmgstep.hh mmgstep.cc multigridstep.hh multigridstep.cc \
+ projectedblockgsstep.hh projectedblockgsstep.cc truncatedblockgsstep.hh \
+ truncatedsaddlepointgsstep.hh trustregiongsstep.cc trustregiongsstep.hh
+
+include $(top_srcdir)/am/global-rules
\ No newline at end of file
diff --git a/dune-solvers/norms/Makefile.am b/dune-solvers/norms/Makefile.am
new file mode 100644
index 00000000..2ec1e08d
--- /dev/null
+++ b/dune-solvers/norms/Makefile.am
@@ -0,0 +1,7 @@
+SUBDIRS =
+
+normsdir = $(includedir)/dune/dune-solvers/norms
+norms_HEADERS = blocknorm.hh diagnorm.hh energynorm.hh fullnorm.hh h1seminorm.hh \
+ norm.hh pnorm.hh sumnorm.hh
+
+include $(top_srcdir)/am/global-rules
\ No newline at end of file
diff --git a/dune-solvers/solvers/Makefile.am b/dune-solvers/solvers/Makefile.am
new file mode 100644
index 00000000..46b64fb0
--- /dev/null
+++ b/dune-solvers/solvers/Makefile.am
@@ -0,0 +1,7 @@
+SUBDIRS =
+
+solversdir = $(includedir)/dune/dune-solvers/solvers
+solvers_HEADERS = cgsolver.cc cgsolver.hh iterativesolver.cc iterativesolver.hh \
+ loopsolver.cc loopsolver.hh solver.hh tcgsolver.cc tcgsolver.hh
+
+include $(top_srcdir)/am/global-rules
\ No newline at end of file
diff --git a/dune-solvers/solvers/quadraticipopt.hh b/dune-solvers/solvers/quadraticipopt.hh
new file mode 100644
index 00000000..6cea09a3
--- /dev/null
+++ b/dune-solvers/solvers/quadraticipopt.hh
@@ -0,0 +1,605 @@
+#ifndef DUNE_QUADRATIC_IPOPT_SOLVER_HH
+#define DUNE_QUADRATIC_IPOPT_SOLVER_HH
+
+/** \file
+ \brief A wrapper for the IPOpt solver for quadratic optimization problem
+ with box constraints.
+*/
+#include <algorithm>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune-solvers/boxconstraint.hh>
+#include <dune-solvers/solvers/iterativesolver.hh>
+
+/** \todo Undef those because they clash. I think this is an IPOPT bug.
+ (their config.h is included by a header). */
+#undef PACKAGE_BUGREPORT
+#undef PACKAGE_NAME
+#undef PACKAGE_STRING
+#undef PACKAGE_TARNAME
+#undef PACKAGE_VERSION
+
+#include "coin/IpTNLP.hpp"
+#include "coin/IpIpoptApplication.hpp"
+
+template <class MatrixType, class VectorType>
+class QuadraticIPOptProblem : public Ipopt::TNLP
+{
+ enum {blocksize = VectorType::block_type::size};
+
+ typedef typename VectorType::field_type field_type;
+
+public:
+ /** default constructor */
+ QuadraticIPOptProblem(const MatrixType* hessian, VectorType* x, const VectorType* rhs,
+ const Dune::BitSetVector<blocksize>* dirichletNodes,
+ std::vector<BoxConstraint<field_type, blocksize> >* obstacles)
+ : hessian_(hessian), x_(x), rhs_(rhs),
+ dirichletNodes_(dirichletNodes), obstacles_(obstacles)
+ {}
+
+ /** default destructor */
+ virtual ~QuadraticIPOptProblem() {};
+
+ /**@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
+ // /////////////////////////////////
+
+ const MatrixType* hessian_;
+
+ VectorType* x_;
+
+ const VectorType* rhs_;
+
+ const Dune::BitSetVector<blocksize>* dirichletNodes_;
+
+ std::vector<BoxConstraint<field_type, blocksize> >* obstacles_;
+
+private:
+ /**@name Methods to block default compiler methods.
+ */
+ //@{
+ // QuadraticIPOptProblem();
+ QuadraticIPOptProblem(const QuadraticIPOptProblem&);
+ QuadraticIPOptProblem& operator=(const QuadraticIPOptProblem&);
+ //@}
+};
+
+// returns the size of the problem
+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)
+{
+ // The problem described in HS071_NLP.hpp has 4 variables, x[0] through x[3]
+ n = x_->dim();
+
+ // No equality constraints: Dirichlet conditions are actually inequality constraints
+ m = 0;
+
+ // hence the constraint jacobian is empty
+ nnz_jac_g = 0;
+
+ // We only need the lower left corner of the Hessian (since it is symmetric)
+ if (!hessian_)
+ DUNE_THROW(SolverError, "No Hessian matrix has been supplied!");
+
+ nnz_h_lag = 0;
+
+ for (int rowIdx=0; rowIdx<hessian_->N(); rowIdx++) {
+
+ const typename MatrixType::row_type& row = (*hessian_)[rowIdx];
+
+ typename MatrixType::row_type::ConstIterator cIt = row.begin();
+ typename MatrixType::row_type::ConstIterator cEndIt = row.end();
+
+ for (; cIt!=cEndIt; ++cIt) {
+
+ if (rowIdx > cIt.index()) {
+
+ nnz_h_lag += blocksize*blocksize;
+
+ } else if (rowIdx == cIt.index()) {
+
+ nnz_h_lag += ((blocksize+1)*blocksize) / 2;
+
+ }
+
+ }
+
+ }
+
+ // use the C style indexing (0-based)
+ index_style = Ipopt::TNLP::C_STYLE;
+
+ return true;
+}
+
+
+// returns the variable bounds
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptProblem<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)
+{
+ // 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 == x_->dim());
+ assert(m == 0);
+
+ if (obstacles_) {
+
+ for (int 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 (int 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 (int 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];
+
+ }
+
+ }
+
+ }
+
+ return true;
+}
+
+// returns the initial point for the problem
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptProblem<MatrixType,VectorType>::
+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 (int 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 MatrixType, class VectorType>
+bool QuadraticIPOptProblem<MatrixType,VectorType>::
+eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value)
+{
+// std::cout << "x:" << std::endl;
+// for (int i=0; i<n; i++)
+// std::cout << x[i] << std::endl;
+
+ // Init return value
+ obj_value = 0;
+
+ ////////////////////////////////////
+ // Compute x^T*A*x
+ ////////////////////////////////////
+
+ for (int rowIdx=0; rowIdx<hessian_->N(); rowIdx++) {
+
+ const typename MatrixType::row_type& row = (*hessian_)[rowIdx];
+
+ typename MatrixType::row_type::ConstIterator cIt = row.begin();
+ typename MatrixType::row_type::ConstIterator cEndIt = row.end();
+
+ for (; cIt!=cEndIt; ++cIt) {
+
+ // Multiply block
+ for (int i=0; i<blocksize; i++)
+ for (int j=0; j<blocksize; j++)
+ obj_value += x[rowIdx*blocksize+i] * x[cIt.index()*blocksize+j] * (*cIt)[i][j];
+
+ }
+
+ }
+
+ //std::cout << "bilinear form: " << *f << std::endl;
+
+ // because J = 1/2 * A(x,x) - b(x)
+ obj_value *= 0.5;
+
+ // -b(x)
+ for (int i=0; i<rhs_->size(); i++)
+ for (int j=0; j<blocksize; j++)
+ obj_value -= x[i*blocksize+j] * (*rhs_)[i][j];
+
+ //std::cout << "IPOPT Energy: " << obj_value << std::endl;
+ //exit(0);
+
+ return true;
+}
+
+// return the gradient of the objective function grad_{x} f(x)
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptProblem<MatrixType,VectorType>::
+eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f)
+{
+ //std::cout << "### eval_grad_f ###" << std::endl;
+
+ // \nabla J = A(x,.) - b(x)
+ for (int i=0; i<n; i++)
+ grad_f[i] = 0;
+
+ for (int i=0; i<hessian_->N(); i++) {
+
+ const typename MatrixType::row_type& row = (*hessian_)[i];
+
+ typename MatrixType::row_type::ConstIterator cIt = row.begin();
+ typename MatrixType::row_type::ConstIterator cEndIt = row.end();
+
+ for (; cIt!=cEndIt; ++cIt) {
+
+ for (int j=0; j<blocksize; j++)
+ for (int k=0; k<blocksize; k++)
+ grad_f[i*blocksize+j] += (*cIt)[j][k] * x[cIt.index()*blocksize+k];
+
+ }
+
+ }
+
+ // g = g - b
+ for (int i=0; i<rhs_->size(); i++)
+ for (int j=0; j<blocksize; j++)
+ grad_f[i*blocksize+j] -= (*rhs_)[i][j];
+
+// for (int i=0; i<n; i++)
+// printf("%d: x = %g, rhs = %g, g = %g\n", i, x[i], (*rhs_)[i/blocksize][i%blocksize], grad_f[i]);
+ return true;
+}
+
+// return the value of the constraints: g(x)
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptProblem<MatrixType,VectorType>::
+eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
+{
+ return true;
+}
+
+// return the structure or values of the jacobian
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptProblem<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)
+{
+ return true;
+}
+
+//return the structure or values of the hessian
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptProblem<MatrixType,VectorType>::
+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)
+{
+ if (values == NULL) {
+
+ // return the structure. This is a symmetric matrix, fill the lower left
+ // triangle only.
+
+ int idx = 0;
+
+ // Loop over all matrix rows
+ for (int rowIdx=0; rowIdx<hessian_->N(); rowIdx++) {
+
+ const typename MatrixType::row_type& row = (*hessian_)[rowIdx];
+
+ typename MatrixType::row_type::ConstIterator cIt = row.begin();
+ typename MatrixType::row_type::ConstIterator cEndIt = row.end();
+
+ // Loop over all columns in this row
+ for (; cIt!=cEndIt; ++cIt) {
+
+ // Current entry on the lower left triangular matrix and not on the diagonal
+ // --> copy it completely
+ if (rowIdx > cIt.index()) {
+
+ for (int i=0; i<blocksize; i++) {
+
+ for (int j=0; j<blocksize; j++) {
+
+ iRow[idx] = rowIdx*blocksize + i;
+ jCol[idx] = cIt.index()*blocksize + j;
+ idx++;
+
+ }
+
+ }
+
+ } else if (rowIdx == cIt.index()) {
+
+ // Current entry is on the diagonal
+ // --> copy only its lower triangular half
+ for (int i=0; i<blocksize; i++) {
+
+ for (int j=0; j<i+1; j++) {
+
+ iRow[idx] = rowIdx*blocksize + i;
+ jCol[idx] = cIt.index()*blocksize + j;
+ idx++;
+
+ }
+
+ }
+
+ }
+
+ }
+
+ }
+
+ assert(idx == nele_hess);
+// for (int i=0; i<nele_hess; i++)
+// printf("%d %d \n", iRow[i], jCol[i]);
+
+ } else {
+
+ // return the values. This is a symmetric matrix, fill the lower left
+ // triangle only
+
+ int idx = 0;
+
+ // Loop over all matrix rows
+ for (int rowIdx=0; rowIdx<hessian_->N(); rowIdx++) {
+
+ const typename MatrixType::row_type& row = (*hessian_)[rowIdx];
+
+ typedef typename MatrixType::row_type::ConstIterator ColumnIterator;
+
+ ColumnIterator cIt = row.begin();
+ ColumnIterator cEndIt = row.end();
+
+ // Loop over all columns in this row
+ for (; cIt!=cEndIt; ++cIt) {
+
+ // Current entry on the lower left triangular matrix and not on the diagonal
+ // --> copy it completely
+ if (rowIdx > cIt.index()) {
+
+ for (int i=0; i<blocksize; i++)
+ for (int j=0; j<blocksize; j++) {
+ //printf("%d %d %g\n", rowIdx*blocksize + i, cIt.index()*blocksize+j,(*cIt)[i][j]);
+ values[idx++] = obj_factor * (*cIt)[i][j];
+ }
+
+ } else if (rowIdx == cIt.index()) {
+
+ // Current entry is on the diagonal
+ // --> copy only its lower triangular half
+ for (int i=0; i<blocksize; i++)
+ for (int j=0; j<i+1; j++) {
+ //printf("%d %d %g\n", rowIdx*blocksize + i, cIt.index()*blocksize+j,(*cIt)[i][j]);
+ values[idx++] = obj_factor * (*cIt)[i][j];
+ }
+ }
+
+ }
+
+ }
+
+ assert(idx == nele_hess);
+
+// for (int i=0; i<nele_hess; i++)
+// printf("%g\n", values[i]);
+
+ }
+
+ return true;
+}
+
+template <class MatrixType, class VectorType>
+void QuadraticIPOptProblem<MatrixType,VectorType>::
+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 (int i=0; i<x_->size(); i++)
+ for (int j=0; j<blocksize; j++)
+ (*x_)[i][j] = x[i*blocksize+j];
+
+}
+
+
+template <class MatrixType, class VectorType>
+class QuadraticIPOptSolver : public IterativeSolver<VectorType>
+{
+
+ enum {blocksize = VectorType::block_type::size};
+
+ 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 */
+ QuadraticIPOptSolver () : IterativeSolver<VectorType>(1e-8, 100, NumProc::FULL),
+ hessian_(NULL), rhs_(NULL), ignoreNodes_(NULL),
+ obstacles_(NULL)
+ {}
+
+ /** \brief Constructor for a linear problem */
+ QuadraticIPOptSolver (const MatrixType& hessian,
+ VectorType& x,
+ const VectorType& rhs)
+ : IterativeSolver<VectorType>(1e-8, 100, NumProc::FULL),
+ hessian_(&hessian), rhs_(&rhs), ignoreNodes_(NULL), obstacles_(NULL)
+ {
+ this->x_ = &x;
+ }
+
+ void setProblem(const MatrixType& hessian,
+ VectorType& x,
+ const VectorType& rhs) {
+ hessian_ = &hessian;
+ x_ = &x;
+ rhs_ = &rhs;
+ }
+
+ virtual void solve();
+
+ ///////////////////////////////////////////////////////
+ // Data
+ ///////////////////////////////////////////////////////
+
+ const MatrixType* hessian_;
+
+ VectorType* x_;
+
+ const VectorType* rhs_;
+
+ const Dune::BitSetVector<blocksize>* ignoreNodes_;
+
+ std::vector<BoxConstraint<field_type,blocksize> >* obstacles_;
+
+};
+
+template <class MatrixType, class VectorType>
+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_,
+ ignoreNodes_, obstacles_);
+
+ // 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");
+ app->Options()->SetStringValue("output_file", "ipopt.out");
+ switch (this->verbosity_) {
+ case NumProc::QUIET:
+ app->Options()->SetIntegerValue("print_level", -2);
+ 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);
+
+ // As the SmartPtrs go out of scope, the reference count
+ // will be decremented and the objects will automatically
+ // be deleted.
+
+}
+
+#endif
diff --git a/dune-solvers/transferoperators/Makefile.am b/dune-solvers/transferoperators/Makefile.am
new file mode 100644
index 00000000..a486dab3
--- /dev/null
+++ b/dune-solvers/transferoperators/Makefile.am
@@ -0,0 +1,10 @@
+SUBDIRS =
+
+transferoperatorsdir = $(includedir)/dune/dune-solvers/transferoperators
+transferoperators_HEADERS = compressedmultigridtransfer.hh multigridtransfer.hh \
+ truncatedcompressedmgtransfer.hh truncatedmgtransfer.hh densemultigridtransfer.hh \
+ mandelobsrestrictor.cc obstaclerestrictor.hh truncateddensemgtransfer.cc \
+ genericmultigridtransfer.hh mandelobsrestrictor.hh truncatedcompressedmgtransfer.cc \
+ truncateddensemgtransfer.hh
+
+include $(top_srcdir)/am/global-rules
\ No newline at end of file
diff --git a/m4/ipopt.m4 b/m4/ipopt.m4
new file mode 100644
index 00000000..44b5007f
--- /dev/null
+++ b/m4/ipopt.m4
@@ -0,0 +1,97 @@
+# searches for ipopt headers and libs
+
+AC_DEFUN([DUNE_IPOPT],[
+ AC_REQUIRE([AC_PROG_CXX])
+
+ AC_ARG_WITH(ipopt,
+ AC_HELP_STRING([--with-ipopt=PATH],[directory with IPOpt inside]))
+
+# store values
+ac_save_LDFLAGS="$LDFLAGS"
+ac_save_CPPFLAGS="$CPPFLAGS"
+ac_save_LIBS="$LIBS"
+
+## do nothing if --without-amiramesh is used
+if test x$with_ipopt != xno ; then
+
+# is --with-ipopt=bla used?
+if test "x$with_ipopt" != x ; then
+ if ! test -d $with_ipopt; then
+ AC_MSG_WARN([IPOpt directory $with_ipopt does not exist])
+ else
+ # expand tilde / other stuff
+ IPOPTROOT=`cd $with_ipopt && pwd`
+ fi
+fi
+if test "x$IPOPTROOT" = x; then
+ # use some default value...
+ IPOPTROOT="/usr/local/ipopt"
+fi
+
+IPOPT_LIB_PATH="$IPOPTROOT/lib"
+IPOPT_INCLUDE_PATH="$IPOPTROOT/include"
+
+LDFLAGS="$LDFLAGS -L$IPOPT_LIB_PATH"
+CPPFLAGS="$CPPFLAGS -I$IPOPT_INCLUDE_PATH"
+
+AC_LANG_PUSH([C++])
+
+# check for header
+AC_CHECK_HEADER([coin/IpIpoptApplication.hpp],
+ [IPOPT_CPPFLAGS="-I$IPOPT_INCLUDE_PATH"
+ HAVE_IPOPT="1"],
+ AC_MSG_WARN([IpIpoptApplication.hpp not found in $IPOPT_INCLUDE_PATH/coin]))
+
+CPPFLAGS="$IPOPT_CPPFLAGS"
+
+# if header is found...
+if test x$HAVE_IPOPT = x1 ; then
+ LIBS="$LIBS -lipopt -llapack -lblas -lgfortran"
+
+ AC_LINK_IFELSE(AC_LANG_PROGRAM(
+ [#include "coin/IpIpoptApplication.hpp"],
+ [Ipopt::SmartPtr<Ipopt::IpoptApplication> app = new Ipopt::IpoptApplication();]),
+ [IPOPT_LIBS="-lipopt -llapack -lblas -lgfortran"
+ IPOPT_LDFLAGS="-L$IPOPT_LIB_PATH"
+ LIBS="$LIBS $IPOPT_LIBS"],
+ [HAVE_IPOPT="0"
+ AC_MSG_WARN(IPOpt not found!)])
+fi
+
+AC_LANG_POP([C++])
+
+## end of ipopt check (--without wasn't set)
+fi
+
+with_ipopt="no"
+# survived all tests?
+if test x$HAVE_IPOPT = x1 ; then
+ AC_SUBST(IPOPT_LIBS, $IPOPT_LIBS)
+ AC_SUBST(IPOPT_LDFLAGS, $IPOPT_LDFLAGS)
+ AC_SUBST(IPOPT_CPPFLAGS, $IPOPT_CPPFLAGS)
+ AC_DEFINE(HAVE_IPOPT, 1, [Define to 1 if IPOpt library is found])
+
+ # add to global list
+ DUNE_PKG_LDFLAGS="$DUNE_PKG_LDFLAGS $IPOPT_LDFLAGS"
+ DUNE_PKG_LIBS="$DUNE_PKG_LIBS $IPOPT_LIBS"
+ DUNE_PKG_CPPFLAGS="$DUNE_PKG_CPPFLAGS $IPOPT_CPPFLAGS"
+
+ # set variable for summary
+ with_ipopt="yes"
+else
+ AC_SUBST(IPOPT_LIBS, "")
+ AC_SUBST(IPOPT_LDFLAGS, "")
+ AC_SUBST(IPOPT_CPPFLAGS, "")
+fi
+
+# also tell automake
+AM_CONDITIONAL(IPOPT, test x$HAVE_IPOPT = x1)
+
+# reset old values
+LIBS="$ac_save_LIBS"
+CPPFLAGS="$ac_save_CPPFLAGS"
+LDFLAGS="$ac_save_LDFLAGS"
+
+DUNE_ADD_SUMMARY_ENTRY([IPOpt],[$with_ipopt])
+
+])
--
GitLab