Operator.h

Go to the documentation of this file.

//
// Created by Nikhil Chandra Admal on 7/5/23.
//
 
#ifndef OILAB_OPERATOR_H
#define OILAB_OPERATOR_H
 
#include "../Lattices/LatticeModule.h"
#include <unsupported/Eigen/CXX11/Tensor>
 
namespace oILAB {
template <typename Derived, int dim> class Operator {
public:
  using Scalar = double;
  const Derived &derivedOperator;
  const Lattice<dim> L;
  const Eigen::array<Eigen::Index, dim> n;
 
  Operator(const Eigen::Matrix<double, dim, dim> &A,
           const Eigen::array<Eigen::Index, dim> &n_)
      : derivedOperator(static_cast<const Derived &>(*this)), L(A), n(n_) {}
 
  auto domain() const { return L.latticeBasis; }
 
  Eigen::Index rows() const {
    return std::accumulate(begin(n), end(n), 1, std::multiplies<>());
  }
  Eigen::Index cols() const {
    return std::accumulate(begin(n), end(n), 1, std::multiplies<>());
  }
 
  void perform_op(const double *x_in, double *y_out) const {
    derivedOperator.perform_op(x_in, y_out);
  }
 
  // implement expression templates for operators
 
    };
 
    template<typename E1, typename E2, int dim>
    class OperatorSum : public Operator<OperatorSum<E1,E2,dim>,dim> {
        const E1 o1;
        const E2 o2;
    public:
        OperatorSum(const E1 &o1_, const E2 &o2_) : Operator<OperatorSum<E1,E2,dim>,dim>(o1_.L.latticeBasis,o1_.n),
                                                    o1(o1_), o2(o2_)
 
        {
            // ensure that the domains of the two operators and their discretizations match
            assert(o1.n == o2.n && o1.domain().isApprox(o2.domain()));
        }
 
        void perform_op(const double *x_in, double *y_out) const {
            int nx = this->rows();
            Eigen::VectorXd y1(nx);
            Eigen::VectorXd y2(nx);
            o1.perform_op(x_in, y1.data());
            o2.perform_op(x_in, y2.data());
            Eigen::Map<Eigen::VectorXd> y(y_out, nx);
            y = y1 + y2;
        }
    };
 
    template<typename E1, typename E2, int dim>
    OperatorSum<E1, E2, dim> operator+(const Operator<E1, dim>& u, const Operator<E2,dim>& v) {
        return OperatorSum<E1, E2, dim>(*static_cast<const E1 *>(&u), *static_cast<const E2 *>(&v));
    }
 
    template<typename T, typename E, int dim>
    class OperatorScalarProduct : public Operator<OperatorScalarProduct<T,E,dim>,dim> {
        double s;
        const E op;
    public:
        OperatorScalarProduct(const T& s_, const E& op_) : Operator<OperatorScalarProduct<T,E,dim>,dim>(op_.L.latticeBasis,op_.n),
                                                           s(s_), op(op_)
 
        {}
 
        void perform_op(const double* x_in, double* y_out) const {
            int nx = this->rows();
            op.perform_op(x_in, y_out);
            Eigen::Map<Eigen::VectorXd> y(y_out, nx);
            y = s * y;
        }
    };
 
    template<typename T, typename E, int dim>
    OperatorScalarProduct<T, E, dim> operator*(const T &s, const Operator<E, dim> &u) {
        return OperatorScalarProduct<T, E, dim>(s, *static_cast<const E *>(&u));
    }
 
    } // namespace oILAB
#endif //OILAB_OPERATOR_H