// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2015-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

#define EIGEN_TEST_NO_LONGDOUBLE
#define EIGEN_DEFAULT_DENSE_INDEX_TYPE int

#define EIGEN_USE_GPU
#include "main.h"
#include "gpu_common.h"

// Check that dense modules can be properly parsed by nvcc
#include <Eigen/Dense>

// struct Foo{
//   EIGEN_DEVICE_FUNC
//   void operator()(int i, const float* mats, float* vecs) const {
//     using namespace Eigen;
//   //   Matrix3f M(data);
//   //   Vector3f x(data+9);
//   //   Map<Vector3f>(data+9) = M.inverse() * x;
//     Matrix3f M(mats+i/16);
//     Vector3f x(vecs+i*3);
//   //   using std::min;
//   //   using std::sqrt;
//     Map<Vector3f>(vecs+i*3) << x.minCoeff(), 1, 2;// / x.dot(x);//(M.inverse() *  x) / x.x();
//     //x = x*2 + x.y() * x + x * x.maxCoeff() - x / x.sum();
//   }
// };

template <typename MatrixType, int UpLo>
struct selfadjoint_permutation {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename MatrixType::Scalar* in,
                                    typename MatrixType::Scalar* out) const {
    using namespace Eigen;
    MatrixType matrix = Map<const MatrixType>(in + i);
    PermutationMatrix<4> permutation;
    permutation.indices() << 1, 3, 0, 2;
    matrix.template selfadjointView<UpLo>().applyPermutation(permutation);
    Map<MatrixType>(out + 48 * i) = matrix;
    matrix.template selfadjointView<UpLo>().applyPermutation(permutation.inverse());
    Map<MatrixType>(out + 48 * i + 16) = matrix;
    matrix.template selfadjointView<UpLo>().applyTransposition(3, 0);
    Map<MatrixType>(out + 48 * i + 32) = matrix;
  }
};

template <typename Scalar>
struct scaled_permutation_product {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    using MatrixType = Eigen::Matrix<Scalar, 3, 3>;
    Eigen::ScaledPermutationMatrix<Scalar, 3> scaled;
    scaled.setIdentity();
    scaled.indices() << 1, 2, 0;
    scaled.scales() << Scalar(2), Scalar(3), Scalar(-1);
    const MatrixType matrix = Eigen::Map<const Eigen::Matrix3f>(in + i).template cast<Scalar>();
    MatrixType result = scaled * matrix;
    result.noalias() += matrix * scaled;
    result.noalias() -= scaled * matrix;
    result = scaled * result;
    result = result * scaled;
    result += scaled;
    result -= scaled;
    result += scaled.toDenseMatrix();
    result = result + scaled;
    result.noalias() += scaled * matrix.template triangularView<Eigen::Lower>();
    result.noalias() += matrix.template selfadjointView<Eigen::Upper>() * scaled;
    Eigen::Map<Eigen::Matrix3f>(out + 9 * i) = result.template cast<float>();
  }
};

template <typename T>
struct coeff_wise {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    T x1(in + i);
    T x2(in + i + 1);
    T x3(in + i + 2);
    Map<T> res(out + i * T::MaxSizeAtCompileTime);

    res.array() += (in[0] * x1 + x2).array() * x3.array();
  }
};

// The static Map() helpers of the plain types are device-callable.
template <typename T>
struct static_map {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    T::Map(out + i * T::SizeAtCompileTime) =
        2 * T::Map(in + i) + T::Map(in + i, T::RowsAtCompileTime, T::ColsAtCompileTime);
  }
};

// A transform times a diagonal matrix scales the linear part through MatrixBase::operator*=.
struct transform_times_diagonal {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    using namespace Eigen;
    Affine3f T;
    T.matrix() = Matrix4f(in + i);
    Map<Matrix4f>(out + i * 16) = (T * DiagonalMatrix<float, 3>(in[i], in[i + 1], in[i + 2])).matrix();
  }
};

// The Scaling() helpers and UniformScaling compose with transforms on the device.
struct transform_scaling {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    using namespace Eigen;
    Affine3f T = Translation3f(Vector3f(in + i)) * Scaling(in[i + 3]);
    T *= Scaling(in[i + 4]).inverse();
    Map<Matrix4f>(out + i * 16) = (T * Scaling(in[i + 5], in[i + 6], in[i + 7])).matrix();
  }
};

template <int Order>
struct scaled_structured_product {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    using Mat = Eigen::Matrix<float, 3, 3, Order>;
    const Eigen::Map<const Mat> matrix(in + i);
    Eigen::Map<Mat> result(out + i * 9);
    result.noalias() = 3.0f * (matrix.template triangularView<Eigen::UnitLower>() * matrix.diagonal().asDiagonal());
    result += 2.0f * (matrix.diagonal().asDiagonal() * matrix) + Mat::Zero();
  }
};

template <int Order>
struct scaled_outer_product {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    // More than 16 rows select the scaled outer-product functor.
    using Lhs = Eigen::Matrix<float, 17, 1>;
    using Rhs = Eigen::RowVector3f;
    using ProductImpl =
        Eigen::internal::generic_product_impl<Lhs, Rhs, Eigen::DenseShape, Eigen::DenseShape, Eigen::OuterProduct>;
    const Lhs lhs(in + i);
    const Rhs rhs(in + i + 17);
    Eigen::Map<Eigen::Matrix<float, 17, 3, Order>> result(out + i * 51);
    result.setZero();
    ProductImpl::scaleAndAddTo(result, lhs, rhs, 2.0f);
  }
};

template <int Order>
struct mixed_outer_product {
  EIGEN_DEVICE_FUNC void operator()(int i, const std::complex<float>* in, std::complex<float>* out) const {
    using Mat = Eigen::Matrix<std::complex<float>, 4, 4, Order>;
    const Eigen::Map<const Eigen::Vector4cf> complex(in + i);
    const Eigen::Vector4f real = complex.real();
    Eigen::Map<Mat> result(out + i * 16);
    if (Order == Eigen::RowMajor) {
      result.noalias() = complex * real.transpose();
      result += complex.lazyProduct(real.transpose());
    } else {
      result.noalias() = real * complex.transpose();
      result += real.lazyProduct(complex.transpose());
    }
  }
};

