test/realview.cpp

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2025 The Eigen Authors
//
// 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

#include "main.h"

// wrapper that disables array-oriented access to real and imaginary components
struct TestComplex : public std::complex<float> {
  TestComplex() = default;
  TestComplex(std::complex<float> x) : std::complex<float>(x){};
  TestComplex(float x) : std::complex<float>(x){};
};
template <>
struct Eigen::NumTraits<TestComplex> : Eigen::NumTraits<std::complex<float>> {};
template <>
struct Eigen::internal::random_impl<TestComplex> : Eigen::internal::random_impl<std::complex<float>> {};

template <typename T>
void test_realview_readonly(const T&) {
  using Scalar = typename T::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;

  constexpr Index minRows = T::RowsAtCompileTime == Dynamic ? 1 : T::RowsAtCompileTime;
  constexpr Index maxRows = T::MaxRowsAtCompileTime == Dynamic ? (EIGEN_TEST_MAX_SIZE / 2) : T::MaxRowsAtCompileTime;
  constexpr Index minCols = T::ColsAtCompileTime == Dynamic ? 1 : T::ColsAtCompileTime;
  constexpr Index maxCols = T::MaxColsAtCompileTime == Dynamic ? (EIGEN_TEST_MAX_SIZE / 2) : T::MaxColsAtCompileTime;

  constexpr Index rowFactor = (NumTraits<Scalar>::IsComplex && !T::IsRowMajor) ? 2 : 1;
  constexpr Index colFactor = (NumTraits<Scalar>::IsComplex && T::IsRowMajor) ? 2 : 1;
  constexpr Index sizeFactor = NumTraits<Scalar>::IsComplex ? 2 : 1;

  Index rows = internal::random<Index>(minRows, maxRows);
  Index cols = internal::random<Index>(minCols, maxCols);

  T A(rows, cols), B(rows, cols);

  VERIFY(A.realView().rows() == rowFactor * A.rows());
  VERIFY(A.realView().cols() == colFactor * A.cols());
  VERIFY(A.realView().size() == sizeFactor * A.size());

  A.setRandom();
  VERIFY_IS_APPROX(A.matrix().cwiseAbs2().sum(), A.realView().matrix().cwiseAbs2().sum());
  for (Index col = 0; col < cols; ++col) {
    VERIFY_IS_APPROX(A.colwise().squaredNorm()(col), A.col(col).matrix().squaredNorm());
  }
  for (Index row = 0; row < rows; ++row) {
    VERIFY_IS_APPROX(A.rowwise().squaredNorm()(row), A.row(row).matrix().squaredNorm());
  }

  RealScalar alpha = internal::random(RealScalar(1), RealScalar(2));

  // B = A * alpha
  for (Index r = 0; r < rows; r++) {
    for (Index c = 0; c < cols; c++) {
      B.coeffRef(r, c) = A.coeff(r, c) * Scalar(alpha);
    }
  }
  VERIFY_IS_CWISE_APPROX(B.realView(), A.realView() * alpha);

  // B = A / alpha
  for (Index r = 0; r < rows; r++) {
    for (Index c = 0; c < cols; c++) {
      B.coeffRef(r, c) = A.coeff(r, c) / Scalar(alpha);
    }
  }
  VERIFY_IS_CWISE_APPROX(B.realView(), A.realView() / alpha);
}

