Пользовательский скалярный тип в Eigen

Question

В настоящее время я пытаюсь настроить пользовательский скалярный тип для использования с библиотекой Eigen3 (атм, это простая оболочка для double). Я следовал https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html, насколько мне известно, и основы c работают отлично.

Мне нужно решить задачи на собственные значения для матриц с помощью моего пользовательского типа, и именно здесь все начинает разваливаться. Мой компилятор выдает мне следующее сообщение об ошибке:

/Eigen3/Eigen/src/Householder/Householder.h:131:18: error: no viable overloaded '-='
    this->row(0) -= tau * tmp;
    ~~~~~~~~~~~~ ^  ~~~~~~~~~
/Eigen3/Eigen/src/Eigenvalues/HessenbergDecomposition.h:314:10: note: in instantiation of function template specialization 'Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<MyDouble, 2,
      2, 0, 2, 2>, -1, -1, false> >::applyHouseholderOnTheLeft<Eigen::VectorBlock<Eigen::Block<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2>, 2, 1, true>, -1> >' requested here
        .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
         ^
/Eigen3/Eigen/src/Eigenvalues/HessenbergDecomposition.h:161:7: note: in instantiation of member function 'Eigen::HessenbergDecomposition<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2>
      >::_compute' requested here
      _compute(m_matrix, m_hCoeffs, m_temp);
      ^
/Eigen3/Eigen/src/Eigenvalues/./RealSchur.h:271:10: note: in instantiation of function template specialization 'Eigen::HessenbergDecomposition<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2>
      >::compute<Eigen::CwiseBinaryOp<Eigen::internal::scalar_quotient_op<MyDouble, MyDouble>, const Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<MyDouble>, const
      Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2> > > >' requested here
  m_hess.compute(matrix.derived()/scale);
         ^
/Eigen3/Eigen/src/Eigenvalues/EigenSolver.h:389:15: note: in instantiation of function template specialization 'Eigen::RealSchur<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2>
      >::compute<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2> >' requested here
  m_realSchur.compute(matrix.derived(), computeEigenvectors);
              ^
/Eigen3/Eigen/src/Eigenvalues/EigenSolver.h:156:7: note: in instantiation of function template specialization 'Eigen::EigenSolver<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2>
      >::compute<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2> >' requested here
      compute(matrix.derived(), computeEigenvectors);
      ^
/home/adam/Documents/Git/spin_phonon_coupling/src/test.cpp:205:52: note: in instantiation of function template specialization 'Eigen::EigenSolver<Eigen::Matrix<MyDouble, 2, 2, 0, 2, 2> >::EigenSolver<Eigen::Matrix<MyDouble, 2, 2, 0, 2,
      2> >' requested here
        Eigen::EigenSolver<Eigen::Matrix<MyDouble, 2, 2>> solver(test);
                                                          ^
/Eigen3/Eigen/src/Core/DenseBase.h:298:14: note: candidate template ignored: could not match 'EigenBase<type-parameter-0-0>' against 'MyDouble'
    Derived& operator-=(const EigenBase<OtherDerived> &other);
             ^
/Eigen3/Eigen/src/Core/MatrixBase.h:161:14: note: candidate template ignored: could not match 'MatrixBase<type-parameter-0-0>' against 'MyDouble'
    Derived& operator-=(const MatrixBase<OtherDerived>& other);
             ^
/Eigen3/Eigen/src/Core/MatrixBase.h:495:46: note: candidate template ignored: could not match 'ArrayBase<type-parameter-0-0>' against 'MyDouble'
    template<typename OtherDerived> Derived& operator-=(const ArrayBase<OtherDerived>& )

Таким образом, основной проблемой является отсутствующая перегрузка для оператора -=. Проблема в том, что я понятия не имею, почему. Насколько я вижу, я определил перегруженный оператор для своего типа, но если я правильно понимаю ошибку, мне кажется, что внутренний перегруженный тип пропускает эту перегрузку ...

У кого-нибудь есть идея, в чем может быть проблема здесь?

---

Вот код, выдающий ошибку выше:

#include <iostream>
#include <cmath>

#include <Eigen/Dense>

class MyDouble {
    public:
        double value;
        MyDouble() : value() {};
        MyDouble(double val) : value(val) {};

