Skip to content
Snippets Groups Projects
Commit 680d5ce2 authored by patavirt's avatar patavirt
Browse files

common: fix binops and Eigen binops for CppAD

parent fe5d2f90
No related branches found
No related tags found
No related merge requests found
......@@ -10,23 +10,66 @@
#include <iostream>
#include <Eigen/Core>
#include <cppad/cppad.hpp>
#include <cppad/example/cppad_eigen.hpp>
#include <Eigen/Core>
#include "array.hpp"
typedef std::complex<double> Complex;
typedef CppAD::AD<double> ADDouble;
typedef std::complex<ADDouble> ADComplex;
ADComplex adcomplex(const ADComplex& cpx)
template <typename Scalar_, int Rows_, int Cols_>
using Matrix = Eigen::Matrix<Scalar_, Rows_, Cols_, Eigen::RowMajor>;
ADComplex adcomplex(const ADComplex&& cpx)
{
return {cpx};
}
ADComplex adcomplex(const std::complex<double>& cpx)
ADComplex adcomplex(const std::complex<double>&& cpx)
{
return {cpx};
}
namespace Eigen
{
/* Eigen binop maps: (double, ADDouble) -> ADDouble */
template<typename BinaryOp> struct ScalarBinaryOpTraits<ADDouble,double,BinaryOp> { typedef ADDouble ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<double,ADDouble,BinaryOp> { typedef ADDouble ReturnType; };
/* Eigen binop maps: (double/complex, ADComplex) -> ADComplex */
template<typename BinaryOp> struct ScalarBinaryOpTraits<ADComplex,double,BinaryOp> { typedef ADComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<double,ADComplex,BinaryOp> { typedef ADComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<ADComplex,Complex,BinaryOp> { typedef ADComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<Complex,ADComplex,BinaryOp> { typedef ADComplex ReturnType; };
namespace usadelndsoc {
ADComplex operator*(const ADComplex& a, const double& b) { return a * adcomplex(b); }
ADComplex operator*(const double& a, const ADComplex& b) { return adcomplex(a) * b; }
ADComplex operator*(const ADComplex& a, const Complex& b) { return a * adcomplex(b); }
ADComplex operator*(const Complex& a, const ADComplex& b) { return adcomplex(a) * b; }
ADComplex operator+(const ADComplex& a, const double& b) { return a + adcomplex(b); }
ADComplex operator+(const double& a, const ADComplex& b) { return adcomplex(a) + b; }
ADComplex operator+(const ADComplex& a, const Complex& b) { return a + adcomplex(b); }
ADComplex operator+(const Complex& a, const ADComplex& b) { return adcomplex(a) + b; }
ADComplex operator-(const ADComplex& a, const double& b) { return a - adcomplex(b); }
ADComplex operator-(const double& a, const ADComplex& b) { return adcomplex(a) - b; }
ADComplex operator-(const ADComplex& a, const Complex& b) { return a - adcomplex(b); }
ADComplex operator-(const Complex& a, const ADComplex& b) { return adcomplex(a) - b; }
ADComplex operator/(const ADComplex& a, const double& b) { return a / adcomplex(b); }
ADComplex operator/(const double& a, const ADComplex& b) { return adcomplex(a) / b; }
ADComplex operator/(const ADComplex& a, const Complex& b) { return a / adcomplex(b); }
ADComplex operator/(const Complex& a, const ADComplex& b) { return adcomplex(a) / b; }
}
using namespace usadelndsoc;
}
using namespace Eigen::usadelndsoc;
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment