Eigen Wrapper working

mfds/CMakeLists.txt

@@ -21,7 +21,7 @@ TARGET_LINK_LIBRARIES(mfds PUBLIC ${MFEM_LIBRARIES} hilbert_transformation_lib v
 target_include_directories(mfds
         PUBLIC ${MFEM_INCLUDE_DIRS}
         INTERFACE ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..)
-target_compile_features(mfds PRIVATE cxx_std_17)
+target_compile_features(mfds PRIVATE cxx_std_20)
 
 if(ENABLE_BASIC_TESTS)
     add_subdirectory(tests)

mfds/active_set_operator.cpp

@@ -115,7 +115,8 @@ namespace mfds {
         return space_time_boundary_dofs;
     }
 
-    Eigen::VectorXd heat_analytic_linear_operator_t::operator*(const Eigen::VectorXd &v) const {
+    Eigen::VectorXd heat_analytic_linear_operator_t::mult_internal(const Eigen::VectorXd &v) const {
+        std::cout << "heat: mult internal " << std::endl << std::flush;
         auto const vb = vtint::split_to_block_vector(v, t_points.size()-1);
         auto const res = this->operator*(vb);
         return res.to_vector();

mfds/active_set_operator.hpp

@@ -9,236 +9,134 @@
 #include "analytic_eigenfunctions.hpp"
 #include "vtint.hpp"
 
+#include "Eigen/Core"
 #include "Eigen/Sparse"
 #include "Eigen/Dense"
 
 #include "conjugate_gradient.hpp"
 #include <optional>
 
-namespace mfds {
 
+
+namespace mfds {
     template<class Op_t>
     class boundary_condition_wrapper_t;
+    class heat_analytic_linear_operator_t;
+}
+using Eigen::SparseMatrix;
+
 
+namespace Eigen::internal {
+    using namespace mfds;
+    // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
+    template<>
+    struct traits<heat_analytic_linear_operator_t>
+            : public Eigen::internal::traits<Eigen::SparseMatrix<double> > {
+    };
+    // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
     template<class Op_t>
-    class active_set_lin_op_t {
-    public:
-        active_set_lin_op_t(Op_t *op_, std::vector<bool> const &is_active_) : op(op_), is_active(is_active_) {};
-
-        template<class Vec_t>
-        auto
-        operator*(Vec_t const &v) const {
-            Vec_t v_mod = v;
-            for (size_t i = 0; i < is_active.size(); ++i) {
-                if (is_active[i]) {
-                    v_mod[i] = 0.;
-                }
-            }
-            Eigen::VectorXd res = (*op) * v_mod;
-            for (size_t i = 0; i < is_active.size(); ++i) {
-                if (is_active[i]) {
-                    res[i] = v[i];
-                }
-            }
-            return res;
+    struct traits<boundary_condition_wrapper_t<Op_t>>
+            : public Eigen::internal::traits<Eigen::SparseMatrix<double> > {
+    };
 }
 
-        Eigen::VectorXd modify_rhs(Eigen::VectorXd const &rhs) const {
-            Eigen::VectorXd rhs_mod = rhs;
-            for (size_t i = 0; i < is_active.size(); ++i) {
-                if (is_active[i]) {
-                    rhs_mod[i] *= op->operator()(i, i); //do this to avoid exploding condition number
-                }
-            }
-            return rhs_mod;
-        }
 
-    private:
-        Op_t *op;
-        std::vector<bool> const &is_active;
-    };
+namespace mfds {
 
     template<class Op_t>
-    class active_set_solver_t {
+    class boundary_condition_wrapper_t : public Eigen::EigenBase<boundary_condition_wrapper_t<Op_t>> {
     public:
-        active_set_solver_t(Op_t *op_, Eigen::VectorXd const &rhs_, std::vector<bool> const &isBoundary_,
-                            Eigen::VectorXd const &lowerBound_, Eigen::VectorXd const &upperBound_)
-                : op(op_), rhs(rhs_), isBoundary(isBoundary_), lowerBound(lowerBound_), upperBound(upperBound_) {}
-
-        Eigen::VectorXd solve(double const cg_tol, int const cg_max_iter, double const alpha = 0.5,
-                              std::optional<Eigen::VectorXd> const &u0 = std::nullopt,
-                              std::optional<Eigen::VectorXd> const &lambda0 = std::nullopt) {
-            Eigen::VectorXd u = u0.value_or(0.5*(lowerBound + upperBound));
-            Eigen::VectorXd lambdas = lambda0.value_or(Eigen::VectorXd::Zero(rhs.size()));
-            std::vector<bool> is_active_lower(rhs.size(), false);
-            std::vector<bool> is_active_upper(rhs.size(), false);
-            std::vector<bool> is_active(rhs.size(), false);
-            int n_changed = rhs.size() + 20;
-            int iter = 0;
-            while(iter < 1000){
-                //first go through all the conditions to see, which values are active
-                auto const is_active_old = is_active;
-                for (size_t i = 0; i < rhs.size(); ++i) {
-                    //default all to false, except condition is met
-                    is_active_lower[i] = false;
-                    is_active_upper[i] = false;
-                    if (lambdas[i] + alpha * (lowerBound[i] - u[i]) > 0) {
-                        is_active_lower[i] = !isBoundary[i];
-                    } else if (lambdas[i] + alpha * (upperBound[i] - u[i]) < 0) {
-                        is_active_upper[i] = !isBoundary[i];
-                    }
-                    is_active[i] = is_active_lower[i] || is_active_upper[i];
-                }
-                std::cout << "Changed indices: ";
-                n_changed = [is_active, is_active_old]() {
-                    int n = 0;
-                    for (size_t i = 0; i < is_active.size(); ++i) {
-                        if (is_active[i] != is_active_old[i]) {
-                            n++;
-                            std::cout << i << " ";
-                        }
-                    }
-                    return n;
-                }();
-                //rhs adaption has to be done different, as we do not fix values
-                Eigen::VectorXd cur_rhs = rhs;
-                {
-                    Eigen::VectorXd correction = Eigen::VectorXd::Zero(rhs.size());
-                    for (size_t i = 0; i < rhs.size(); ++i) {
-                        if (is_active[i]) {
-                            if (is_active_lower[i]) {
-                                correction[i] = lowerBound[i];
-                            } else {
-                                correction[i] = upperBound[i];
-                            }
-                        } else {
-                            correction[i] = 0.;
+        // Required typedefs, constants, and methods to be usable as Eigen Type
+        typedef double Scalar;
+        typedef double RealScalar;
+        typedef int StorageIndex;
+        enum {
+            ColsAtCompileTime = Eigen::Dynamic,
+            MaxColsAtCompileTime = Eigen::Dynamic,
+            IsRowMajor = false
+        };
+
+        Eigen::Index rows() const { return Op.rows(); }
+
+        Eigen::Index cols() const { return Op.cols(); }
+
+        template<typename Rhs>
+        Eigen::Product<boundary_condition_wrapper_t, Rhs, Eigen::AliasFreeProduct>
+        operator*(const Eigen::MatrixBase<Rhs> &x) const {
+            return Eigen::Product<boundary_condition_wrapper_t, Rhs, Eigen::AliasFreeProduct>(*this,
+                                                                                                       x.derived());
         }
+
+    public:
+        boundary_condition_wrapper_t(Op_t& Op_, std::vector<int> const& boundary_dofs_) : Op(Op_), boundary_dofs(boundary_dofs_){
+            is_boundary = std::vector<bool>(Op.cols(), false);
+            for(auto const dof: boundary_dofs){
+                is_boundary[dof] = true;
             }
-                    cur_rhs -= (*op) * correction;
         }
 
-                //now we (ab)use the boundary value wrapper to enforce the active set conditions
-                active_set_lin_op_t<Op_t> active_set_wrapper(op, is_active);
-
-                //solve system with cg
-                Eigen::VectorXd u0 = Eigen::VectorXd::Zero(rhs.size());
-                //cur_rhs = active_set_wrapper.modify_rhs(cur_rhs);
-                auto const u_active = conjugate_gradient(active_set_wrapper, cur_rhs, u0, cg_tol, cg_max_iter);
-
-                auto const lambdas_old = lambdas;
-                auto const u_old = u;
-                //distribute the solution to u and lambdas
-                for (size_t i = 0; i < rhs.size(); ++i) {
-                    if (is_active[i]) {
-                        lambdas[i] = -u_active[i]; //as we have the negative lambdas in the vector for CG
-                        if (is_active_lower[i]) {
-                            u[i] = lowerBound[i];
-                        } else {
-                            u[i] = upperBound[i];
+        [[nodiscard]] Eigen::VectorXd modify_rhs(Eigen::VectorXd const& rhs, Eigen::VectorXd const& u) const {
+            if(rhs.size() != u.size()){
+                throw std::runtime_error("Sizes of rhs and u do not match");
             }
-                    } else {
-                        u[i] = u_active[i];
-                        lambdas[i] = 0.;
+            if(rhs.size() != Op.cols()){
+                throw std::runtime_error("Sizes of rhs and operator do not match");
             }
+            Eigen::VectorXd rhs_mod = rhs;
+            for(Eigen::Index i=0; i<is_boundary.size(); ++i){
+                if(is_boundary[i]){
+                    rhs_mod(i) = u(i);
                 }
-
-                std::cout << std::endl;
-                std::cout << "n_changed: " << n_changed << std::endl;
-                if( (n_changed == 0 && iter > 0) || iter > 100){
-                    if(iter>100){
-                        std::cout << "Max Newton iterations reached" << std::endl;
             }
-                    break;
+            return rhs_mod;
         }
 
-                //calculate difference to old solution
-                auto const diff_u = u - u_old;
-                auto const diff_lambdas = lambdas - lambdas_old;
-                std::cout << "Diff u: " << diff_u.norm() << ", Diff lambdas: " << diff_lambdas.norm() << std::endl;
-                iter++;
-                double const omega = 0.5;
-                u = u_old + omega*diff_u;
-                lambdas = lambdas_old + omega*diff_lambdas;
+        [[nodiscard]] Eigen::VectorXd mult_internal(Eigen::VectorXd const &v) const {
+            std::cout << "Bwrap mult" << std::endl << std::flush;
+            Eigen::VectorXd res = Op*v;
+            for(auto const& dof: boundary_dofs){
+                res(dof) = v(dof);
+            };
+            return res;
         }
 
-            return u;
-        }
 
     private:
-        Op_t *op;
-        std::vector<bool> isBoundary; //we do not apply the active set strategy on boundary dofs
-        Eigen::VectorXd rhs;
-        Eigen::VectorXd lowerBound;
-        Eigen::VectorXd upperBound;
-    };
+        Op_t & Op;
+        std::vector<bool> is_boundary;
+        std::vector<int> boundary_dofs;
 
-    //we assume an operator, that supports the operator*
-    template<class Op_t>
-    class boundary_condition_wrapper_t {
-    public:
-        boundary_condition_wrapper_t(Op_t *op_, std::vector<bool> const &isBoundary_, std::vector<double> const &bVals_)
-                : op(op_),
-                  isBoundary(isBoundary_), bVals(bVals_) {
-            assert(isBoundary.size() == bVals.size());
-        }
 
-        template<class Vec_t>
-        auto operator*(Vec_t const &v) const {
-            auto rv = (*op) * v;
-            for (size_t i = 0; i < isBoundary.size(); ++i) {
-                if (isBoundary[i]) {
-                    rv[i] = v[i];
-                }
-            }
-            return rv;
-        }
+    };
 
-        template<class Vec_t>
-        auto adapt_rhs(Vec_t const &rhs) const {
-            size_t const n_dofs = isBoundary.size();
-            Vec_t u0(n_dofs);
-            for (size_t i = 0; i < n_dofs; ++i) {
-                if (isBoundary[i]) {
-                    u0[i] = bVals[i];
-                } else {
-                    u0[i] = 0.;
 
-                }
-            }
-            Vec_t rhs_adapted = rhs - (*op) * u0;
-            //set rhs to boundary value if we are at a boundary dof
-            for (size_t i = 0; i < n_dofs; ++i) {
-                if (isBoundary[i]) {
-                    rhs_adapted[i] = bVals[i];
-                }
-            }
-            return rhs_adapted;
-        }
+    class heat_analytic_linear_operator_t : public Eigen::EigenBase<heat_analytic_linear_operator_t> {
+    public:
+        // Required typedefs, constants, and methods to be usable as Eigen Type
+        typedef double Scalar;
+        typedef double RealScalar;
+        typedef int StorageIndex;
+        enum {
+            ColsAtCompileTime = Eigen::Dynamic,
+            MaxColsAtCompileTime = Eigen::Dynamic,
+            IsRowMajor = false
+        };
 
-        double operator()(size_t i, size_t j) const {
-            if (isBoundary[i] && i == j) {
-                return 1.;
-            } else {
-                return (*op)(i, j);
-            }
-        }
+        Index rows() const { return mass_x.rows() * mass_t.rows(); }
 
-        Eigen::VectorXd const& get_diagonal_preconditioner() const{
-            return op->get_diagonal_preconditioner();
-        }
+        Index cols() const { return mass_x.cols() * mass_t.cols(); }
 
-    private:
-        Op_t *op; //just use pointer and leave responsibility to user
-        std::vector<bool> isBoundary;
-        std::vector<double> bVals;
-    };
+        template<typename Rhs>
+        Eigen::Product<mfds::heat_analytic_linear_operator_t, Rhs, Eigen::AliasFreeProduct>
+        operator*(const Eigen::MatrixBase<Rhs> &x) const {
+            return Eigen::Product<heat_analytic_linear_operator_t, Rhs, Eigen::AliasFreeProduct>(*this,
+                                                                                                       x.derived());
+        }
 
-    class heat_analytic_linear_operator_t {
     public:
         vtint::block_vector_t operator*(vtint::block_vector_t const &v) const;
 
-        Eigen::VectorXd operator*(Eigen::VectorXd const &v) const;
+        Eigen::VectorXd mult_internal(Eigen::VectorXd const &v) const;
 
         heat_analytic_linear_operator_t(mfem::FiniteElementSpace &X, int const nt, double const T, double rho_);
 
@@ -269,4 +167,40 @@ namespace mfds {
 
 } // mfds
 
+
+//we have to put the specialization of matrix free vector multiplications here
+// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
+
+namespace Eigen::internal {
+
+    using namespace mfds;
+    template<typename Rhs>
+    struct generic_product_impl<heat_analytic_linear_operator_t, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
+            : generic_product_impl_base<heat_analytic_linear_operator_t, Rhs, generic_product_impl<heat_analytic_linear_operator_t, Rhs> > {
+        typedef typename Product<heat_analytic_linear_operator_t, Rhs>::Scalar Scalar;
+
+        template<typename Dest>
+        static void scaleAndAddTo(Dest &dst, const heat_analytic_linear_operator_t &lhs, const Rhs &rhs,
+                                  const Scalar &alpha) {
+
+            dst.noalias() += alpha * lhs.mult_internal(rhs);
+        }
+    };
+
+    template<typename Rhs, typename Op_t>
+    struct generic_product_impl<boundary_condition_wrapper_t<Op_t>, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
+            : generic_product_impl_base<boundary_condition_wrapper_t<Op_t>, Rhs, generic_product_impl<boundary_condition_wrapper_t<Op_t>, Rhs> > {
+        typedef typename Product<boundary_condition_wrapper_t<Op_t>, Rhs>::Scalar Scalar;
+
+        template<typename Dest>
+        static void scaleAndAddTo(Dest &dst, const boundary_condition_wrapper_t<Op_t> &lhs, const Rhs &rhs,
+                                  const Scalar &alpha) {
+
+            dst.noalias() += alpha * lhs.mult_internal(rhs);
+        }
+    };
+
+}
+
+
 #endif //MFDS_ACTIVE_SET_OPERATOR_HPP

playground/example_matrix_free_operator_2d.cpp

@@ -11,6 +11,7 @@
 #include "conjugate_gradient.hpp"
 #include "homogenization.hpp"
 #include "time_dependent_gridfunction.hpp"
+#include <Eigen/IterativeLinearSolvers>
 
 static double target_function(mfem::Vector const &p, double const t) {
     double const x = p(0);
@@ -28,7 +29,7 @@ static double upper_bound_fun(double const t, mfem::Vector const &p) {
 }
 
 int main(int argc, char *argv[]) {
-    int const refinement = 3;
+    int const refinement = 6;
     int const nx = 2 * std::pow(2, refinement);
     int const ny = 2 * std::pow(2, refinement);
     int const nt = 2 * std::pow(2, refinement);
@@ -38,7 +39,7 @@ int main(int argc, char *argv[]) {
     double const newton_relax_param = 1;
 
 
-    mfem::Mesh spatial_mesh = mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::TRIANGLE, true, 1.0, 1.0);
+    mfem::Mesh spatial_mesh = mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::QUADRILATERAL, true, 1.0, 1.0);
     int const dimX = spatial_mesh.Dimension();
     mfem::Mesh temporal_mesh = mfem::Mesh::MakeCartesian1D(nt, T);
 
@@ -64,6 +65,23 @@ int main(int argc, char *argv[]) {
     //calculate inhomogeneous rhs for target
     auto const rhs_inhom = mfds::compute_spatially_inhom_rhs(target_function, T_test, X);
 
+    Eigen::VectorXd const u_bdr = Eigen::VectorXd::Zero(X.GetTrueVSize() * nt);
+
+    //create bdr condition wrapper
+    mfds::boundary_condition_wrapper_t boundary_condition_wrapper(heat_operator, space_time_boundary_dofs);
+    auto const rhs_bdr = boundary_condition_wrapper.modify_rhs(rhs_inhom.to_vector(), u_bdr);
+
+
+
+    Eigen::ConjugateGradient<mfds::boundary_condition_wrapper_t<mfds::heat_analytic_linear_operator_t>,
+            Eigen::Lower | Eigen::Upper, Eigen::IdentityPreconditioner> cg;
+    cg.compute(boundary_condition_wrapper);
+    Eigen::VectorXd const u_vec = cg.solve(rhs_bdr);
+    std::cout << "CG:       #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
+
+
+    /*
+
     //wrap operator so that boundary conditions can be applied
     //with zero boundary condition
     std::vector<double> const bVals(X.GetTrueVSize() * nt, 0.);
@@ -87,6 +105,9 @@ int main(int argc, char *argv[]) {
 
     auto const u_vec = active_set_solver.solve(cg_tol, cg_max_iter, newton_relax_param, u0, lambda0);
 
+
+     */
+
     //add initial condition
     auto const state_block = vtint::split_to_block_vector(u_vec, nt);
     mfem::Vector u0_vec(X.GetTrueVSize());