template <typename T>
void test_realview(const T&) {
  using Scalar = typename T::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;

  constexpr Index minRows = T::RowsAtCompileTime == Dynamic ? 1 : T::RowsAtCompileTime;
  constexpr Index maxRows = T::MaxRowsAtCompileTime == Dynamic ? (EIGEN_TEST_MAX_SIZE / 2) : T::MaxRowsAtCompileTime;
  constexpr Index minCols = T::ColsAtCompileTime == Dynamic ? 1 : T::ColsAtCompileTime;
  constexpr Index maxCols = T::MaxColsAtCompileTime == Dynamic ? (EIGEN_TEST_MAX_SIZE / 2) : T::MaxColsAtCompileTime;

  constexpr Index rowFactor = (NumTraits<Scalar>::IsComplex && !T::IsRowMajor) ? 2 : 1;
  constexpr Index colFactor = (NumTraits<Scalar>::IsComplex && T::IsRowMajor) ? 2 : 1;
  constexpr Index sizeFactor = NumTraits<Scalar>::IsComplex ? 2 : 1;

  const Index rows = internal::random<Index>(minRows, maxRows);
  const Index cols = internal::random<Index>(minCols, maxCols);
  const Index realViewRows = rowFactor * rows;
  const Index realViewCols = colFactor * cols;

  const T A = T::Random(rows, cols);
  T B;

  VERIFY_IS_EQUAL(A.realView().rows(), rowFactor * A.rows());
  VERIFY_IS_EQUAL(A.realView().cols(), colFactor * A.cols());
  VERIFY_IS_EQUAL(A.realView().size(), sizeFactor * A.size());

  VERIFY_IS_APPROX(A.matrix().cwiseAbs2().sum(), A.realView().matrix().cwiseAbs2().sum());

  // test re-sizing realView during assignment
  B.realView() = A.realView();
  VERIFY_IS_APPROX(A, B);
  VERIFY_IS_APPROX(A.realView(), B.realView());

  const RealScalar alpha = internal::random(RealScalar(1), RealScalar(2));

  // B = A * alpha
  for (Index r = 0; r < rows; r++) {
    for (Index c = 0; c < cols; c++) {
      B.coeffRef(r, c) = A.coeff(r, c) * Scalar(alpha);
    }
  }
  VERIFY_IS_APPROX(B.realView(), A.realView() * alpha);

  B = A;
  B.realView() *= alpha;
  VERIFY_IS_APPROX(B.realView(), A.realView() * alpha);

  // B = A / alpha
  for (Index r = 0; r < rows; r++) {
    for (Index c = 0; c < cols; c++) {
      B.coeffRef(r, c) = A.coeff(r, c) / Scalar(alpha);
    }
  }
  VERIFY_IS_APPROX(B.realView(), A.realView() / alpha);

  B = A;
  B.realView() /= alpha;
  VERIFY_IS_APPROX(B.realView(), A.realView() / alpha);

  // force some usual access patterns
  Index malloc_size = (rows * cols * sizeof(Scalar)) + sizeof(RealScalar);
  void* data1 = internal::aligned_malloc(malloc_size);
  void* data2 = internal::aligned_malloc(malloc_size);
  Scalar* ptr1 = static_cast<Scalar*>(static_cast<void*>(static_cast<uint8_t*>(data1) + sizeof(RealScalar)));
  Scalar* ptr2 = static_cast<Scalar*>(static_cast<void*>(static_cast<uint8_t*>(data2) + sizeof(RealScalar)));
  Map<T> C(ptr1, rows, cols), D(ptr2, rows, cols);

  C.setRandom();
  D.setRandom();
  for (Index r = 0; r < realViewRows; r++) {
    for (Index c = 0; c < realViewCols; c++) {
      C.realView().coeffRef(r, c) = D.realView().coeff(r, c);
    }
  }
  VERIFY_IS_CWISE_EQUAL(C, D);

  // swap() through the views: these maps start at an odd real offset, and a reversed view has no direct access.
  C.setRandom();
  T C0 = C, D0 = D;
  C.realView().swap(D.realView());
  VERIFY_IS_CWISE_EQUAL(C, D0);
  VERIFY_IS_CWISE_EQUAL(D, C0);
  C.reverse().realView().swap(D.realView());
  VERIFY_IS_CWISE_EQUAL(C, T(C0.reverse()));
  VERIFY_IS_CWISE_EQUAL(D, T(D0.reverse()));

  C = A;

  for (Index c = 0; c < realViewCols - 1; c++) {
    B.realView().row(0).coeffRef(realViewCols - 1 - c) = C.realView().row(0).coeff(c + 1);
  }
  D.realView().row(0).tail(realViewCols - 1) = C.realView().row(0).tail(realViewCols - 1).reverse();
  VERIFY_IS_CWISE_EQUAL(B.realView().row(0).tail(realViewCols - 1), D.realView().row(0).tail(realViewCols - 1));

  for (Index r = 0; r < realViewRows - 1; r++) {
    B.realView().col(0).coeffRef(realViewRows - 1 - r) = C.realView().col(0).coeff(r + 1);
  }
  D.realView().col(0).tail(realViewRows - 1) = C.realView().col(0).tail(realViewRows - 1).reverse();
  VERIFY_IS_CWISE_EQUAL(B.realView().col(0).tail(realViewRows - 1), D.realView().col(0).tail(realViewRows - 1));
}

