Newton-Raphson with Known Types

Include

#include "pressio/solvers_linear.hpp"
#include "pressio/solvers_nonlinear.hpp"
#include <Eigen/Core>

Problem class

struct MySystem
{
  using scalar_type = double;
  using state_type = Eigen::VectorXd;
  using residual_type = state_type;
  using jacobian_type = Eigen::SparseMatrix<scalar_type>;

  residual_type createResidual() const {
    return residual_type(2);
  }

  jacobian_type createJacobian() const {
    return jacobian_type(2, 2);
  }

  void residual(const state_type& x,
                residual_type& res) const
  {
    res(0) =  x(0)*x(0)*x(0) + x(1) - 1.0;
    res(1) = -x(0) + x(1)*x(1)*x(1) + 1.0;
  }

  void jacobian(const state_type& x, jacobian_type& jac) const {
    jac.coeffRef(0, 0) = 3.0*x(0)*x(0);
    jac.coeffRef(0, 1) =  1.0;
    jac.coeffRef(1, 0) = -1.0;
    jac.coeffRef(1, 1) = 3.0*x(1)*x(1);
  }
};

Main

int main()
{
  namespace plog   = pressio::log;
  namespace pls    = pressio::linearsolvers;
  namespace pnonls = pressio::nonlinearsolvers;

  plog::initialize(pressio::logto::terminal);
  plog::setVerbosity({plog::level::info});

  using problem_t  = MySystem;
  problem_t problemObj;

  using state_t    = problem_t::state_type;
  state_t y(2);
  y(0) = 0.001; y(1) = -0.1;

  // linear solver
  using jacobian_t = problem_t::jacobian_type;
  using lin_solver_t = pls::Solver<pls::iterative::LSCG, jacobian_t>;
  lin_solver_t linearSolverObj;
  // nonlinear solvers
  auto nonLinSolver = pnonls::create_newton_raphson(problemObj, y, linearSolverObj);
  nonLinSolver.solve(problemObj, y);

  // check solution
  std::cout << "Computed solution: ["
            << y(0) << " " << y(1) << " " << "] "
            << "Expected solution: [1., 0.] "
            << std::endl;

  plog::finalize();
  return 0;
}

Full Code

The full code is available TODO.