template <int Order>
struct scaled_selfadjoint_diagonal_product {
  EIGEN_DEVICE_FUNC void operator()(int i, const std::complex<float>* in, std::complex<float>* out) const {
    using Scalar = std::complex<float>;
    using Mat = Eigen::Matrix<Scalar, 3, 3, Order>;
    const Eigen::Map<const Mat> matrix(in + i);
    const Eigen::Vector3cf diagonal = Eigen::Vector3cf::Constant(Scalar(2, 1));
    const Scalar alpha(2, 3);
    Eigen::Map<Mat> result(out + i * 9);
    result = (alpha * matrix.template selfadjointView<Eigen::Lower>()) * diagonal.asDiagonal() + Mat::Zero();
    result += diagonal.asDiagonal() * (alpha * matrix.template selfadjointView<Eigen::Upper>()) + Mat::Zero();
  }
};

struct make_householder_small_tail {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* /*in*/, float* out) const {
    Eigen::Vector3f vector;
    vector << 0.0f, 1e-20f, -2e-20f;
    Eigen::Vector2f essential;
    float tau;
    float beta;
    vector.makeHouseholder(essential, tau, beta);
    out[4 * i] = tau;
    out[4 * i + 1] = beta;
    out[4 * i + 2] = essential[0];
    out[4 * i + 3] = essential[1];
  }
};

struct make_householder_complex_zero_tail {
  EIGEN_DEVICE_FUNC void operator()(int i, const std::complex<float>* /*in*/, std::complex<float>* out) const {
    Eigen::Vector2cf vector;
    vector << std::complex<float>(0.0f, 1e-20f), std::complex<float>(0.0f, 0.0f);
    Eigen::Matrix<std::complex<float>, 1, 1> essential;
    std::complex<float> tau;
    float beta;
    vector.makeHouseholder(essential, tau, beta);
    out[3 * i] = tau;
    out[3 * i + 1] = std::complex<float>(beta, 0.0f);
    out[3 * i + 2] = essential[0];
  }
};

// Covers each unitOrthogonal path: fixed size 2, both branches of fixed size 3, and the generic one.
struct unit_orthogonal {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    using namespace Eigen;
    Map<Vector2f>(out + 12 * i) = Map<const Vector2f>(in + i).unitOrthogonal();
    Map<Vector3f>(out + 12 * i + 2) = Vector3f(in[i] + 2.0f, in[i + 1], in[i + 2]).unitOrthogonal();
    Map<Vector3f>(out + 12 * i + 5) = Vector3f(0.0f, 1e-20f * in[i], in[i + 1] + 2.0f).unitOrthogonal();
    Map<Vector4f>(out + 12 * i + 8) = Map<const Vector4f>(in + i).unitOrthogonal();
  }
};

// Applies the complex operators inside Eigen's own templates, which see the device overloads only through
// Eigen/Core's include order; complex_operators below finds them through its using-directive.
template <typename ComplexType>
struct complex_internal_operators {
  EIGEN_DEVICE_FUNC void operator()(int i, const ComplexType* in, ComplexType* out) const {
    const int num_operators = 8;
    int out_idx = i * num_operators;
    const ComplexType a = in[i];
    const ComplexType b = in[i + 1];

    out[out_idx++] = numext::negate(a);
    out[out_idx++] = numext::conj(a);
    out[out_idx++] = internal::padd(a, b);
    out[out_idx++] = internal::psub(a, b);
    out[out_idx++] = internal::pmul(a, b);
    out[out_idx++] = internal::pdiv(a, b);
    out[out_idx++] = internal::pnegate(a);
    out[out_idx++] = internal::pconj(a);
  }
};

template <typename T>
struct complex_sqrt {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    typedef typename T::Scalar ComplexType;
    typedef typename T::Scalar::value_type ValueType;
    const int num_special_inputs = 18;

    if (i == 0) {
      const ValueType nan = std::numeric_limits<ValueType>::quiet_NaN();
      typedef Eigen::Vector<ComplexType, num_special_inputs> SpecialInputs;
      SpecialInputs special_in;
      special_in.setZero();
      int idx = 0;
      special_in[idx++] = ComplexType(0, 0);
      special_in[idx++] = ComplexType(-0, 0);
      special_in[idx++] = ComplexType(0, -0);
      special_in[idx++] = ComplexType(-0, -0);
// GCC's fallback sqrt implementation fails for inf inputs.
// It is called when _GLIBCXX_USE_C99_COMPLEX is false or if
// clang includes the GCC header (which temporarily disables
// _GLIBCXX_USE_C99_COMPLEX)
#if !defined(_GLIBCXX_COMPLEX) || (_GLIBCXX_USE_C99_COMPLEX && !defined(__CLANG_CUDA_WRAPPERS_COMPLEX))
      const ValueType inf = std::numeric_limits<ValueType>::infinity();
      special_in[idx++] = ComplexType(1.0, inf);
      special_in[idx++] = ComplexType(nan, inf);
      special_in[idx++] = ComplexType(1.0, -inf);
      special_in[idx++] = ComplexType(nan, -inf);
      special_in[idx++] = ComplexType(-inf, 1.0);
      special_in[idx++] = ComplexType(inf, 1.0);
      special_in[idx++] = ComplexType(-inf, -1.0);
      special_in[idx++] = ComplexType(inf, -1.0);
      special_in[idx++] = ComplexType(-inf, nan);
      special_in[idx++] = ComplexType(inf, nan);
#endif
      special_in[idx++] = ComplexType(1.0, nan);
      special_in[idx++] = ComplexType(nan, 1.0);
      special_in[idx++] = ComplexType(nan, -1.0);
      special_in[idx++] = ComplexType(nan, nan);

      Map<SpecialInputs> special_out(out);
      special_out = special_in.cwiseSqrt();
    }