        template<typename T>
            MyDouble operator+(T other) const {
                return value + other;
            }
        template<>
            MyDouble operator+(MyDouble other) const {
                return MyDouble(value + other.value);
            }

        template<typename T>
            MyDouble operator-(T other) const {
                return value - other;
            }
        template<>
            MyDouble operator-(MyDouble other) const {
                return MyDouble(value - other.value);
            }

        template<typename T>
            MyDouble operator*(T other) const {
                return value * other;
            }
        template<>
            MyDouble operator*(MyDouble other) const {
                return MyDouble(value * other.value);
            }

        template<typename T>
            MyDouble operator/(T other) const {
                return value / other;
            }
        template<>
            MyDouble operator/(MyDouble other) const {
                return MyDouble(value / other.value);
            }

        template<typename T>
            MyDouble& operator+=(T other) {
                value += other;
                return *this;
            }
        template<>
            MyDouble& operator+=(MyDouble other) {
                value += other.value;
                return *this;
            }

        template<typename T>
            MyDouble& operator-=(const T &other) {
                value -= other;
                return *this;
            }
        template<>
            MyDouble& operator-=(const MyDouble &other) {
                value -= other.value;
                return *this;
            }

        template<typename T>
            MyDouble& operator*=(T other) {
                value *= other.value;
                return *this;
            }
        template<>
            MyDouble& operator*=(MyDouble other) {
                value *= other.value;
                return *this;
            }

        template<typename T>
            MyDouble& operator/=(T other) {
                value /= other;
                return *this;
            }
        template<>
            MyDouble& operator/=(MyDouble other) {
                value /= other.value;
                return *this;
            }

        MyDouble operator-() const {
            return -value;
        }

        template<typename T>
            bool operator<(T other) const {
                return value < other;
            }
        template<>
            bool operator<(MyDouble other) const {
                return value < other.value;
            }

        template<typename T>
            bool operator>(T other) const {
                return value > other;
            }
        template<>
            bool operator>(MyDouble other) const {
                return value > other.value;
            }

        template<typename T>
            bool operator<=(T other) const {
                return value <= other;
            }
        template<>
            bool operator <=(MyDouble other) const {
                return value <= other.value;
            }

        template<typename T>
            bool operator>=(T other) const {
                return value >= other;
            }
        template<>
            bool operator>=(MyDouble other) const {
                return value >= other.value;
            }

        template<typename T>
            bool operator==(T other) const {
                return value == other;
            }
        template<>
            bool operator==(MyDouble other) const {
                return value == other.value;
            }

        template<typename T>
            bool operator!=(T other) const {
                return value != other;
            }
        template<>
            bool operator!=(MyDouble other) const {
                return value != other.value;
            }

        friend std::ostream& operator<<(std::ostream& out, const MyDouble& val) {
            out << val.value << " m";
            return out;
        }

        operator double() {
            return value;
        }
};

MyDouble sqrt(MyDouble val) {
    return std::sqrt(val.value);
}
MyDouble abs(MyDouble val) {
    return std::abs(val.value);
}
MyDouble abs2(MyDouble val) {
    return val * val;
}

namespace Eigen {
    template<> struct NumTraits<MyDouble>
        : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
        {
            typedef MyDouble Real;
            typedef MyDouble NonInteger;
            typedef MyDouble Nested;
            enum {
                IsComplex = 0,
                IsInteger = 0,
                IsSigned = 1,
                RequireInitialization = 1,
                ReadCost = 1,
                AddCost = 3,
                MulCost = 3
            };
        };
    template<typename BinaryOp>
    struct ScalarBinaryOpTraits<MyDouble,double,BinaryOp> { typedef MyDouble ReturnType;  };

    template<typename BinaryOp>
    struct ScalarBinaryOpTraits<double,MyDouble,BinaryOp> { typedef MyDouble ReturnType;  };
}

int main() {
    Eigen::Matrix<MyDouble, 2, 2> test;
    test << 1, 2, 3, 4;

    Eigen::Matrix<double, 2, 2> test2;
    test2 << 5, 6, 7, 8;

    MyDouble a = 3;
    a -= 2;
    a -= MyDouble(3);

    Eigen::EigenSolver<Eigen::Matrix<MyDouble, 2, 2>> solver(test);
    std::cout << test.trace() << std::endl;
    std::cout << solver.eigenvalues() << std::endl;
}