template <typename ComplexScalar, bool Enable = internal::packet_traits<ComplexScalar>::Vectorizable>
struct test_edge_cases_impl {
  static void run() {
    using namespace internal;
    using RealScalar = typename NumTraits<ComplexScalar>::Real;
    using ComplexPacket = typename packet_traits<ComplexScalar>::type;
    using RealPacket = typename unpacket_traits<ComplexPacket>::as_real;
    constexpr int ComplexSize = unpacket_traits<ComplexPacket>::size;
    constexpr int RealSize = 2 * ComplexSize;
    VectorX<ComplexScalar> a_data(2 * ComplexSize);
    Map<const VectorX<RealScalar>> a_data_asreal(reinterpret_cast<const RealScalar*>(a_data.data()), 2 * a_data.size());
    VectorX<RealScalar> b_data(RealSize);

    a_data.setRandom();
    evaluator<RealView<VectorX<ComplexScalar>>> eval(a_data.realView());

    for (Index offset = 0; offset < RealSize; offset++) {
      for (Index begin = 0; offset + begin < RealSize; begin++) {
        for (Index count = 0; begin + count < RealSize; count++) {
          b_data.setRandom();
          RealPacket res = eval.template packetSegment<Unaligned, RealPacket>(offset, begin, count);
          pstoreSegment(b_data.data(), res, begin, count);
          VERIFY_IS_CWISE_EQUAL(a_data_asreal.segment(offset + begin, count), b_data.segment(begin, count));
        }
      }
    }
  }
};

template <typename ComplexScalar>
struct test_edge_cases_impl<ComplexScalar, false> {
  static void run() {}
};

// Complex sequence c(i) = (i, -i) whose scalar path records the largest index it is asked for.
template <typename ComplexScalar>
struct index_recording_complex_op {
  Index* maxIndex;
  ComplexScalar operator()(Index i) const {
    *maxIndex = numext::maxi(*maxIndex, i);
    return value(i);
  }
  template <typename Packet>
  Packet packetOp(Index i) const {
    ComplexScalar values[internal::unpacket_traits<Packet>::size];
    for (int k = 0; k < internal::unpacket_traits<Packet>::size; ++k) values[k] = value(i + k);
    return internal::ploadu<Packet>(values);
  }
  static ComplexScalar value(Index i) {
    using RealScalar = typename NumTraits<ComplexScalar>::Real;
    return ComplexScalar(RealScalar(i), RealScalar(-i));
  }
};
template <typename ComplexScalar>
struct Eigen::internal::functor_traits<index_recording_complex_op<ComplexScalar>> {
  static constexpr int Cost = 1;
  static constexpr bool PacketAccess = true;
  static constexpr bool IsRepeatable = true;
};

// A segment starting at an odd (imaginary) position must not read the coefficient past its end.
template <typename ComplexScalar>
void test_odd_segment_bounds() {
  using namespace internal;
  using RealScalar = typename NumTraits<ComplexScalar>::Real;
  using ComplexPacket = typename packet_traits<ComplexScalar>::type;
  using RealPacket = typename unpacket_traits<ComplexPacket>::as_real;
  constexpr Index RealSize = 2 * unpacket_traits<ComplexPacket>::size;
  const Index n = RealSize;
  Index maxIndex = -1;
  const auto sequence = VectorX<ComplexScalar>::NullaryExpr(n, index_recording_complex_op<ComplexScalar>{&maxIndex});
  evaluator<RealView<const std::decay_t<decltype(sequence)>>> eval(sequence.realView());
  VectorX<RealScalar> expected(2 * n), actual(RealSize);
  for (Index i = 0; i < n; ++i) expected.template segment<2>(2 * i) << RealScalar(i), RealScalar(-i);
  for (Index count = 1; count <= RealSize; ++count) {
    const Index offset = 2 * n - count;
    RealPacket res = eval.template packetSegment<Unaligned, RealPacket>(offset, 0, count);
    pstoreSegment(actual.data(), res, 0, count);
    VERIFY_IS_CWISE_EQUAL(actual.head(count), expected.segment(offset, count));
    VERIFY(maxIndex < n);
  }
}