    T x1(in + i);
    Map<T> res(out + num_special_inputs + i * T::MaxSizeAtCompileTime);
    res = x1.cwiseSqrt();
  }
};

namespace complex_scalar_operators_test {
#if defined(EIGEN_GPUCC) && defined(EIGEN_GPU_COMPILE_PHASE)
// Shadow the global `using namespace Eigen;` from main.h so that unqualified
// lookup in this namespace only finds the device std::complex operators via
// EIGEN_USING_STD_COMPLEX_OPERATORS.
struct shadow_tag {};
EIGEN_DEVICE_FUNC void operator+(shadow_tag);
EIGEN_DEVICE_FUNC void operator-(shadow_tag);
EIGEN_DEVICE_FUNC void operator*(shadow_tag, shadow_tag);
EIGEN_DEVICE_FUNC void operator/(shadow_tag, shadow_tag);
EIGEN_DEVICE_FUNC void operator+=(shadow_tag, shadow_tag);
EIGEN_DEVICE_FUNC void operator-=(shadow_tag, shadow_tag);
EIGEN_DEVICE_FUNC void operator*=(shadow_tag, shadow_tag);
EIGEN_DEVICE_FUNC void operator/=(shadow_tag, shadow_tag);
EIGEN_DEVICE_FUNC void operator==(shadow_tag, shadow_tag);
EIGEN_DEVICE_FUNC void operator!=(shadow_tag, shadow_tag);
EIGEN_USING_STD_COMPLEX_OPERATORS
#endif

template <typename ComplexType>
EIGEN_DEVICE_FUNC void run(const ComplexType& a, const ComplexType& b, ComplexType* out, int& out_idx) {
  typedef typename ComplexType::value_type ValueType;
  const ValueType real_a = Eigen::numext::real(a);
  const ValueType real_b = Eigen::numext::real(b);

  out[out_idx++] = +a;
  out[out_idx++] = -a;

  out[out_idx++] = a + b;
  out[out_idx++] = a + real_b;
  out[out_idx++] = real_a + b;
  out[out_idx++] = a - b;
  out[out_idx++] = a - real_b;
  out[out_idx++] = real_a - b;
  out[out_idx++] = a * b;
  out[out_idx++] = a * real_b;
  out[out_idx++] = real_a * b;
  out[out_idx++] = a / b;
  out[out_idx++] = a / real_b;
  out[out_idx++] = real_a / b;

#if !EIGEN_COMP_MSVC
  out[out_idx] = a;
  out[out_idx++] += b;
  out[out_idx] = a;
  out[out_idx++] -= b;
  out[out_idx] = a;
  out[out_idx++] *= b;
  out[out_idx] = a;
  out[out_idx++] /= b;
#endif

  const ComplexType true_value = ComplexType(ValueType(1), ValueType(0));
  const ComplexType false_value = ComplexType(ValueType(0), ValueType(0));
  out[out_idx++] = (a == b ? true_value : false_value);
  out[out_idx++] = (a == real_b ? true_value : false_value);
  out[out_idx++] = (real_a == b ? true_value : false_value);
  out[out_idx++] = (a != b ? true_value : false_value);
  out[out_idx++] = (a != real_b ? true_value : false_value);
  out[out_idx++] = (real_a != b ? true_value : false_value);
}
}  // namespace complex_scalar_operators_test

template <typename T>
struct complex_operators {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    typedef typename T::Scalar ComplexType;
    typedef typename T::Scalar::value_type ValueType;
    const int num_scalar_operators = 24;
    const int num_vector_operators = 23;  // no unary + operator.
    int out_idx = i * (num_scalar_operators + num_vector_operators * T::MaxSizeAtCompileTime);

    // Scalar operators.
    const ComplexType a = in[i];
    const ComplexType b = in[i + 1];
    complex_scalar_operators_test::run(a, b, out, out_idx);

    const ComplexType true_value = ComplexType(ValueType(1), ValueType(0));
    const ComplexType false_value = ComplexType(ValueType(0), ValueType(0));

    // Vector versions.
    T x1(in + i);
    T x2(in + i + 1);
    const int res_size = T::MaxSizeAtCompileTime * num_scalar_operators;
    const int size = T::MaxSizeAtCompileTime;
    int block_idx = 0;

    Map<VectorX<ComplexType>> res(out + out_idx, res_size);
    res.segment(block_idx, size) = -x1;
    block_idx += size;

    res.segment(block_idx, size) = x1 + x2;
    block_idx += size;
    res.segment(block_idx, size) = x1 + x2.real();
    block_idx += size;
    res.segment(block_idx, size) = x1.real() + x2;
    block_idx += size;
    res.segment(block_idx, size) = x1 - x2;
    block_idx += size;
    res.segment(block_idx, size) = x1 - x2.real();
    block_idx += size;
    res.segment(block_idx, size) = x1.real() - x2;
    block_idx += size;
    res.segment(block_idx, size) = x1.array() * x2.array();
    block_idx += size;
    res.segment(block_idx, size) = x1.array() * x2.real().array();
    block_idx += size;
    res.segment(block_idx, size) = x1.real().array() * x2.array();
    block_idx += size;
    res.segment(block_idx, size) = x1.array() / x2.array();
    block_idx += size;
    res.segment(block_idx, size) = x1.array() / x2.real().array();
    block_idx += size;
    res.segment(block_idx, size) = x1.real().array() / x2.array();
    block_idx += size;

#if !EIGEN_COMP_MSVC
    res.segment(block_idx, size) = x1;
    res.segment(block_idx, size) += x2;
    block_idx += size;
    res.segment(block_idx, size) = x1;
    res.segment(block_idx, size) -= x2;
    block_idx += size;
    res.segment(block_idx, size) = x1;
    res.segment(block_idx, size).array() *= x2.array();
    block_idx += size;
    res.segment(block_idx, size) = x1;
    res.segment(block_idx, size).array() /= x2.array();
    block_idx += size;
#endif

