test/complex_components.cpp - eigen - Git at Google

 // SPDX-FileCopyrightText: The Eigen Authors
 // SPDX-License-Identifier: MPL-2.0

 #include "main.h"
 #include <Eigen/Core>

 // Component access must not assume that a user-defined complex scalar has std::complex's layout.
 struct PaddedComplex : std::complex<double> {
   using std::complex<double>::complex;
   double padding = 123;
 };

 template <>
 struct Eigen::NumTraits<PaddedComplex> : Eigen::NumTraits<std::complex<double>> {};

 template <typename Derived>
 void check_component(const DenseBase<Derived>& component) {
   using Real = typename Derived::Scalar;
   Matrix<Real, Dynamic, Dynamic, Derived::IsRowMajor ? RowMajor : ColMajor> copy = component;
   Real sum = 0, minimum = component(0, 0), maximum = minimum, abs_maximum = numext::abs(minimum);
   for (Index col = 0; col < component.cols(); ++col) {
     for (Index row = 0; row < component.rows(); ++row) {
       const Real value = component(row, col);
       VERIFY_IS_EQUAL(copy(row, col), value);
       sum += value;
       if (value < minimum) minimum = value;
       if (value > maximum) maximum = value;
       if (numext::abs(value) > abs_maximum) abs_maximum = numext::abs(value);
     }
   }
   VERIFY_IS_EQUAL(component.sum(), sum);
   VERIFY_IS_EQUAL(component.template minCoeff<PropagateNaN>(), minimum);
   VERIFY_IS_EQUAL(component.template maxCoeff<PropagateNaN>(), maximum);
   VERIFY_IS_EQUAL(component.derived().cwiseAbs().template maxCoeff<PropagateNaN>(), abs_maximum);
   VERIFY_IS_CWISE_EQUAL(component.derived().colwise().sum().eval(), copy.colwise().sum().eval());
   VERIFY_IS_CWISE_EQUAL(component.derived().rowwise().sum().eval(), copy.rowwise().sum().eval());
 }

 template <typename MatrixType>
 void component_expressions(Index rows, Index cols) {
   using Complex = typename MatrixType::Scalar;
   using Real = typename NumTraits<Complex>::Real;
   MatrixType values(rows, cols);
   for (Index i = 0; i < values.size(); ++i) values.data()[i] = Complex(Real(i % 13 - 6), Real(i % 17 + 1));
   const MatrixType& const_values = values;
   check_component(values.real());
   check_component(values.imag());
   check_component(const_values.real());
   check_component(const_values.imag());
   check_component(values.transpose().real());
   check_component(values.transpose().imag());
   check_component(values.conjugate().real());
   check_component(values.conjugate().imag());
   check_component(values.block(0, 0, rows - 1, cols - 1).real());
   check_component(values.block(0, 0, rows - 1, cols - 1).imag());
   check_component(const_values.block(0, 0, rows - 1, cols - 1).real());
   check_component(const_values.block(0, 0, rows - 1, cols - 1).imag());
   const auto real = values.real();
   const auto imag = values.imag();
   check_component(real.block(0, 0, rows - 1, cols - 1));
   check_component(imag.block(0, 0, rows - 1, cols - 1));

   MatrixType original = values;
   values.real() = original.imag();
   values.imag() = original.real();
   for (Index i = 0; i < values.size(); ++i) {
     VERIFY_IS_EQUAL(values.data()[i], Complex(original.data()[i].imag(), original.data()[i].real()));
   }
   values.real().block(0, 0, rows - 1, cols - 1).setConstant(Real(3));
   values.imag().array() += Real(2);
   for (Index col = 0; col < cols; ++col) {
     for (Index row = 0; row < rows; ++row) {
       VERIFY_IS_EQUAL(values(row, col).imag(), original(row, col).real() + Real(2));
       VERIFY_IS_EQUAL(values(row, col).real(), row < rows - 1 && col < cols - 1 ? Real(3) : original(row, col).imag());
     }
   }
 }

 template <typename Real>
 void component_reductions() {
   using Complex = std::complex<Real>;
   using Vector = Matrix<Complex, Dynamic, 1>;
   constexpr bool vectorizable = (std::is_same<Real, float>::value || std::is_same<Real, double>::value) &&
                                 internal::packet_traits<Real>::Vectorizable;
   using View = decltype(std::declval<Vector&>().real());
   using ConstOp = decltype(std::declval<const Vector&>().imag());
   STATIC_CHECK(bool(internal::redux_evaluator<View>::Flags & PacketAccessBit) == vectorizable);
   STATIC_CHECK(bool(internal::evaluator<ConstOp>::Flags & PacketAccessBit) == vectorizable);
   STATIC_CHECK(!(internal::evaluator<View>::Flags & PacketAccessBit));
   using Strided = Map<Vector, Unaligned, InnerStride<2>>;
   using StridedView = decltype(std::declval<Strided&>().real());
   STATIC_CHECK(!(internal::redux_evaluator<StridedView>::Flags & PacketAccessBit));
   using DynamicStride = Map<Vector, Unaligned, InnerStride<Dynamic>>;
   using DynamicView = decltype(std::declval<const DynamicStride&>().imag());
   STATIC_CHECK(!(internal::evaluator<DynamicView>::Flags & PacketAccessBit));
   using SumView = decltype((std::declval<Vector>() + std::declval<Vector>()).real());
   STATIC_CHECK(!(internal::evaluator<SumView>::Flags & PacketAccessBit));

   const Real nan = NumTraits<Real>::quiet_NaN();
   const Real inf = NumTraits<Real>::infinity();
   constexpr Index packet_size = internal::packet_traits<Real>::size;
   for (Index size : {Index(1), packet_size, packet_size + 1, 2 * packet_size - 1, Index(65)}) {
     Vector buffer(2 * size + 1);
     buffer.setConstant(Complex(Real(7), Real(-3)));
     Map<Vector> values(buffer.data() + 1, size);
     const auto& const_values = values;
     check_component(values.real());
     check_component(values.imag());
     Strided strided(buffer.data() + 1, size);
     check_component(strided.real());
     check_component(strided.imag());
     DynamicStride dynamic(buffer.data() + 1, size, InnerStride<Dynamic>(2));
     check_component(dynamic.real());
     check_component(dynamic.imag());
     for (Index i = 0; i < size; ++i) {
       for (Real special : {nan, inf, -inf, Real(0), -Real(0), std::numeric_limits<Real>::denorm_min()}) {
         values.setConstant(Complex(Real(7), Real(-3)));
         values(i) = Complex(special, Real(-3));
         Matrix<Real, Dynamic, 1> real = const_values.real();
         if ((numext::isnan)(special)) {
           VERIFY((numext::isnan)(real(i)));
           VERIFY((numext::isnan)(values.real().template maxCoeff<PropagateNaN>()));
           VERIFY((numext::isnan)(const_values.real().template minCoeff<PropagateNaN>()));
           if (size > 1) VERIFY_IS_EQUAL(values.real().template maxCoeff<PropagateNumbers>(), Real(7));
         } else {
           VERIFY_IS_EQUAL(real(i), special);
           VERIFY_IS_EQUAL((std::signbit)(real(i)), (std::signbit)(special));
         }
         VERIFY_IS_EQUAL(values.imag().template maxCoeff<PropagateNaN>(), Real(-3));
         values(i) = Complex(Real(7), special);
         Matrix<Real, Dynamic, 1> imag = const_values.imag();
         if ((numext::isnan)(special)) {
           VERIFY((numext::isnan)(imag(i)));
           VERIFY((numext::isnan)(values.imag().template minCoeff<PropagateNaN>()));
           VERIFY((numext::isnan)(const_values.imag().template maxCoeff<PropagateNaN>()));
           if (size > 1) VERIFY_IS_EQUAL(values.imag().template minCoeff<PropagateNumbers>(), Real(-3));
         } else {
           VERIFY_IS_EQUAL(imag(i), special);
           VERIFY_IS_EQUAL((std::signbit)(imag(i)), (std::signbit)(special));
         }
         VERIFY_IS_EQUAL(values.real().template minCoeff<PropagateNaN>(), Real(7));
       }
     }
   }
 }

 void component_custom_scalar() {
   using Vector = Matrix<PaddedComplex, Dynamic, 1>;
   Vector values(17);
   for (Index i = 0; i < values.size(); ++i) values(i) = PaddedComplex(double(i), double(2 * i + 1));
   const Vector& const_values = values;
   using RealView = decltype(const_values.real());
   using ImagView = decltype(const_values.imag());
   STATIC_CHECK(!(internal::evaluator<RealView>::Flags & PacketAccessBit));
   STATIC_CHECK(!(internal::evaluator<ImagView>::Flags & PacketAccessBit));
   check_component(const_values.real());
   check_component(const_values.imag());
 }

 EIGEN_DECLARE_TEST(complex_components) {
   CALL_SUBTEST_1(component_reductions<float>());
   CALL_SUBTEST_2(component_reductions<double>());
   CALL_SUBTEST_3(component_reductions<long double>());
   CALL_SUBTEST_3(component_custom_scalar());
   for (Index rows : {Index(2), Index(7), Index(17)}) {
     for (Index cols : {Index(3), Index(9)}) {
       EIGEN_UNUSED_VARIABLE(rows);
       EIGEN_UNUSED_VARIABLE(cols);
       CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, Dynamic, Dynamic, ColMajor>>(rows, cols)));
       CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>(rows, cols)));
       CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, Dynamic, Dynamic, ColMajor>>(rows, cols)));
       CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>(rows, cols)));
     }
   }
   CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, 3, 5>>(3, 5)));
   CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, 2, 2>>(2, 2)));
 }
	// SPDX-FileCopyrightText: The Eigen Authors
	// SPDX-License-Identifier: MPL-2.0

	#include "main.h"
	#include <Eigen/Core>

	// Component access must not assume that a user-defined complex scalar has std::complex's layout.
	struct PaddedComplex : std::complex<double> {
	using std::complex<double>::complex;
	double padding = 123;
	};

	template <>
	struct Eigen::NumTraits<PaddedComplex> : Eigen::NumTraits<std::complex<double>> {};

	template <typename Derived>
	void check_component(const DenseBase<Derived>& component) {
	using Real = typename Derived::Scalar;
	Matrix<Real, Dynamic, Dynamic, Derived::IsRowMajor ? RowMajor : ColMajor> copy = component;
	Real sum = 0, minimum = component(0, 0), maximum = minimum, abs_maximum = numext::abs(minimum);
	for (Index col = 0; col < component.cols(); ++col) {
	for (Index row = 0; row < component.rows(); ++row) {
	const Real value = component(row, col);
	VERIFY_IS_EQUAL(copy(row, col), value);
	sum += value;
	if (value < minimum) minimum = value;
	if (value > maximum) maximum = value;
	if (numext::abs(value) > abs_maximum) abs_maximum = numext::abs(value);
	}
	}
	VERIFY_IS_EQUAL(component.sum(), sum);
	VERIFY_IS_EQUAL(component.template minCoeff<PropagateNaN>(), minimum);
	VERIFY_IS_EQUAL(component.template maxCoeff<PropagateNaN>(), maximum);
	VERIFY_IS_EQUAL(component.derived().cwiseAbs().template maxCoeff<PropagateNaN>(), abs_maximum);
	VERIFY_IS_CWISE_EQUAL(component.derived().colwise().sum().eval(), copy.colwise().sum().eval());
	VERIFY_IS_CWISE_EQUAL(component.derived().rowwise().sum().eval(), copy.rowwise().sum().eval());
	}

	template <typename MatrixType>
	void component_expressions(Index rows, Index cols) {
	using Complex = typename MatrixType::Scalar;
	using Real = typename NumTraits<Complex>::Real;
	MatrixType values(rows, cols);
	for (Index i = 0; i < values.size(); ++i) values.data()[i] = Complex(Real(i % 13 - 6), Real(i % 17 + 1));
	const MatrixType& const_values = values;
	check_component(values.real());
	check_component(values.imag());
	check_component(const_values.real());
	check_component(const_values.imag());
	check_component(values.transpose().real());
	check_component(values.transpose().imag());
	check_component(values.conjugate().real());
	check_component(values.conjugate().imag());
	check_component(values.block(0, 0, rows - 1, cols - 1).real());
	check_component(values.block(0, 0, rows - 1, cols - 1).imag());
	check_component(const_values.block(0, 0, rows - 1, cols - 1).real());
	check_component(const_values.block(0, 0, rows - 1, cols - 1).imag());
	const auto real = values.real();
	const auto imag = values.imag();
	check_component(real.block(0, 0, rows - 1, cols - 1));
	check_component(imag.block(0, 0, rows - 1, cols - 1));

	MatrixType original = values;
	values.real() = original.imag();
	values.imag() = original.real();
	for (Index i = 0; i < values.size(); ++i) {
	VERIFY_IS_EQUAL(values.data()[i], Complex(original.data()[i].imag(), original.data()[i].real()));
	}
	values.real().block(0, 0, rows - 1, cols - 1).setConstant(Real(3));
	values.imag().array() += Real(2);
	for (Index col = 0; col < cols; ++col) {
	for (Index row = 0; row < rows; ++row) {
	VERIFY_IS_EQUAL(values(row, col).imag(), original(row, col).real() + Real(2));
	VERIFY_IS_EQUAL(values(row, col).real(), row < rows - 1 && col < cols - 1 ? Real(3) : original(row, col).imag());
	}
	}
	}

	template <typename Real>
	void component_reductions() {
	using Complex = std::complex<Real>;
	using Vector = Matrix<Complex, Dynamic, 1>;
	constexpr bool vectorizable = (std::is_same<Real, float>::value \|\| std::is_same<Real, double>::value) &&
	internal::packet_traits<Real>::Vectorizable;
	using View = decltype(std::declval<Vector&>().real());
	using ConstOp = decltype(std::declval<const Vector&>().imag());
	STATIC_CHECK(bool(internal::redux_evaluator<View>::Flags & PacketAccessBit) == vectorizable);
	STATIC_CHECK(bool(internal::evaluator<ConstOp>::Flags & PacketAccessBit) == vectorizable);
	STATIC_CHECK(!(internal::evaluator<View>::Flags & PacketAccessBit));
	using Strided = Map<Vector, Unaligned, InnerStride<2>>;
	using StridedView = decltype(std::declval<Strided&>().real());
	STATIC_CHECK(!(internal::redux_evaluator<StridedView>::Flags & PacketAccessBit));
	using DynamicStride = Map<Vector, Unaligned, InnerStride<Dynamic>>;
	using DynamicView = decltype(std::declval<const DynamicStride&>().imag());
	STATIC_CHECK(!(internal::evaluator<DynamicView>::Flags & PacketAccessBit));
	using SumView = decltype((std::declval<Vector>() + std::declval<Vector>()).real());
	STATIC_CHECK(!(internal::evaluator<SumView>::Flags & PacketAccessBit));

	const Real nan = NumTraits<Real>::quiet_NaN();
	const Real inf = NumTraits<Real>::infinity();
	constexpr Index packet_size = internal::packet_traits<Real>::size;
	for (Index size : {Index(1), packet_size, packet_size + 1, 2 * packet_size - 1, Index(65)}) {
	Vector buffer(2 * size + 1);
	buffer.setConstant(Complex(Real(7), Real(-3)));
	Map<Vector> values(buffer.data() + 1, size);
	const auto& const_values = values;
	check_component(values.real());
	check_component(values.imag());
	Strided strided(buffer.data() + 1, size);
	check_component(strided.real());
	check_component(strided.imag());
	DynamicStride dynamic(buffer.data() + 1, size, InnerStride<Dynamic>(2));
	check_component(dynamic.real());
	check_component(dynamic.imag());
	for (Index i = 0; i < size; ++i) {
	for (Real special : {nan, inf, -inf, Real(0), -Real(0), std::numeric_limits<Real>::denorm_min()}) {
	values.setConstant(Complex(Real(7), Real(-3)));
	values(i) = Complex(special, Real(-3));
	Matrix<Real, Dynamic, 1> real = const_values.real();
	if ((numext::isnan)(special)) {
	VERIFY((numext::isnan)(real(i)));
	VERIFY((numext::isnan)(values.real().template maxCoeff<PropagateNaN>()));
	VERIFY((numext::isnan)(const_values.real().template minCoeff<PropagateNaN>()));
	if (size > 1) VERIFY_IS_EQUAL(values.real().template maxCoeff<PropagateNumbers>(), Real(7));
	} else {
	VERIFY_IS_EQUAL(real(i), special);
	VERIFY_IS_EQUAL((std::signbit)(real(i)), (std::signbit)(special));
	}
	VERIFY_IS_EQUAL(values.imag().template maxCoeff<PropagateNaN>(), Real(-3));
	values(i) = Complex(Real(7), special);
	Matrix<Real, Dynamic, 1> imag = const_values.imag();
	if ((numext::isnan)(special)) {
	VERIFY((numext::isnan)(imag(i)));
	VERIFY((numext::isnan)(values.imag().template minCoeff<PropagateNaN>()));
	VERIFY((numext::isnan)(const_values.imag().template maxCoeff<PropagateNaN>()));
	if (size > 1) VERIFY_IS_EQUAL(values.imag().template minCoeff<PropagateNumbers>(), Real(-3));
	} else {
	VERIFY_IS_EQUAL(imag(i), special);
	VERIFY_IS_EQUAL((std::signbit)(imag(i)), (std::signbit)(special));
	}
	VERIFY_IS_EQUAL(values.real().template minCoeff<PropagateNaN>(), Real(7));
	}
	}
	}
	}

	void component_custom_scalar() {
	using Vector = Matrix<PaddedComplex, Dynamic, 1>;
	Vector values(17);
	for (Index i = 0; i < values.size(); ++i) values(i) = PaddedComplex(double(i), double(2 * i + 1));
	const Vector& const_values = values;
	using RealView = decltype(const_values.real());
	using ImagView = decltype(const_values.imag());
	STATIC_CHECK(!(internal::evaluator<RealView>::Flags & PacketAccessBit));
	STATIC_CHECK(!(internal::evaluator<ImagView>::Flags & PacketAccessBit));
	check_component(const_values.real());
	check_component(const_values.imag());
	}

	EIGEN_DECLARE_TEST(complex_components) {
	CALL_SUBTEST_1(component_reductions<float>());
	CALL_SUBTEST_2(component_reductions<double>());
	CALL_SUBTEST_3(component_reductions<long double>());
	CALL_SUBTEST_3(component_custom_scalar());
	for (Index rows : {Index(2), Index(7), Index(17)}) {
	for (Index cols : {Index(3), Index(9)}) {
	EIGEN_UNUSED_VARIABLE(rows);
	EIGEN_UNUSED_VARIABLE(cols);
	CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, Dynamic, Dynamic, ColMajor>>(rows, cols)));
	CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>(rows, cols)));
	CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, Dynamic, Dynamic, ColMajor>>(rows, cols)));
	CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>(rows, cols)));
	}
	}
	CALL_SUBTEST_1((component_expressions<Matrix<std::complex<float>, 3, 5>>(3, 5)));
	CALL_SUBTEST_2((component_expressions<Matrix<std::complex<double>, 2, 2>>(2, 2)));
	}