template <typename ComplexScalar, bool Enable = internal::packet_traits<ComplexScalar>::Vectorizable>
struct test_odd_segment_bounds_impl {
  static void run() { test_odd_segment_bounds<ComplexScalar>(); }
};
template <typename ComplexScalar>
struct test_odd_segment_bounds_impl<ComplexScalar, false> {
  static void run() {}
};

template <typename ComplexScalar>
void test_edge_cases(const ComplexScalar&) {
  test_edge_cases_impl<ComplexScalar>::run();
  test_odd_segment_bounds_impl<ComplexScalar>::run();
}

// The parts of non-adjacent complex coefficients have no uniform stride, so their RealView has no direct access.
template <typename ComplexScalar>
void test_strided_realview(const ComplexScalar&) {
  using RealScalar = typename NumTraits<ComplexScalar>::Real;
  using ComplexMatrix = Matrix<ComplexScalar, Dynamic, Dynamic>;
  using RealRowVector = Matrix<RealScalar, 1, Dynamic>;
  const ComplexMatrix a = ComplexMatrix::Random(6, 8);
  VERIFY((int(internal::traits<std::decay_t<decltype(a.row(2).realView())>>::Flags) & DirectAccessBit) == 0);
  VERIFY((int(internal::traits<std::decay_t<decltype(a.col(2).realView())>>::Flags) & DirectAccessBit) != 0);

  RealRowVector expected(2 * a.cols());
  for (Index j = 0; j < a.cols(); ++j) {
    expected(2 * j) = numext::real(a(2, j));
    expected(2 * j + 1) = numext::imag(a(2, j));
  }
  const Ref<const RealRowVector, 0, InnerStride<>> row = a.row(2).realView();
  VERIFY_IS_CWISE_EQUAL(row, expected);

  // A writable view keeps LvalueBit, so its coeff() and wrapping evaluators return references into the complex
  // storage.
  ComplexMatrix b = a;
  auto bRow = b.row(2);
  auto bView = bRow.realView();
  const auto& constView = bView;
  for (Index j = 0; j < constView.size(); ++j) VERIFY_IS_EQUAL(constView.coeff(j), expected(j));
  VERIFY_IS_CWISE_EQUAL(RealRowVector(b.row(2).realView().array().matrix()), expected);
  VERIFY_IS_EQUAL(b.row(2).realView().segment(1, 5).sum(), expected.segment(1, 5).sum());

  // Products read a direct-access vector operand through data() and innerStride().
  using RealMatrix = Matrix<RealScalar, Dynamic, Dynamic>;
  const RealMatrix e = RealMatrix::Random(3, 2 * a.cols());
  VERIFY_IS_APPROX(e * a.row(2).realView().transpose(), e * expected.transpose());

  // data() keeps the constness of the viewed storage, including through an IndexedView of a const matrix.
  const RealScalar* aCol1 = reinterpret_cast<const RealScalar*>(a.data() + a.rows());
  VERIFY_IS_EQUAL(a(placeholders::all, seq(1, 7, 3)).realView().data(), aCol1);
  auto aCols = a(placeholders::all, seq(1, 7, 3));
  const auto aColsView = aCols.realView();
  VERIFY_IS_EQUAL(aColsView.data(), aCol1);
}

template <typename Scalar, int Rows, int Cols, int MaxRows = Rows, int MaxCols = Cols>
void test_realview_readonly() {
  // if Rows == 1, don't test ColMajor as it is not a valid array
  using ColMajorMatrixType = Matrix<Scalar, Rows, Cols, Rows == 1 ? RowMajor : ColMajor, MaxRows, MaxCols>;
  using ColMajorArrayType = Array<Scalar, Rows, Cols, Rows == 1 ? RowMajor : ColMajor, MaxRows, MaxCols>;
  // if Cols == 1, don't test RowMajor as it is not a valid array
  using RowMajorMatrixType = Matrix<Scalar, Rows, Cols, Cols == 1 ? ColMajor : RowMajor, MaxRows, MaxCols>;
  using RowMajorArrayType = Array<Scalar, Rows, Cols, Cols == 1 ? ColMajor : RowMajor, MaxRows, MaxCols>;
  test_realview_readonly(ColMajorMatrixType());
  test_realview_readonly(ColMajorArrayType());
  test_realview_readonly(RowMajorMatrixType());
  test_realview_readonly(RowMajorArrayType());
}