    const T true_vector = T::Constant(true_value);
    const T false_vector = T::Constant(false_value);
    res.segment(block_idx, size) = (x1 == x2 ? true_vector : false_vector);
    block_idx += size;
    // Mixing types in equality comparison does not work.
    // res.segment(block_idx, size) = (x1 == x2.real() ? true_vector : false_vector);
    // block_idx += size;
    // res.segment(block_idx, size) = (x1.real() == x2 ? true_vector : false_vector);
    // block_idx += size;
    res.segment(block_idx, size) = (x1 != x2 ? true_vector : false_vector);
    block_idx += size;
    // res.segment(block_idx, size) = (x1 != x2.real() ? true_vector : false_vector);
    // block_idx += size;
    // res.segment(block_idx, size) = (x1.real() != x2 ? true_vector : false_vector);
    // block_idx += size;
  }
};

template <typename T>
struct replicate {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    T x1(in + i);
    int step = x1.size() * 4;
    int stride = 3 * step;

    typedef Map<Array<typename T::Scalar, Dynamic, Dynamic>> MapType;
    MapType(out + i * stride + 0 * step, x1.rows() * 2, x1.cols() * 2) = x1.replicate(2, 2);
    MapType(out + i * stride + 1 * step, x1.rows() * 3, x1.cols()) = in[i] * x1.colwise().replicate(3);
    MapType(out + i * stride + 2 * step, x1.rows(), x1.cols() * 3) = in[i] * x1.rowwise().replicate(3);
  }
};

template <typename T>
struct alloc_new_delete {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    int offset = 2 * i * T::MaxSizeAtCompileTime;
    T* x = new T(in + offset);
    Eigen::Map<T> u(out + offset);
    u = *x;
    delete x;

    offset += T::MaxSizeAtCompileTime;
    T* y = new T[1];
    y[0] = T(in + offset);
    Eigen::Map<T> v(out + offset);
    v = y[0];
    delete[] y;
  }
};

template <typename T>
struct redux {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    int N = 10;
    T x1(in + i);
    out[i * N + 0] = x1.minCoeff();
    out[i * N + 1] = x1.maxCoeff();
    out[i * N + 2] = x1.sum();
    out[i * N + 3] = x1.prod();
    out[i * N + 4] = x1.matrix().squaredNorm();
    out[i * N + 5] = x1.matrix().norm();
    out[i * N + 6] = x1.colwise().sum().maxCoeff();
    out[i * N + 7] = x1.rowwise().maxCoeff().sum();
    out[i * N + 8] = x1.matrix().colwise().squaredNorm().sum();
  }
};

template <typename T1, typename T2>
struct prod_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T1::Scalar* in, typename T1::Scalar* out) const {
    using namespace Eigen;
    typedef Matrix<typename T1::Scalar, T1::RowsAtCompileTime, T2::ColsAtCompileTime> T3;
    T1 x1(in + i);
    T2 x2(in + i + 1);
    Map<T3> res(out + i * T3::MaxSizeAtCompileTime);
    res += in[i] * x1 * x2;
  }
};

template <typename T>
struct compound_product {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using Scalar = typename T::Scalar;
    using ProdEval =
        Eigen::internal::product_evaluator<Eigen::Product<T, T, Eigen::LazyProduct>, Eigen::CoeffBasedProductMode,
                                           Eigen::DenseShape, Eigen::DenseShape>;
    using AddOp = Eigen::internal::add_assign_op<Scalar, Scalar>;
    static_assert(!Eigen::internal::product_packet_cascade_traits<ProdEval, T, AddOp>::Enable,
                  "product packet cascade must stay disabled when ActualPacketAccessBit == 0");
    constexpr int size = T::SizeAtCompileTime;
    const T lhs(in + i);
    const T rhs(in + i + 1);
    const T addend(in + i + 2);
    T accum = T::Zero();
    accum.noalias() += lhs * rhs;
    accum.noalias() -= lhs * rhs;
    accum += addend + lhs * rhs;
#if defined(EIGEN_GPU_COMPILE_PHASE)
    using Packet = typename Eigen::internal::packet_traits<Scalar>::type;
    if (size >= Eigen::internal::unpacket_traits<Packet>::size) {
      const Packet zero = Eigen::internal::pzero(Packet());
      AddOp().template assignPacket<Eigen::Unaligned, Packet>(accum.data(), zero);
    }
#endif
    Eigen::Map<T>(out + i * size) = accum;
  }
};

template <typename Scalar, int Options>
void test_compound_product() {
  constexpr int n = 4;
  Eigen::Array<Scalar, Eigen::Dynamic, 1> in(n * 512), out(n * 512);
  in.setRandom();
  out.setZero();

  run_and_compare_to_gpu(compound_product<Eigen::Matrix<Scalar, 3, 3, Options>>(), n, in, out);
  run_and_compare_to_gpu(compound_product<Eigen::Matrix<Scalar, 5, 5, Options>>(), n, in, out);
  run_and_compare_to_gpu(compound_product<Eigen::Matrix<Scalar, 9, 9, Options>>(), n, in, out);
}

template <typename T1, typename T2>
struct diagonal {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T1::Scalar* in, typename T1::Scalar* out) const {
    using namespace Eigen;
    T1 x1(in + i);
    Map<T2> res(out + i * T2::MaxSizeAtCompileTime);
    res += x1.diagonal();
  }
};

// The copy constructors of the structured plain types are device-callable.
struct structured_copies {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    using namespace Eigen;
    const DiagonalMatrix<float, 3> d(in[i], in[i + 1], in[i + 2]);
    const DiagonalMatrix<float, 3> diagonalCopy(d);
    const SkewSymmetricMatrix3<float> s(in[i + 3], in[i + 4], in[i + 5]);
    const SkewSymmetricMatrix3<float> skewCopy(s);
    Map<Vector3f>(out + i * 6) = diagonalCopy.diagonal();
    Map<Vector3f>(out + i * 6 + 3) = skewCopy.vector();
  }
};

