Newton-Raphson with Custom Types

Custom data types

template<class ScalarType>
struct CustomVector
{
  using value_type = ScalarType;
  using size_type = std::size_t;

  CustomVector(std::size_t ext) : d_(ext){}

  ScalarType & operator()(int i){ return d_[i]; }
  const ScalarType & operator()(int i)const { return d_[i]; }
  ScalarType & operator[](int i){ return d_[i]; }
  const ScalarType & operator[](int i)const { return d_[i]; }

  std::size_t extent(int k)const { return (k==0) ? d_.size() : 0; }

  void fill(ScalarType value){
    std::for_each(d_.begin(), d_.end(), [](ScalarType & v){ v= 0.; });
  }

private:
  std::vector<ScalarType> d_ = {};
};

template<class ScalarType>
struct CustomMatrix
{
  using value_type = ScalarType;
  using size_type = std::size_t;

  CustomMatrix(std::size_t nr, std::size_t nc)
    : num_rows_(nr), num_cols_(nc), d_(nr*nc){}

  std::size_t extent(int k)const { return (k==0) ? num_rows_ : num_cols_; }

  ScalarType & operator()(int i, int j){ return d_[num_cols_*i+j]; }
  const ScalarType & operator()(int i, int j) const { return d_[num_cols_*i+j]; }

  void fill(ScalarType value){
    std::for_each(d_.begin(), d_.end(), [=](ScalarType & v){ v= value; });
  }

private:
  std::size_t num_rows_ = {};
  std::size_t num_cols_ = {};
  std::vector<ScalarType> d_ = {};
};

Problem class

struct MySystem
{
  using scalar_type = double;
  using state_type = CustomVector<scalar_type>;
  using residual_type = state_type;
  using jacobian_type = CustomMatrix<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& J) const
  {
    J(0,0) = 3.0*x[0]*x[0];
    J(0,1) =  1.0;
    J(1,0) = -1.0;
    J(1,1) = 3.0*x[1]*x[1];
  }
};

Specialize trait and ops

Because we are working with custom data types, we need to provide the necessary operations to do the algebra that pressio needs. This is done via specialization as follows:

namespace pressio{
template<class T> struct Traits<CustomVector<T>>{
  using scalar_type = double;
  using size_type = std::size_t;
};

template<class T> struct Traits<CustomMatrix<T>>{
  using scalar_type = double;
  using size_type = std::size_t;
};

namespace ops{
template<class T> 
CustomVector<T> clone(const CustomVector<T> & src){ return src; }

template<class T>
CustomMatrix<T> clone(const CustomMatrix<T> & src){ return src; }

template<class T>
void set_zero(CustomVector<T> & o){ o.fill(0); }

template<class T>
void set_zero(CustomMatrix<T> & o){ o.fill(0); }

template<class T>
T norm2(const CustomVector<T> & v){
  return std::sqrt(v[0]*v[0] + v[1]*v[1]);
}

template<class T>
void update(CustomVector<T> & v, double a,
        const CustomVector<T> & v1, double b){
  v[0] = v[0]*a + b*v1[0];
  v[1] = v[1]*a + b*v1[1];
}

template<class T>
void scale(CustomVector<T> & v, double factor){
  v[0] = v[0]*factor;
  v[1] = v[1]*factor;
}
}}//end namespace pressio::ops

Custom Linear solver

In this example, we use a custom linear solver.

template<class ScalarType>
struct LinearSolver{
  using matrix_type = CustomMatrix<ScalarType>;

  void solve(const matrix_type & M, 
             const CustomVector<ScalarType> & rhs, 
             CustomVector<ScalarType> & x)
  {
    const auto a = M(0,0);
    const auto b = M(0,1);
    const auto c = M(1,0);
    const auto d = M(1,1);
    const auto det = a*d - b*c;

    x[0] = (d*rhs[0]  - b*rhs[1])/det;
    x[1] = (-c*rhs[0] + a*rhs[1])/det;
  }
};

Main

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

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

  using problem_t  = MySystem;
  using state_t    = problem_t::state_type;
  using scalar_t   = problem_t::scalar_type;
  problem_t sys;
  state_t y(2);
  y[0] = 0.001;
  y[1] = -0.1;

  LinearSolver<scalar_t> linearSolverObj;
  auto nonLinSolver = pnonls::create_newton_raphson(sys, y, linearSolverObj);
  nonLinSolver.solve(sys, 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.