template <typename Scalar, int Rows, int Cols, int MaxRows = Rows, int MaxCols = Cols>
void test_realview_readwrite() {
  // if Rows == 1, don't test ColMajor as it is not a valid array
  using ColMajorMatrixType = Matrix<Scalar, Rows, Cols, Rows == 1 ? RowMajor : ColMajor, MaxRows, MaxCols>;
  using ColMajorArrayType = Array<Scalar, Rows, Cols, Rows == 1 ? RowMajor : ColMajor, MaxRows, MaxCols>;
  // if Cols == 1, don't test RowMajor as it is not a valid array
  using RowMajorMatrixType = Matrix<Scalar, Rows, Cols, Cols == 1 ? ColMajor : RowMajor, MaxRows, MaxCols>;
  using RowMajorArrayType = Array<Scalar, Rows, Cols, Cols == 1 ? ColMajor : RowMajor, MaxRows, MaxCols>;
  test_realview(ColMajorMatrixType());
  test_realview(ColMajorArrayType());
  test_realview(RowMajorMatrixType());
  test_realview(RowMajorArrayType());
}

template <int Rows, int Cols, int MaxRows = Rows, int MaxCols = Cols>
void test_realview() {
  test_realview_readwrite<float, Rows, Cols, MaxRows, MaxCols>();
  test_realview_readwrite<std::complex<float>, Rows, Cols, MaxRows, MaxCols>();
  test_realview_readwrite<double, Rows, Cols, MaxRows, MaxCols>();
  test_realview_readwrite<std::complex<double>, Rows, Cols, MaxRows, MaxCols>();
  test_realview_readonly<TestComplex, Rows, Cols, MaxRows, MaxCols>();
}

// packet_traits<long double>::size is 1 on every supported target, so RealView's
// vectorized paths are unreachable for it and repeating it over the shape matrix
// only re-instantiates the same scalar fallback. Cover it on the fully dynamic
// and fully fixed shapes, which are the two storage layouts it can distinguish.
template <int Rows, int Cols, int MaxRows = Rows, int MaxCols = Cols>
void test_realview_long_double() {
  test_realview_readwrite<long double, Rows, Cols, MaxRows, MaxCols>();
  test_realview_readwrite<std::complex<long double>, Rows, Cols, MaxRows, MaxCols>();
}

EIGEN_DECLARE_TEST(realview) {
  for (int i = 0; i < g_repeat; i++) {
    CALL_SUBTEST_1((test_realview<Dynamic, Dynamic, Dynamic, Dynamic>()));
    CALL_SUBTEST_2((test_realview<Dynamic, Dynamic, 17, Dynamic>()));
    CALL_SUBTEST_3((test_realview<Dynamic, Dynamic, Dynamic, 19>()));
    CALL_SUBTEST_4((test_realview<Dynamic, Dynamic, 17, 19>()));
    CALL_SUBTEST_5((test_realview<17, Dynamic, 17, Dynamic>()));
    CALL_SUBTEST_6((test_realview<Dynamic, 19, Dynamic, 19>()));
    CALL_SUBTEST_7((test_realview<17, 19, 17, 19>()));
    CALL_SUBTEST_8((test_realview<Dynamic, 1>()));
    CALL_SUBTEST_9((test_realview<1, Dynamic>()));
    CALL_SUBTEST_10((test_realview<1, 1>()));
    CALL_SUBTEST_11(test_edge_cases(std::complex<float>()));
    CALL_SUBTEST_11(test_strided_realview(std::complex<float>()));
    CALL_SUBTEST_12(test_edge_cases(std::complex<double>()));
    CALL_SUBTEST_12(test_strided_realview(std::complex<double>()));
    CALL_SUBTEST_13((test_realview_long_double<Dynamic, Dynamic, Dynamic, Dynamic>()));
    CALL_SUBTEST_13((test_realview_long_double<17, 19, 17, 19>()));
  }
}