wip matrix free galerkin

include/pressio/rom/impl/galerkin_steady_system_default.hpp

@@ -65,6 +65,29 @@ class GalerkinSteadyDefaultSystem
     return impl::CreateGalerkinJacobian<jacobian_type>()(trialSubspace_.get().dimension());
   }
 
+  template<class _FomSystemType = FomSystemType>
+  std::enable_if_t<
+    std::is_same_v<
+      void,
+      decltype(
+	       std::declval<_FomSystemType const>().residual(
+							     std::declval<state_type const>(),
+							     std::declval<residual_type &>()
+							     )
+	       )
+    >
+  >
+  residual(const state_type & reducedState,
+	   residual_type & reducedResidual) const
+  {
+    const auto & phi = trialSubspace_.get().basisOfTranslatedSpace();
+    trialSubspace_.get().mapFromReducedState(reducedState, fomState_);
+
+    fomSystem_.get().residual(fomState_, fomResidual_);
+    ::pressio::ops::product(::pressio::transpose(),
+			    1, phi, fomResidual_, 0, reducedResidual);
+  }
+
   void residualAndJacobian(const state_type & reducedState,
 			   residual_type & reducedResidual,
 			   std::optional<jacobian_type*> reducedJacobian) const
@@ -94,6 +117,41 @@ class GalerkinSteadyDefaultSystem
     }
   }
 
+  template<class OperandT, class ResultT, class _FomSystemType = FomSystemType>
+  std::enable_if_t<
+    std::is_same_v<
+      void,
+      decltype(
+	       std::declval<_FomSystemType const>().jacobianAction(
+							     std::declval<state_type const &>(),
+							     std::declval<state_type const &>(),
+							     std::declval<state_type &>()
+							     )
+	       )
+      >
+  >
+  applyJacobian(const state_type & reducedState,
+		OperandT const & reducedOperand,
+		ResultT & out) const
+  {
+    trialSubspace_.get().mapFromReducedState(reducedState, fomState_);
+
+    auto auxVec = pressio::ops::clone(fomState_);
+    auto auxVec2 = pressio::ops::clone(fomState_);
+    pressio::ops::set_zero(auxVec);
+    trialSubspace_.get().mapFromReducedState(reducedOperand, auxVec2);
+
+    fomSystem_.get().jacobianAction(fomState_, reducedOperand, auxVec);
+
+    using scalar_t = typename ::pressio::Traits<basis_matrix_type>::scalar_type;
+    scalar_t alpha{1};
+    scalar_t beta{0};
+
+    const auto & phi = trialSubspace_.get().basisOfTranslatedSpace();
+    ::pressio::ops::product(::pressio::transpose(),
+			    alpha, phi, auxVec2, beta, out);
+  }
+
 private:
   std::reference_wrapper<const TrialSubspaceType> trialSubspace_;
   std::reference_wrapper<const FomSystemType> fomSystem_;

tests/functional_small/rom/CMakeLists.txt

@@ -29,6 +29,7 @@ endif()
 if(PRESSIO_ENABLE_TPL_EIGEN)
   set(SOURCES_GALERKIN_STEADY
     ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main1.cc
+    ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main2.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main3.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main4.cc)
   add_serial_utest(${TESTING_LEVEL}_rom_galerkin_steady ${SOURCES_GALERKIN_STEADY})

tests/functional_small/rom/galerkin_steady/main2.cc

@@ -0,0 +1,74 @@
+
+#include <gtest/gtest.h>
+#include "pressio/rom_subspaces.hpp"
+#include "pressio/rom_galerkin_steady.hpp"
+
+namespace{
+
+struct MyFom
+{
+  using state_type        = Eigen::VectorXd;
+  using residual_type     = state_type;
+  int N_ = {};
+
+  MyFom(int N): N_(N){}
+  residual_type createResidual() const{ return residual_type(N_); }
+
+  Eigen::MatrixXd createResultOfJacobianActionOn(const Eigen::MatrixXd & B) const{
+    Eigen::MatrixXd A(N_, B.cols());
+    return A;
+  }
+
+  void residual(const state_type & u,
+		residual_type & r) const
+  {}
+
+  template<class OperandType, class ResultType>
+  void jacobianAction(const state_type &,
+		      OperandType const &,
+		      ResultType &) const
+  {}
+
+  void residualAndJacobianAction(const state_type & u,
+				 residual_type & r,
+				 const Eigen::MatrixXd & B,
+				 std::optional<Eigen::MatrixXd *> Ain) const
+  {
+  }
+
+};
+
+}
+
+TEST(rom_galerkin_steady, default_matrix_free)
+{
+  using namespace pressio;
+
+  /*
+    just supposed to compile for now
+  */
+
+  log::initialize(logto::terminal);
+  log::setVerbosity({log::level::debug});
+
+  constexpr int N = 8;
+  using fom_t = MyFom;
+  fom_t fomSystem(N);
+
+  using phi_t = Eigen::Matrix<double, -1,-1>;
+  phi_t phi(N, 3);
+
+  using reduced_state_type = Eigen::VectorXd;
+  typename fom_t::state_type shift(N);
+  auto space = rom::create_trial_column_subspace<reduced_state_type>(phi, shift, false);
+
+  auto problem = rom::galerkin::create_steady_problem(space, fomSystem);
+
+  using tag = linearsolvers::iterative::GMRES;
+  auto nonLinSolver = experimental::create_matrixfree_newtonkrylov_solver<tag>(problem);
+
+  auto romState = space.createReducedState();
+  nonLinSolver.solve(romState);
+
+  log::finalize();
+}