template <typename T>
struct jacobi_rotations {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using Scalar = typename T::Scalar;
    constexpr int size = T::SizeAtCompileTime;
    const Eigen::JacobiRotation<Scalar> rotation(Scalar(3) / Scalar(5), Scalar(4) / Scalar(5));

    T fixed(in + i);
    fixed.applyOnTheLeft(0, 1, rotation);
    fixed.applyOnTheRight(0, 1, rotation);
    Eigen::Map<T>(out + 2 * i * size) = fixed;

    using DynamicMatrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, T::Options>;
    Eigen::Map<DynamicMatrix> dynamic(out + (2 * i + 1) * size, T::RowsAtCompileTime, T::ColsAtCompileTime);
    dynamic = T(in + i);
    dynamic.applyOnTheLeft(0, 1, rotation);
    dynamic.applyOnTheRight(0, 1, rotation);
  }
};

template <typename Scalar, int Options>
void test_jacobi_rotations() {
  constexpr int n = 4;
  Eigen::Array<Scalar, Eigen::Dynamic, 1> in(n * 512), out(n * 512);
  in.setRandom();
  out.setZero();

  // Sizes 3, 4 and 9 cover scalar fallback, aligned packets, and dynamic packets with a tail.
  run_and_compare_to_gpu(jacobi_rotations<Eigen::Matrix<Scalar, 3, 3, Options>>(), n, in, out);
  run_and_compare_to_gpu(jacobi_rotations<Eigen::Matrix<Scalar, 4, 4, Options>>(), n, in, out);
  run_and_compare_to_gpu(jacobi_rotations<Eigen::Matrix<Scalar, 9, 9, Options>>(), n, in, out);
}

template <typename T>
struct eigenvalues_direct {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    typedef Matrix<typename T::Scalar, T::RowsAtCompileTime, 1> Vec;
    T M(in + i);
    Map<Vec> res(out + i * Vec::MaxSizeAtCompileTime);
    T A = M * M.adjoint();
    SelfAdjointEigenSolver<T> eig;
    eig.computeDirect(A);
    res = eig.eigenvalues();
  }
};

template <typename T>
struct eigenvalues {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    typedef Matrix<typename T::Scalar, T::RowsAtCompileTime, 1> Vec;
    T M(in + i);
    Map<Vec> res(out + i * Vec::MaxSizeAtCompileTime);
    T A = M * M.adjoint();
    SelfAdjointEigenSolver<T> eig;
    eig.compute(A);
    res = eig.eigenvalues();
  }
};

template <typename T, int UpLo>
struct selfadjoint_rank2_update {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    typedef Matrix<typename T::Scalar, T::RowsAtCompileTime, 1> Vec;
    T M(in + i);
    Vec u(in + i + T::MaxSizeAtCompileTime);
    Vec v(in + i + T::MaxSizeAtCompileTime + Vec::MaxSizeAtCompileTime);
    Map<T> res(out + i * T::MaxSizeAtCompileTime);
    res = M;
    res.template selfadjointView<UpLo>().rankUpdate(u, v, typename T::Scalar(0.25));
  }
};

template <typename T, int UpLo>
struct selfadjoint_l1_norm {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    T M(in + i);
    // l1Norm() has a separate device implementation, so the host result is the reference.
    out[i] = M.template selfadjointView<UpLo>().l1Norm();
  }
};

template <typename T>
struct matrix_inverse {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    using namespace Eigen;
    T M(in + i);
    Map<T> res(out + i * T::MaxSizeAtCompileTime);
    res = M.inverse();
  }
};

template <typename T>
struct numeric_limits_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    EIGEN_UNUSED_VARIABLE(in);
    int out_idx = i * 5;
    out[out_idx++] = numext::numeric_limits<float>::epsilon();
    out[out_idx++] = (numext::numeric_limits<float>::max)();
    out[out_idx++] = (numext::numeric_limits<float>::min)();
    out[out_idx++] = numext::numeric_limits<float>::infinity();
    out[out_idx++] = numext::numeric_limits<float>::quiet_NaN();
  }
};

struct custom_less_scalar {
  int value;

  EIGEN_DEVICE_FUNC explicit custom_less_scalar(int x = 0) : value(x) {}
};

EIGEN_DEVICE_FUNC bool operator<(const custom_less_scalar& x, const custom_less_scalar& y) { return x.value < y.value; }

struct custom_less_scalar_minmax_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const int* in, int* out) const {
    EIGEN_UNUSED_VARIABLE(i);
    const custom_less_scalar x(in[0]);
    const custom_less_scalar y(in[1]);
    out[0] = Eigen::numext::mini(x, y).value;
    out[1] = Eigen::numext::maxi(x, y).value;
  }
};

void test_custom_less_scalar_minmax() {
  Eigen::ArrayXi in(2), out_ref(2), out_gpu(2);
  in << 1, 2;
  out_ref.setConstant(-1);
  out_gpu.setConstant(-1);

  run_on_cpu(custom_less_scalar_minmax_test(), 1, in, out_ref);
  run_on_gpu(custom_less_scalar_minmax_test(), 1, in, out_gpu);

#if !defined(EIGEN_GPU_COMPILE_PHASE)
  VERIFY_IS_EQUAL(out_ref(0), out_gpu(0));
  VERIFY_IS_EQUAL(out_ref(1), out_gpu(1));
#endif
}

struct float_nan_minmax_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    EIGEN_UNUSED_VARIABLE(i);
    const float nan = in[0];
    const float one = in[1];
    out[0] = Eigen::numext::mini(nan, one);
    out[1] = Eigen::numext::mini(one, nan);
    out[2] = Eigen::numext::maxi(nan, one);
    out[3] = Eigen::numext::maxi(one, nan);
    out[4] = Eigen::numext::mini(nan, nan);
    out[5] = Eigen::numext::maxi(nan, nan);
  }
};

