devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Remove assumption of std::complex for complex scalar types.

Browse Source
This commit is contained in:
Antonio Sanchez
2025-02-12 11:21:44 -08:00
parent 6b4881ad48
commit 22cd7307dd
21 changed files with 273 additions and 115 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
View File

@@ -23,8 +23,10 @@ namespace internal {
*
* This struct is used by CwiseUnaryOp to scale a matrix by \f$ 2^{-s} \f$.
*/
template <typename RealScalar>
template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct MatrixExponentialScalingOp {
using RealScalar = typename NumTraits<Scalar>::Real;
/** \brief Constructor.
*
* \param[in] squarings The integer \f$ s \f$ in this document.
@@ -35,20 +37,30 @@ struct MatrixExponentialScalingOp {
*
* \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$.
*/
inline const RealScalar operator()(const RealScalar& x) const {
inline const Scalar operator()(const Scalar& x) const {
using std::ldexp;
return ldexp(x, -m_squarings);
return Scalar(ldexp(Eigen::numext::real(x), -m_squarings), ldexp(Eigen::numext::imag(x), -m_squarings));
}
typedef std::complex<RealScalar> ComplexScalar;
private:
int m_squarings;
};
template <typename Scalar>
struct MatrixExponentialScalingOp<Scalar, /*IsComplex=*/false> {
/** \brief Constructor.
*
* \param[in] squarings The integer \f$ s \f$ in this document.
*/
MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) {}
/** \brief Scale a matrix coefficient.
*
* \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$.
*/
inline const ComplexScalar operator()(const ComplexScalar& x) const {
inline const Scalar operator()(const Scalar& x) const {
using std::ldexp;
return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
return ldexp(x, -m_squarings);
}
private:
@@ -220,6 +232,7 @@ struct matrix_exp_computeUV {
template <typename MatrixType>
struct matrix_exp_computeUV<MatrixType, float> {
using Scalar = typename traits<MatrixType>::Scalar;
template <typename ArgType>
static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
using std::frexp;
@@ -234,7 +247,7 @@ struct matrix_exp_computeUV<MatrixType, float> {
const float maxnorm = 3.925724783138660f;
frexp(l1norm / maxnorm, &squarings);
if (squarings < 0) squarings = 0;
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings));
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
matrix_exp_pade7(A, U, V);
}
}
@@ -242,12 +255,12 @@ struct matrix_exp_computeUV<MatrixType, float> {
template <typename MatrixType>
struct matrix_exp_computeUV<MatrixType, double> {
typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
using Scalar = typename traits<MatrixType>::Scalar;
template <typename ArgType>
static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
using std::frexp;
using std::pow;
const RealScalar l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
squarings = 0;
if (l1norm < 1.495585217958292e-002) {
matrix_exp_pade3(arg, U, V);
@@ -258,10 +271,10 @@ struct matrix_exp_computeUV<MatrixType, double> {
} else if (l1norm < 2.097847961257068e+000) {
matrix_exp_pade9(arg, U, V);
} else {
const RealScalar maxnorm = 5.371920351148152;
const double maxnorm = 5.371920351148152;
frexp(l1norm / maxnorm, &squarings);
if (squarings < 0) squarings = 0;
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<RealScalar>(squarings));
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
matrix_exp_pade13(A, U, V);
}
}
@@ -271,6 +284,7 @@ template <typename MatrixType>
struct matrix_exp_computeUV<MatrixType, long double> {
template <typename ArgType>
static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
using Scalar = typename traits<MatrixType>::Scalar;
#if LDBL_MANT_DIG == 53 // double precision
matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
@@ -295,7 +309,7 @@ struct matrix_exp_computeUV<MatrixType, long double> {
const long double maxnorm = 4.0246098906697353063L;
frexp(l1norm / maxnorm, &squarings);
if (squarings < 0) squarings = 0;
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
matrix_exp_pade13(A, U, V);
}
@@ -315,7 +329,7 @@ struct matrix_exp_computeUV<MatrixType, long double> {
const long double maxnorm = 3.2579440895405400856599663723517L;
frexp(l1norm / maxnorm, &squarings);
if (squarings < 0) squarings = 0;
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
matrix_exp_pade17(A, U, V);
}
@@ -335,7 +349,7 @@ struct matrix_exp_computeUV<MatrixType, long double> {
const long double maxnorm = 2.884233277829519311757165057717815L;
frexp(l1norm / maxnorm, &squarings);
if (squarings < 0) squarings = 0;
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
matrix_exp_pade17(A, U, V);
}
@@ -382,9 +396,7 @@ template <typename ArgType, typename ResultType>
void matrix_exp_compute(const ArgType& arg, ResultType& result, false_type) // default
{
typedef typename ArgType::PlainObject MatrixType;
typedef typename traits<MatrixType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename std::complex<RealScalar> ComplexScalar;
typedef make_complex_t<typename traits<MatrixType>::Scalar> ComplexScalar;
result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
}

4
unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
View File

@@ -382,7 +382,7 @@ struct matrix_function_compute<MatrixType, 0> {
static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime;
static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime;
typedef std::complex<Scalar> ComplexScalar;
typedef internal::make_complex_t<Scalar> ComplexScalar;
typedef Matrix<ComplexScalar, Rows, Cols, 0, MaxRows, MaxCols> ComplexMatrix;
ComplexMatrix CA = A.template cast<ComplexScalar>();
@@ -476,7 +476,7 @@ class MatrixFunctionReturnValue : public ReturnByValue<MatrixFunctionReturnValue
typedef typename internal::nested_eval<Derived, 10>::type NestedEvalType;
typedef internal::remove_all_t<NestedEvalType> NestedEvalTypeClean;
typedef internal::traits<NestedEvalTypeClean> Traits;
typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
typedef internal::make_complex_t<Scalar> ComplexScalar;
typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, Traits::RowsAtCompileTime, Traits::ColsAtCompileTime>
DynMatrixType;

2
unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
View File

@@ -330,7 +330,7 @@ class MatrixLogarithmReturnValue : public ReturnByValue<MatrixLogarithmReturnVal
typedef typename internal::nested_eval<Derived, 10>::type DerivedEvalType;
typedef internal::remove_all_t<DerivedEvalType> DerivedEvalTypeClean;
typedef internal::traits<DerivedEvalTypeClean> Traits;
typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
typedef internal::make_complex_t<Scalar> ComplexScalar;
typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, Traits::RowsAtCompileTime, Traits::ColsAtCompileTime>
DynMatrixType;
typedef internal::MatrixLogarithmAtomic<DynMatrixType> AtomicType;

8
unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
View File

@@ -91,7 +91,7 @@ class MatrixPowerAtomic : internal::noncopyable {
enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime };
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef std::complex<RealScalar> ComplexScalar;
typedef internal::make_complex_t<Scalar> ComplexScalar;
typedef Block<MatrixType, Dynamic, Dynamic> ResultType;
const MatrixType& m_A;
@@ -380,7 +380,7 @@ class MatrixPower : internal::noncopyable {
Index cols() const { return m_A.cols(); }
private:
typedef std::complex<RealScalar> ComplexScalar;
typedef internal::make_complex_t<Scalar> ComplexScalar;
typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime>
ComplexMatrix;
@@ -628,7 +628,7 @@ template <typename Derived>
class MatrixComplexPowerReturnValue : public ReturnByValue<MatrixComplexPowerReturnValue<Derived> > {
public:
typedef typename Derived::PlainObject PlainObject;
typedef typename std::complex<typename Derived::RealScalar> ComplexScalar;
typedef internal::make_complex_t<typename Derived::Scalar> ComplexScalar;
/**
* \brief Constructor.
@@ -685,7 +685,7 @@ const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar&
}
template <typename Derived>
const MatrixComplexPowerReturnValue<Derived> MatrixBase<Derived>::pow(const std::complex<RealScalar>& p) const {
const MatrixComplexPowerReturnValue<Derived> MatrixBase<Derived>::pow(const internal::make_complex_t<Scalar>& p) const {
return MatrixComplexPowerReturnValue<Derived>(derived(), p);
}