Answer 1

Март c Глисс указал, что мои перегрузки операторов выглядят довольно странно, поэтому я переписал их, и первоначальная проблема исчезла. Я понятия не имею, почему именно это так.

Тогда у меня была проблема только в том, что функция Eigen isfinite не была определена для моего пользовательского типа, поэтому я продолжил и добавил ее реализацию (хотя я не уверен, что это на самом деле solid один).

В любом случае мой код компилируется сейчас. Вот модифицированная версия моего кода:

#include <iostream>
#include <cmath>
#include <complex>

#include <Eigen/Dense>

class MyDouble {
    public:
        double value;
        MyDouble() : value() {};
        MyDouble(double val) : value(val) {};

        template<typename T>
            MyDouble& operator+=(T rhs) {
                value = static_cast<double>(value + rhs);
                return *this;
            }

        template<typename T>
            MyDouble& operator-=(const T &rhs) {
                value = static_cast<double>(value - rhs);
                return *this;
            }

        template<typename T>
            MyDouble& operator*=(T rhs) {
                value = static_cast<double>(value * rhs);
                return *this;
            }

        template<typename T>
            MyDouble& operator/=(T rhs) {
                value = static_cast<double>(value / rhs);
                return *this;
            }

        MyDouble operator-() const {
            return -value;
        }

        friend std::ostream& operator<<(std::ostream& out, const MyDouble& val) {
            out << val.value << " m";
            return out;
        }

        explicit operator double() {
            return value;
        }
};

#define OVERLOAD_OPERATOR(op,ret) ret operator op(const MyDouble &lhs, const MyDouble &rhs) { \
        return lhs.value op rhs.value; \
    }

OVERLOAD_OPERATOR(+, MyDouble)
OVERLOAD_OPERATOR(-, MyDouble)
OVERLOAD_OPERATOR(*, MyDouble)
OVERLOAD_OPERATOR(/, MyDouble)

OVERLOAD_OPERATOR(>, bool)
OVERLOAD_OPERATOR(<, bool)
OVERLOAD_OPERATOR(>=, bool)
OVERLOAD_OPERATOR(<=, bool)
OVERLOAD_OPERATOR(==, bool)
OVERLOAD_OPERATOR(!=, bool)


MyDouble sqrt(MyDouble val) {
    return std::sqrt(val.value);
}
MyDouble abs(MyDouble val) {
    return std::abs(val.value);
}
MyDouble abs2(MyDouble val) {
    return val * val;
}
bool isfinite(const MyDouble &) { return true; }

namespace Eigen {
    template<> struct NumTraits<MyDouble>
        : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
        {
            typedef MyDouble Real;
            typedef MyDouble NonInteger;
            typedef MyDouble Nested;
            enum {
                IsComplex = 0,
                IsInteger = 0,
                IsSigned = 1,
                RequireInitialization = 0,
                ReadCost = 1,
                AddCost = 3,
                MulCost = 3
            };
        };

    template<typename BinaryOp>
    struct ScalarBinaryOpTraits<MyDouble,double,BinaryOp> { typedef MyDouble ReturnType;  };

    template<typename BinaryOp>
    struct ScalarBinaryOpTraits<double,MyDouble,BinaryOp> { typedef MyDouble ReturnType;  };
}

int main() {
    Eigen::Matrix<MyDouble, 2, 2> test;
    test << 1, 2, 3, 4;

    Eigen::Matrix<double, 2, 2> reference;
    reference << 1, 2, 3, 4;

    MyDouble a = 3;
    a += 2;
    a = 2 + a;
    a = a + 2;
    a -= 2;
    a -= MyDouble(3);

    a = a / a;

    std::complex<MyDouble> complexTest(3,4);
    complexTest *= 2;

    Eigen::EigenSolver<Eigen::Matrix<MyDouble, 2, 2>> solver(test);
    Eigen::EigenSolver<Eigen::Matrix<double, 2, 2>> refSolver(reference);
    std::cout << "MyDouble:" << std::endl;
    std::cout << test.trace() << std::endl;
    std::cout << solver.eigenvalues() << std::endl;
    std::cout << "\nRefernce:" << std::endl;
    std::cout << reference.trace() << std::endl;
    std::cout << refSolver.eigenvalues() << std::endl;
}