void test_float_nan_minmax() {
  Eigen::ArrayXf in(2), out_ref(6), out_gpu(6);
  in << std::numeric_limits<float>::quiet_NaN(), 1.f;
  out_ref.setConstant(-1.f);
  out_gpu.setConstant(-1.f);

  run_on_cpu(float_nan_minmax_test(), 1, in, out_ref);
  run_on_gpu(float_nan_minmax_test(), 1, in, out_gpu);

#if !defined(EIGEN_GPU_COMPILE_PHASE)
  VERIFY_IS_CWISE_EQUAL(out_ref, out_gpu);
#if defined(EIGEN_CONSTEXPR_ARE_DEVICE_FUNC)
  VERIFY((numext::isnan)(out_ref(0)));
  VERIFY((numext::isnan)(out_ref(2)));
#else
  // Without relaxed constexpr, numext::mini/maxi select the number-preferring fmin/fmax overloads.
  VERIFY_IS_EQUAL(out_ref(0), 1.f);
  VERIFY_IS_EQUAL(out_ref(2), 1.f);
#endif
  VERIFY_IS_EQUAL(out_ref(1), 1.f);
  VERIFY_IS_EQUAL(out_ref(3), 1.f);
  VERIFY((numext::isnan)(out_ref(4)));
  VERIFY((numext::isnan)(out_ref(5)));
#endif
}

// clz, ctz and popcount of the low 8, 16, 32 and 64 bits of each input. The clz test checks the host counts.
struct bit_count_test {
  template <typename T>
  EIGEN_DEVICE_FUNC static void count(T bits, int* out) {
    out[0] = internal::clz(bits);
    out[1] = internal::ctz(bits);
    out[2] = internal::popcount(bits);
  }

  EIGEN_DEVICE_FUNC void operator()(int i, const uint64_t* in, int* out) const {
    const uint64_t bits = in[i];
    count(static_cast<uint8_t>(bits), out + 12 * i);
    count(static_cast<uint16_t>(bits), out + 12 * i + 3);
    count(static_cast<uint32_t>(bits), out + 12 * i + 6);
    count(bits, out + 12 * i + 9);
  }
};

void test_bit_count() {
  const int n = 3 * 64 + 2;
  Eigen::Array<uint64_t, Eigen::Dynamic, 1> in(n);
  Eigen::ArrayXi out_ref(12 * n), out_gpu(12 * n);
  for (int k = 0; k < 64; ++k) {
    in(3 * k) = uint64_t(1) << k;
    in(3 * k + 1) = ~(uint64_t(1) << k);
    in(3 * k + 2) = internal::random<uint64_t>();
  }
  in(n - 2) = 0;
  in(n - 1) = ~uint64_t(0);
  out_ref.setConstant(-1);
  out_gpu.setConstant(-1);

  run_on_cpu(bit_count_test(), n, in, out_ref);
  run_on_gpu(bit_count_test(), n, in, out_gpu);

#if !defined(EIGEN_GPU_COMPILE_PHASE)
  VERIFY_IS_CWISE_EQUAL(out_ref, out_gpu);
#endif
}

// Eigen::half != is true when either operand is NaN, on the device as on the host.
struct half_nan_not_equal_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    EIGEN_UNUSED_VARIABLE(i);
    const Eigen::half nan(in[0]), one(in[1]);
    out[0] = float(nan != nan);
    out[1] = float(nan != one);
    out[2] = float(one != nan);
    out[3] = float(one != one);
  }
};

void test_half_nan_not_equal() {
  Eigen::ArrayXf in(2), out_ref(4), out_gpu(4);
  in << std::numeric_limits<float>::quiet_NaN(), 1.f;
  run_on_cpu(half_nan_not_equal_test(), 1, in, out_ref);
  run_on_gpu(half_nan_not_equal_test(), 1, in, out_gpu);
#if !defined(EIGEN_GPU_COMPILE_PHASE)
  VERIFY_IS_CWISE_EQUAL(out_gpu, Eigen::Array4f(1.f, 1.f, 1.f, 0.f));
  VERIFY_IS_CWISE_EQUAL(out_ref, out_gpu);
#endif
}

// exp of a real-axis complex keeps the zero imaginary part, and its sign, where the real part overflows.
struct complex_exp_overflow_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* in, float* out) const {
    EIGEN_UNUSED_VARIABLE(i);
    const std::complex<float> e = numext::exp(std::complex<float>(in[0], 0.f));
    const std::complex<double> e2 = numext::exp(std::complex<double>(in[0] * 10, -0.0));
    const std::complex<double> p2 = numext::exp2(std::complex<double>(in[0] * 20, 0.0));
    out[0] = e.imag();
    out[1] = float(e2.imag());
    out[2] = float(p2.imag());
  }
};

void test_complex_exp_overflow() {
  Eigen::ArrayXf in(1), out_ref(3), out_gpu(3);
  in << 100.f;
  run_on_cpu(complex_exp_overflow_test(), 1, in, out_ref);
  run_on_gpu(complex_exp_overflow_test(), 1, in, out_gpu);
#if !defined(EIGEN_GPU_COMPILE_PHASE)
  VERIFY_IS_CWISE_EQUAL(out_ref, Eigen::Array3f::Zero());
  VERIFY_IS_CWISE_EQUAL(out_gpu, Eigen::Array3f::Zero());
  VERIFY((std::signbit)(out_ref(1)) && (std::signbit)(out_gpu(1)));
#endif
}

template <typename Type1, typename Type2>
bool verifyIsApproxWithInfsNans(const Type1& a, const Type2& b,
                                typename Type1::Scalar* = 0)  // Enabled for Eigen's type only
{
  if (a.rows() != b.rows()) {
    return false;
  }
  if (a.cols() != b.cols()) {
    return false;
  }
  for (Index r = 0; r < a.rows(); ++r) {
    for (Index c = 0; c < a.cols(); ++c) {
      if (a(r, c) != b(r, c) && !((numext::isnan)(a(r, c)) && (numext::isnan)(b(r, c))) &&
          !test_isApprox(a(r, c), b(r, c))) {
        return false;
      }
    }
  }
  return true;
}

