nonlinear solvers: tutorial 1

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

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

  state_type createState() const {
    state_type a(2);
    a.setZero();
    return a;
  }

  residual_type createResidual() const {
    residual_type a(2);
    a.setZero();
    return a;
  }

  jacobian_type createJacobian() const {
    jacobian_type a(2, 2);
    a.setZero();
    return a;
  }

  void residualAndJacobian(const state_type& x,
			   residual_type& res,
			   std::optional<jacobian_type*> Jin) const
  {
    res(0) = x(0)*x(0) + x(1)*x(1) - 4.0;
    res(1) = x(0)*x(1) - 1.0;

    if (Jin){
     auto & jac = *Jin.value();
      jac.coeffRef(0, 0) = 2*x(0);
      jac.coeffRef(0, 1) = 2*x(1);
      // Have incorrect entries so that line search is required
      jac.coeffRef(1, 0) = 0.1*x(1);
      jac.coeffRef(1, 1) = 0.1*x(0);
    }
  }
};

class MyLinearSolver{
  using A_t = typename MyProblem::jacobian_type;
  using b_t = typename MyProblem::residual_type;
  using x_t = typename MyProblem::state_type;

public:
  void solve(A_t const & A, b_t const & b, x_t & x)
  {
    Eigen::BiCGSTAB< A_t > solver;
    solver.compute(A);
    x = solver.solve(b);
  }
};

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

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

  MyProblem problem;
  auto solver = pressio::create_newton_solver(problem, MyLinearSolver{});
  const auto updateMethod = pnonls::Update::BacktrackStrictlyDecreasingObjective;
  solver.setUpdateCriterion(updateMethod);

  auto x = problem.createState();
  x(0) = 0.3; x(1) = 0.4;
  solver.solve(x);

  // check solution
  std::cout << "Computed solution: \n "
	    << "["
	    << std::setprecision(14)
	    << x(0) << " "
	    << x(1) << " " << "] \n"
	    << "Expected solution: \n "
	    << "[1.93185165, 0.51763809] \n";

  plog::finalize();
  return 0;
}