#if defined(EIGEN_HAS_GPU_FP16) && !defined(EIGEN_GPU_COMPILE_PHASE)
// Host-side check that converting between Eigen::half and the vendor __half type preserves the
// raw bits. This is a regression test for builds where Eigen::half stores a native fp16 type
// (e.g. __fp16 on arm64): the host phase used to perform numeric value conversions instead of
// bit reinterpretations, corrupting every raw-bit constant (NumTraits, numeric_limits, ...).
void test_half_raw_bit_interop() {
  const numext::uint16_t raw_bits[] = {0x0000, 0x3c00 /*1*/, 0x7c00 /*inf*/, 0x7e00 /*qNaN*/, 0xfbff /*lowest*/};
  for (int i = 0; i < 5; ++i) {
    const numext::uint16_t raw = raw_bits[i];
    const Eigen::half h = numext::bit_cast<Eigen::half>(raw);
    // Eigen::half -> __half must preserve the bits (sizeof(__half) == 2 on both CUDA and HIP).
    // Call the conversion operator explicitly: a static_cast would be ambiguous because
    // Eigen::half also converts to __half via operator float() and __half(float).
    const __half v = h.operator __half();
    VERIFY_IS_EQUAL(numext::bit_cast<numext::uint16_t>(v), raw);
    // __half -> Eigen::half must preserve the bits as well.
    const Eigen::half h2(v);
    VERIFY_IS_EQUAL(numext::bit_cast<numext::uint16_t>(h2), raw);
  }
  // Raw-bit constants must survive the host phase of a GPU build.
  VERIFY((numext::isinf)(NumTraits<Eigen::half>::infinity()));
  VERIFY((numext::isnan)(NumTraits<Eigen::half>::quiet_NaN()));
}
#endif

template <typename Kernel, typename Input, typename Output>
void test_with_infs_nans(const Kernel& ker, int n, const Input& in, Output& out) {
  Output out_ref, out_gpu;
#if !defined(EIGEN_GPU_COMPILE_PHASE)
  out_ref = out_gpu = out;
#else
  EIGEN_UNUSED_VARIABLE(in);
  EIGEN_UNUSED_VARIABLE(out);
#endif
  run_on_cpu(ker, n, in, out_ref);
  run_on_gpu(ker, n, in, out_gpu);
#if !defined(EIGEN_GPU_COMPILE_PHASE)
  verifyIsApproxWithInfsNans(out_ref, out_gpu);
#endif
}

// `preverse` falls back to returning the packet unchanged for any packet type
// that does not override it, which is only correct for a packet of one element.
// A GPU translation unit resolves both the host and the device side to the same
// packet type, so running the same reversal on both and comparing them cannot
// catch a missing override: check the reversed values instead.
template <typename T>
struct reverse_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const typename T::Scalar* in, typename T::Scalar* out) const {
    constexpr int size = T::SizeAtCompileTime;
    Eigen::Map<T>(out + i * size) = T(in + i * size).reverse();
  }
};

template <typename Scalar>
struct packet_any_test {
  EIGEN_DEVICE_FUNC void operator()(int i, const float* /*in*/, float* out) const {
#if defined(EIGEN_GPU_COMPILE_PHASE)
    using Packet = typename Eigen::internal::packet_traits<Scalar>::type;
    const Packet zero = Eigen::internal::pzero(Packet());
    const Packet sequence = Eigen::internal::plset<Packet>(Scalar(0));
    const Packet mixed = Eigen::internal::pcmp_eq(sequence, Eigen::internal::pset1<Packet>(Scalar(1)));
    out[3 * i] = float(Eigen::internal::predux_any(zero));
    out[3 * i + 1] = float(Eigen::internal::predux_any(Eigen::internal::ptrue(zero)));
    out[3 * i + 2] = float(Eigen::internal::predux_any(mixed));
#else
    out[3 * i] = 0;
    out[3 * i + 1] = 1;
    out[3 * i + 2] = 1;
#endif
  }
};

template <typename T>
void test_reverse() {
  typedef typename T::Scalar Scalar;
  constexpr int size = T::SizeAtCompileTime;
  constexpr int n = 4;

  Eigen::Array<Scalar, Eigen::Dynamic, 1> in(n * size), out_ref(n * size), out_gpu(n * size);
  for (int i = 0; i < n * size; ++i) in(i) = static_cast<Scalar>(i + 1);
  out_ref.setZero();
  out_gpu.setZero();

  run_on_cpu(reverse_test<T>(), n, in, out_ref);
  run_on_gpu(reverse_test<T>(), n, in, out_gpu);

#if !defined(EIGEN_GPU_COMPILE_PHASE)
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < size; ++j) {
      const Scalar expected = in(i * size + size - 1 - j);
      VERIFY_IS_EQUAL(out_ref(i * size + j), expected);
      VERIFY_IS_EQUAL(out_gpu(i * size + j), expected);
    }
  }
#endif
}

EIGEN_DECLARE_TEST(gpu_basic) {
  ei_test_init_gpu();

  int nthreads = 100;
  Eigen::VectorXf in, out;
  Eigen::VectorXcf cfin, cfout;

#if !defined(EIGEN_GPU_COMPILE_PHASE)
  int data_size = nthreads * 512;
  in.setRandom(data_size);
  out.setConstant(data_size, -1);
  cfin.setRandom(data_size);
  cfout.setConstant(data_size, -1);
#endif

  CALL_SUBTEST(run_and_compare_to_gpu(coeff_wise<Vector3f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(coeff_wise<Array44f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(static_map<Matrix3f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(transform_times_diagonal(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(transform_scaling(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(make_householder_small_tail(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(make_householder_complex_zero_tail(), nthreads, cfin, cfout));

#if !defined(EIGEN_USE_HIP)
  // FIXME
  // These subtests result in a compile failure on the HIP platform
  //
  //  eigen-upstream/Eigen/src/Core/Replicate.h:61:65: error:
  //           base class 'internal::dense_xpr_base<Replicate<Array<float, 4, 1, 0, 4, 1>, -1, -1> >::type'
  //           (aka 'ArrayBase<Eigen::Replicate<Eigen::Array<float, 4, 1, 0, 4, 1>, -1, -1> >') has protected default
  //           constructor
  CALL_SUBTEST(run_and_compare_to_gpu(replicate<Array4f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(replicate<Array33f>(), nthreads, in, out));

  // HIP does not support new/delete on device.
  CALL_SUBTEST(run_and_compare_to_gpu(alloc_new_delete<Vector3f>(), nthreads, in, out));
#endif

  CALL_SUBTEST(run_and_compare_to_gpu(redux<Array4f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(redux<Matrix3f>(), nthreads, in, out));

  CALL_SUBTEST(run_and_compare_to_gpu(prod_test<Matrix3f, Matrix3f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(prod_test<Matrix4f, Vector4f>(), nthreads, in, out));
  CALL_SUBTEST((test_compound_product<float, ColMajor>()));
  CALL_SUBTEST((test_compound_product<float, RowMajor>()));
  CALL_SUBTEST((test_compound_product<double, ColMajor>()));
  CALL_SUBTEST((test_compound_product<double, RowMajor>()));
  CALL_SUBTEST((test_compound_product<Eigen::half, ColMajor>()));
  CALL_SUBTEST((test_compound_product<Eigen::half, RowMajor>()));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_structured_product<RowMajor>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_structured_product<ColMajor>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_outer_product<RowMajor>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_outer_product<ColMajor>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_permutation_product<float>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_permutation_product<double>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(mixed_outer_product<RowMajor>(), nthreads, cfin, cfout));
  CALL_SUBTEST(run_and_compare_to_gpu(mixed_outer_product<ColMajor>(), nthreads, cfin, cfout));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_selfadjoint_diagonal_product<RowMajor>(), nthreads, cfin, cfout));
  CALL_SUBTEST(run_and_compare_to_gpu(scaled_selfadjoint_diagonal_product<ColMajor>(), nthreads, cfin, cfout));

  CALL_SUBTEST(run_and_compare_to_gpu(diagonal<Matrix3f, Vector3f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(diagonal<Matrix4f, Vector4f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(structured_copies(), nthreads, in, out));

  CALL_SUBTEST(run_and_compare_to_gpu(matrix_inverse<Matrix2f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(matrix_inverse<Matrix3f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(matrix_inverse<Matrix4f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(unit_orthogonal(), nthreads, in, out));

  CALL_SUBTEST((test_jacobi_rotations<float, ColMajor>()));
  CALL_SUBTEST((test_jacobi_rotations<float, RowMajor>()));
  CALL_SUBTEST((test_jacobi_rotations<double, ColMajor>()));
  CALL_SUBTEST((test_jacobi_rotations<double, RowMajor>()));

  CALL_SUBTEST(run_and_compare_to_gpu(eigenvalues_direct<Matrix3f>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(eigenvalues_direct<Matrix2f>(), nthreads, in, out));

  // Test std::complex.
  CALL_SUBTEST(run_and_compare_to_gpu(complex_operators<Vector3cf>(), nthreads, cfin, cfout));
  CALL_SUBTEST(run_and_compare_to_gpu(complex_internal_operators<std::complex<float>>(), nthreads, cfin, cfout));
  CALL_SUBTEST(test_with_infs_nans(complex_sqrt<Vector3cf>(), nthreads, cfin, cfout));

  // numeric_limits
  CALL_SUBTEST(test_with_infs_nans(numeric_limits_test<Vector3f>(), 1, in, out));

  // Eigen::half <-> __half raw-bit interop on the host.
#if defined(EIGEN_HAS_GPU_FP16) && !defined(EIGEN_GPU_COMPILE_PHASE)
  CALL_SUBTEST(test_half_raw_bit_interop());
#endif

  CALL_SUBTEST(test_custom_less_scalar_minmax());
  CALL_SUBTEST(test_float_nan_minmax());
  CALL_SUBTEST(test_bit_count());
  CALL_SUBTEST(test_half_nan_not_equal());
  CALL_SUBTEST(test_complex_exp_overflow());

  // `preverse` on every GPU packet type. 32 elements cover several packets of
  // each, including `Packet4h2`, the widest at 8.
  CALL_SUBTEST((test_reverse<Eigen::Array<float, 32, 1>>()));
  CALL_SUBTEST((test_reverse<Eigen::Array<double, 32, 1>>()));
  CALL_SUBTEST((test_reverse<Eigen::Array<Eigen::half, 32, 1>>()));
  CALL_SUBTEST(run_and_compare_to_gpu(packet_any_test<float>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(packet_any_test<double>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(packet_any_test<Eigen::half>(), nthreads, in, out));

  typedef Matrix<float, 6, 6> Matrix6f;
  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_rank2_update<Matrix4f, Lower>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_rank2_update<Matrix4f, Upper>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_rank2_update<Matrix6f, Lower>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_rank2_update<Matrix6f, Upper>(), nthreads, in, out));

  CALL_SUBTEST((run_and_compare_to_gpu(selfadjoint_permutation<Matrix4f, Lower>(), nthreads, in, out)));
  CALL_SUBTEST(
      (run_and_compare_to_gpu(selfadjoint_permutation<Matrix<float, 4, 4, RowMajor>, Upper>(), nthreads, in, out)));
  CALL_SUBTEST((run_and_compare_to_gpu(selfadjoint_permutation<Matrix4cf, Upper>(), nthreads, cfin, cfout)));
  CALL_SUBTEST((run_and_compare_to_gpu(selfadjoint_permutation<Matrix<std::complex<float>, 4, 4, RowMajor>, Lower>(),
                                       nthreads, cfin, cfout)));

  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_l1_norm<Matrix4f, Lower>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_l1_norm<Matrix4f, Upper>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_l1_norm<Matrix6f, Lower>(), nthreads, in, out));
  CALL_SUBTEST(run_and_compare_to_gpu(selfadjoint_l1_norm<Matrix6f, Upper>(), nthreads, in, out));
}