StructuredMatrices: Add Cauchy operator and GKO pivoted LU solver libeigen/eigen!2692 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com> Co-authored-by: Rasmus Munk Larsen <rlarsen@nvidia.com>
diff --git a/.agents/ci-internals.md b/.agents/ci-internals.md index abb828e..d3ad029 100644 --- a/.agents/ci-internals.md +++ b/.agents/ci-internals.md
@@ -121,6 +121,11 @@ collide whenever a regeneration is pending. The failure is not only noisy: `_ko` is `WILL_FAIL`, so a build system that errors for an unrelated reason satisfies it just as well as the compile error it is supposed to assert. +The RISC-V affected tier runs the `failtest` label on an amd64 job with the original cross compiler. Its native +runtime job excludes those compile tests and the nested `buildsystem` scenarios: the runtime image has neither Ninja +nor a compiler, and the cached compiler paths name amd64 executables. The separate `test:linux:buildsystem` job covers +the nested consumers when build-system files change. + ## Clang-Tidy Compilation Database For a source in the compilation database the driver narrows that database first, through
diff --git a/ci/test.linux.gitlab-ci.yml b/ci/test.linux.gitlab-ci.yml index 9074744..65cad49 100644 --- a/ci/test.linux.gitlab-ci.yml +++ b/ci/test.linux.gitlab-ci.yml
@@ -1020,6 +1020,22 @@ test:linux:riscv64:gcc-15:default:affected: extends: [ .test:linux:riscv64:gcc-15:default, .affected:test ] needs: [ build:linux:riscv64:gcc-15:default:affected, select:tests ] + variables: + # Compile tests need the amd64 cross toolchain; nested build-system consumers run in test:linux:buildsystem. + EIGEN_CI_CTEST_EXCLUDE: (_ok|_ko)$|^buildsystem_ + rules: !reference [.rules:libeigen:affected-tests:rvv10, rules] + +test:linux:cross:riscv64:gcc-15:failtest:affected: + extends: [ .test:linux, .affected:test ] + image: ${EIGEN_CI_IMAGE_LINUX_RISCV64_SMOKETEST_BUILD} + needs: [ build:linux:riscv64:gcc-15:default:affected, select:tests ] + variables: + EIGEN_CI_TARGET_ARCH: x86_64 + EIGEN_CI_INSTALL: "" + EIGEN_CI_CTEST_ARGS: -L failtest + EIGEN_CI_TEST_CACHE: "off" + tags: + - saas-linux-medium-amd64 rules: !reference [.rules:libeigen:affected-tests:rvv10, rules] # Three vector lengths: the SVE packet code has VL-dependent fold counts and
diff --git a/cmake/EigenSmokeTestList.cmake b/cmake/EigenSmokeTestList.cmake index 650ec6d..c62b529 100644 --- a/cmake/EigenSmokeTestList.cmake +++ b/cmake/EigenSmokeTestList.cmake
@@ -140,6 +140,7 @@ stdvector_1 stdvector_overload_1 stl_iterators_1 + structured_cauchy structured_matrices_1 structured_matrices_5 structured_matrices_13
diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake index ed47a68..b64b8ae 100644 --- a/cmake/EigenTesting.cmake +++ b/cmake/EigenTesting.cmake
@@ -394,7 +394,7 @@ # cannot tell the compile error it asserts from a build system that failed for # an unrelated reason, so a race there passes vacuously. set_tests_properties(${test_target_ok} ${test_target_ko} PROPERTIES - RESOURCE_LOCK eigen_failtest_build) + RESOURCE_LOCK eigen_failtest_build LABELS failtest) endmacro() # print a summary of the different options
diff --git a/contrib/Eigen/StructuredMatrices b/contrib/Eigen/StructuredMatrices index 3cb6abb..fafcd03 100644 --- a/contrib/Eigen/StructuredMatrices +++ b/contrib/Eigen/StructuredMatrices
@@ -39,14 +39,14 @@ * O(nk) products and O(nk^2) Woodbury solves, closed under inversion; * - \c Vandermonde : a Vandermonde matrix stored as its nodes, with Horner * products and O(n^2) Björck-Pereyra primal/dual solves - * (\c BjorckPereyra). + * (\c BjorckPereyra); + * - \c Cauchy : a Cauchy matrix stored as its node vectors, solved in + * O(n^2) with genuine partial pivoting through the displacement structure + * (\c CauchyLU, Gohberg-Kailath-Olshevsky). * * The operator types derive from \c EigenBase and store only compact generators * or factors. The FFT-backed operators (\c Circulant, \c Toeplitz and \c Hankel) * also keep a precomputed DFT symbol that every product reuses. - * The operator types derive from \c EigenBase and store only their generating - * vectors; the FFT-backed operators (\c Circulant, \c Toeplitz, \c Bccb) also - * keep a precomputed DFT symbol that every product reuses. * Because they expose \c operator* returning an Eigen product expression, they * also plug directly into the matrix-free iterative solvers * (\c ConjugateGradient, \c GMRES, \c MINRES, ...) without forming the dense @@ -88,6 +88,7 @@ #include "src/StructuredMatrices/KroneckerOperator.h" #include "src/StructuredMatrices/DiagonalPlusLowRank.h" #include "src/StructuredMatrices/Vandermonde.h" +#include "src/StructuredMatrices/Cauchy.h" // IWYU pragma: end_exports #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
diff --git a/contrib/Eigen/src/StructuredMatrices/Cauchy.h b/contrib/Eigen/src/StructuredMatrices/Cauchy.h new file mode 100644 index 0000000..e5df904 --- /dev/null +++ b/contrib/Eigen/src/StructuredMatrices/Cauchy.h
@@ -0,0 +1,558 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// 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-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 +// +// References: +// [1] J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W. Stewart, "LINPACK +// Users' Guide", SIAM, 1979. determinant()'s balanced accumulation follows +// the convention of its xGEDI routines, which return determinants as a +// (fraction, exponent) pair to avoid spurious overflow/underflow. +// [2] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974. +// Scaling by a power of two is exact, the property the balanced +// accumulation and the guarded boundary-node evaluations rely on. + +#ifndef EIGEN_STRUCTURED_CAUCHY_H +#define EIGEN_STRUCTURED_CAUCHY_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +template <typename Scalar_, int Rows_ = Dynamic, int Cols_ = Dynamic> +class Cauchy; + +template <typename Scalar_> +class CauchyLU; + +namespace internal { + +template <typename Scalar_, int Rows_, int Cols_> +struct traits<Cauchy<Scalar_, Rows_, Cols_>> { + using Scalar = Scalar_; + using StorageKind = Dense; + using XprKind = MatrixXpr; + using StorageIndex = int; + static constexpr int RowsAtCompileTime = Rows_; + static constexpr int ColsAtCompileTime = Cols_; + static constexpr int MaxRowsAtCompileTime = Rows_; + static constexpr int MaxColsAtCompileTime = Cols_; + // Deliberately no NestByRefBit: transpose(), conjugate() and adjoint() return + // owning temporaries, so Product must nest the operator by value for a + // delayed-evaluated product expression to keep its left factor alive. The copy + // is O(m+n), negligible against the O(mn) product evaluation. + static constexpr unsigned int Flags = 0; +}; + +template <typename Scalar_, int Rows_, int Cols_> +struct evaluator_traits<Cauchy<Scalar_, Rows_, Cols_>> { + using Kind = IndexBased; + using Shape = StructuredShape; +}; + +template <typename Scalar_> +struct traits<CauchyLU<Scalar_>> : traits<Matrix<Scalar_, Dynamic, Dynamic>> { + using XprKind = MatrixXpr; + using StorageKind = SolverStorage; + using StorageIndex = int; + using BaseTraits = traits<Matrix<Scalar_, Dynamic, Dynamic>>; + static constexpr unsigned int Flags = BaseTraits::Flags & RowMajorBit; + static constexpr int CoeffReadCost = Dynamic; +}; + +/** \internal \returns the Cauchy coefficient \c 1 / (a - b), guarded against a + * spurious overflow of the difference (internal::structured_guarded_diff): a + * naively formed coefficient would collapse to 1/Inf = 0, where the true value + * is a representable (possibly subnormal) number. When the guard fires, one + * half is divided directly by the halved difference, so the subnormal result is + * rounded only once. Every coefficient evaluation exposed by Cauchy -- coeff(), + * dense materialization and the product kernel -- goes through this helper. + * CauchyLU and determinant() use the same guarded difference inside their + * balanced accumulations. */ +template <typename Scalar> +Scalar cauchy_reciprocal_diff(const Scalar& a, const Scalar& b) { + using RealScalar = typename NumTraits<Scalar>::Real; + int e; + const Scalar t = structured_guarded_diff(a, b, e); + return Scalar(e == 0 ? RealScalar(1) : RealScalar(0.5)) / t; +} + +/** \internal \returns a guarded node difference as a balanced mantissa, with + * its power of two in \a exponent. */ +template <typename Scalar> +Scalar cauchy_balanced_diff(const Scalar& a, const Scalar& b, Index& exponent) { + int guardedExponent; + const Scalar difference = structured_guarded_diff(a, b, guardedExponent); + exponent = guardedExponent; + return structured_balance(difference, exponent); +} + +/** \internal Evaluates a GKO Schur-complement entry from generators stored as + * mantissa/exponent pairs. The guarded node difference is balanced before the + * division, and the accumulated power of two is applied only to the final + * result, so finite entries do not overflow or underflow in intermediate + * generator arithmetic. */ +template <typename Scalar> +Scalar cauchy_scaled_entry(const Scalar& a, Index aExponent, const Scalar& b, Index bExponent, const Scalar& x, + const Scalar& y) { + // Preserve the single-rounding coefficient path for the initial generators. + if (aExponent == 0 && bExponent == 0 && a == Scalar(1) && b == Scalar(1)) return cauchy_reciprocal_diff(x, y); + Index diffExponent; + const Scalar diff = cauchy_balanced_diff(x, y, diffExponent); + Index exponent = aExponent + bExponent - diffExponent; + const Scalar value = structured_balance(Scalar((a * b) / diff), exponent); + return structured_ldexp_clamped(value, exponent); +} + +/** \internal \returns the balanced mantissa of \c (a - b) / (a - c), with its + * power of two in \a exponent. Separately balancing both guarded differences + * prevents a representable ratio from passing through Inf/Inf or zero. */ +template <typename Scalar> +Scalar cauchy_balanced_diff_ratio(const Scalar& a, const Scalar& b, const Scalar& c, Index& exponent) { + Index numeratorExponent, denominatorExponent; + const Scalar numerator = cauchy_balanced_diff(a, b, numeratorExponent); + const Scalar denominator = cauchy_balanced_diff(a, c, denominatorExponent); + exponent = numeratorExponent - denominatorExponent; + return structured_balance(Scalar(numerator / denominator), exponent); +} + +} // namespace internal + +/** \ingroup StructuredMatrices_Module + * \class Cauchy + * \brief An \c m x \c n Cauchy matrix represented by its two node vectors. + * + * A Cauchy matrix has entry \c (i,j) equal to \f$ 1/(x_i - y_j) \f$; the class + * stores only the \c m + \c n nodes. Products are evaluated directly at O(mn) + * operations -- the same cost as a dense product, but with O(m+n) storage and + * without ever forming the matrix. (Fast approximate products via multipole + * expansions are out of scope.) + * + * The class is closed under transposition: \f$ C(x,y)^T = C(-y,-x) \f$, and the + * determinant of a square Cauchy matrix has the classical closed form + * \f$ \prod_{i<j}(x_j-x_i)(y_i-y_j) \big/ \prod_{i,j}(x_i-y_j) \f$. + * + * Because \c operator* returns an Eigen product expression, a \c Cauchy also + * drops into the matrix-free iterative solvers, and it can be assigned to a + * dense matrix when an explicit representation is needed. As with any + * matrix-free operator, the iterative solvers must be instantiated with + * \c IdentityPreconditioner (e.g. + * \c GMRES<Cauchy<double>,IdentityPreconditioner>): the default preconditioners + * read individual coefficients through \c col() or \c InnerIterator, which the + * structured operators do not expose. + * + * Square systems are solved in O(n^2) by \ref CauchyLU, a partially pivoted LU + * factorization computed through the displacement structure + * (Gohberg-Kailath-Olshevsky): unlike the Toeplitz/Levinson world, Cauchy + * structure survives row permutations, which is what makes fast \em pivoted + * factorization possible. The Hilbert matrix is the Cauchy matrix with + * \c x_i = i+1, \c y_j = -j. + * + * All nodes \c x_i must differ from all nodes \c y_j (entries would be infinite + * otherwise); this is not checked. + * + * \tparam Scalar_ the scalar type, real or complex. + * \tparam Rows_ the number of rows at compile time, or \c Dynamic (the default). + * \tparam Cols_ the number of columns at compile time, or \c Dynamic (the default). + * + * \sa class CauchyLU, makeCauchy() + */ +template <typename Scalar_, int Rows_, int Cols_> +class Cauchy : public EigenBase<Cauchy<Scalar_, Rows_, Cols_>> { + public: + using Scalar = Scalar_; + using RealScalar = typename NumTraits<Scalar>::Real; + using StorageIndex = int; + using RowNodeVector = Matrix<Scalar, Rows_, 1>; + using ColNodeVector = Matrix<Scalar, Cols_, 1>; + + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(bool(RowNodeVector::NeedsToAlign || ColNodeVector::NeedsToAlign)) + + static constexpr int RowsAtCompileTime = Rows_; + static constexpr int ColsAtCompileTime = Cols_; + static constexpr int MaxRowsAtCompileTime = Rows_; + static constexpr int MaxColsAtCompileTime = Cols_; + static constexpr int SizeAtCompileTime = internal::size_at_compile_time(Rows_, Cols_); + static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime; + static constexpr bool IsRowMajor = false; + // Deliberately no IsVectorAtCompileTime: Ref<const Cauchy>'s default StrideType + // argument reads it, so its absence makes internal::is_ref_compatible SFINAE to + // false and keeps the iterative solvers on their matrix-free path. + + /** Builds a Cauchy matrix from the row nodes \a x and the column nodes \a y: + * entry \c (i,j) is \c 1/(x[i] - y[j]). */ + template <typename XDerived, typename YDerived> + Cauchy(const MatrixBase<XDerived>& x, const MatrixBase<YDerived>& y) : m_x(x), m_y(y) { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(XDerived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(YDerived) + eigen_assert(m_x.size() > 0 && m_y.size() > 0 && "Cauchy node vectors must be non-empty"); + m_boundary = computeBoundary(); + } + + EIGEN_DEVICE_FUNC Index rows() const { return m_x.size(); } + EIGEN_DEVICE_FUNC Index cols() const { return m_y.size(); } + + /** \returns the row node vector \c x. */ + const RowNodeVector& rowNodes() const { return m_x; } + /** \returns the column node vector \c y. */ + const ColNodeVector& colNodes() const { return m_y; } + + /** \returns the coefficient at row \a row and column \a col, evaluated + * through the guarded reciprocal: nodes at the overflow boundary yield the + * correctly rounded (possibly subnormal) value instead of a spurious + * 1/Inf = 0. */ + Scalar coeff(Index row, Index col) const { return internal::cauchy_reciprocal_diff(m_x.coeff(row), m_y.coeff(col)); } + + /** \returns the transpose of \c *this, itself a Cauchy operator: + * \f$ C(x,y)^T = C(-y,-x) \f$. */ + Cauchy<Scalar, Cols_, Rows_> transpose() const { return Cauchy<Scalar, Cols_, Rows_>(-m_y, -m_x); } + + /** \returns the complex conjugate of \c *this, itself a Cauchy operator (the + * one with conjugated nodes). */ + Cauchy conjugate() const { return Cauchy(m_x.conjugate(), m_y.conjugate()); } + + /** \returns the adjoint of \c *this, itself a Cauchy operator: + * \f$ C(x,y)^H = C(-\bar y, -\bar x) \f$. */ + Cauchy<Scalar, Cols_, Rows_> adjoint() const { + return Cauchy<Scalar, Cols_, Rows_>(-m_y.conjugate(), -m_x.conjugate()); + } + + /** \returns the determinant of a \b square Cauchy matrix through the classical + * closed form \f$ \prod_{i<j}(x_j-x_i)(y_i-y_j) \big/ \prod_{i,j}(x_i-y_j) \f$, + * in O(n^2) operations. All factors enter a single accumulation kept in the + * balanced form \c m * 2^e (the split fraction/exponent determinant + * convention of LINPACK's xGEDI [1]) -- every factor and the running value + * are renormalized to unit magnitude with the power of two tracked separately + * (exact frexp/ldexp rescaling [2]), numerator factors multiplied in and + * denominator factors divided out -- so no intermediate can overflow or + * underflow when the determinant itself is representable. A node difference + * that overflows even though both nodes are finite (nodes near opposite ends + * of the exponent range) is recomputed from the halved nodes -- exact, since + * such an overflow implies huge normal operands -- with the removed power of + * two entering the same exponent bookkeeping, so it too cannot push the + * accumulation to a spurious Inf. Zero factors (coincident \c x or \c y + * nodes) and genuinely non-finite factors propagate exactly. */ + Scalar determinant() const { + eigen_assert(rows() == cols() && "Cauchy::determinant requires a square matrix"); + const Index n = rows(); + Scalar det(1); + Index exponent = 0; + for (Index j = 1; j < n; ++j) + for (Index i = 0; i < j; ++i) { + Index factorExponent; + const Scalar xDiff = internal::cauchy_balanced_diff(m_x.coeff(j), m_x.coeff(i), factorExponent); + exponent += factorExponent; + det = internal::structured_balance(Scalar(det * xDiff), exponent); + const Scalar yDiff = internal::cauchy_balanced_diff(m_y.coeff(i), m_y.coeff(j), factorExponent); + exponent += factorExponent; + det = internal::structured_balance(Scalar(det * yDiff), exponent); + } + for (Index j = 0; j < n; ++j) + for (Index i = 0; i < n; ++i) { + Index denomExponent = 0; + const Scalar d = internal::cauchy_balanced_diff(m_x.coeff(i), m_y.coeff(j), denomExponent); + exponent -= denomExponent; + det = internal::structured_balance(Scalar(det / d), exponent); + } + return internal::structured_ldexp_clamped(det, exponent); + } + + /** \internal Writes the dense representation into \a dst, one vectorized + * column at a time. Operators whose nodes reach the overflow boundary (see + * computeBoundary()) instead evaluate every entry through the guarded + * reciprocal, staying consistent with coeff(). Invoked through + * \c dense = cauchy; */ + template <typename Dest> + void evalTo(Dest& dst) const { + applyAssignment(dst, internal::assign_op<typename Dest::Scalar, Scalar>()); + } + + /** \internal Computes \c dst += (*this), see evalTo(). */ + template <typename Dest> + void addTo(Dest& dst) const { + applyAssignment(dst, internal::add_assign_op<typename Dest::Scalar, Scalar>()); + } + + /** \internal Computes \c dst -= (*this), see evalTo(). */ + template <typename Dest> + void subTo(Dest& dst) const { + applyAssignment(dst, internal::sub_assign_op<typename Dest::Scalar, Scalar>()); + } + + /** \returns the product expression \c (*this) * \a v, evaluated directly at + * O(mn) operations and O(1) extra storage. The expression carries the default + * product tag, so assigning it behaves like any dense product: a temporary + * resolves aliasing between the destination and \a v, and \c .noalias() skips + * it. */ + template <typename Rhs> + Product<Cauchy, Rhs> operator*(const MatrixBase<Rhs>& v) const { + EIGEN_STATIC_ASSERT(ColsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic || + int(ColsAtCompileTime) == int(Rhs::RowsAtCompileTime), + INVALID_MATRIX_PRODUCT) + eigen_assert(v.rows() == cols() && "invalid product: dimensions do not match"); + return Product<Cauchy, Rhs>(*this, v.derived()); + } + + /** \internal Computes \c dst += alpha * (*this) * rhs. \c ProductScalar is the + * promoted scalar of the product (complex when a real operator is applied to a + * complex right-hand side); the accumulation runs in the promoted type. The + * coefficients enter through the guarded reciprocal, so the product uses + * exactly the values coeff() exposes, boundary nodes included. */ + template <typename Dest, typename Rhs, typename ProductScalar> + void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const { + const Index m = rows(), n = cols(); + eigen_assert(rhs.rows() == n && "invalid product: dimensions do not match"); + for (Index k = 0; k < rhs.cols(); ++k) + for (Index i = 0; i < m; ++i) { + const Scalar xi = m_x.coeff(i); + ProductScalar acc(0); + for (Index j = 0; j < n; ++j) acc += rhs.coeff(j, k) * internal::cauchy_reciprocal_diff(xi, m_y.coeff(j)); + dst.coeffRef(i, k) += alpha * acc; + } + } + + private: + template <typename Dest, typename Assignment> + void applyAssignment(Dest& dst, const Assignment& assignment) const { + if (m_boundary) { + for (Index j = 0; j < cols(); ++j) + for (Index i = 0; i < rows(); ++i) + assignment.assignCoeff(dst.coeffRef(i, j), internal::cauchy_reciprocal_diff(m_x.coeff(i), m_y.coeff(j))); + return; + } + for (Index j = 0; j < cols(); ++j) { + auto dstColumn = dst.col(j); + internal::call_assignment_no_alias(dstColumn, (m_x.array() - m_y.coeff(j)).inverse().matrix(), assignment); + } + } + + /** \internal Whether some node difference \c x_i - y_j could overflow on + * finite nodes (conservative componentwise bound: the largest \c |x| and + * \c |y| components sum past the largest finite value), or a node is + * non-finite. The dense materialization then takes the guarded scalar path + * instead of the vectorized column expression; for moderate nodes the bound + * guarantees the two paths are bit-identical, so the vectorized path is kept. + * The magnitudes are taken componentwise: the modulus of a finite complex + * node near the overflow threshold is not representable. */ + bool computeBoundary() const { + RealScalar mx, my; + EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) { + mx = numext::maxi(m_x.real().cwiseAbs().maxCoeff(), m_x.imag().cwiseAbs().maxCoeff()); + my = numext::maxi(m_y.real().cwiseAbs().maxCoeff(), m_y.imag().cwiseAbs().maxCoeff()); + } else { + mx = m_x.cwiseAbs().maxCoeff(); + my = m_y.cwiseAbs().maxCoeff(); + } + return !(mx + my <= (std::numeric_limits<RealScalar>::max)()); + } + + RowNodeVector m_x; + ColNodeVector m_y; + bool m_boundary; +}; + +/** \ingroup StructuredMatrices_Module + * \returns a \ref Cauchy operator with row nodes \a x and column nodes \a y. The + * compile-time dimensions of the operator are deduced from the node vectors. */ +template <typename XDerived, typename YDerived> +Cauchy<typename XDerived::Scalar, XDerived::SizeAtCompileTime, YDerived::SizeAtCompileTime> makeCauchy( + const MatrixBase<XDerived>& x, const MatrixBase<YDerived>& y) { + return Cauchy<typename XDerived::Scalar, XDerived::SizeAtCompileTime, YDerived::SizeAtCompileTime>(x, y); +} + +/** \ingroup StructuredMatrices_Module + * \class CauchyLU + * \brief Partially pivoted O(n^2) LU solver for square Cauchy systems + * (Gohberg-Kailath-Olshevsky). + * + * A Cauchy matrix satisfies the displacement equation + * \f$ D_x C - C D_y = \mathbf{1}\mathbf{1}^T \f$, and -- crucially -- this + * structure survives row permutations. The GKO algorithm exploits it to compute + * the row-pivoted factorization \c P*C = L*U in O(n^2) operations: at each step + * the current column of the Schur complement is generated from O(n) data, the + * largest entry is chosen as pivot, and the Schur complement stays Cauchy-like + * with updated generators. This gives Cauchy systems what the Levinson world + * lacks: genuine partial pivoting at fast-algorithm cost. + * + * Usage follows the usual decomposition style, including the transposed and + * adjoint solves of \c SolverBase (from the same factorization): + * \code + * CauchyLU<double> lu(C); // or lu.compute(C); + * VectorXd u = lu.solve(b); // solve C * u = b + * VectorXd v = lu.transpose().solve(b); // solve C^T * v = b + * VectorXd w = lu.adjoint().solve(b); // solve C^H * w = b + * \endcode + * + * The factors are stored densely (O(n^2) memory, like the dense LU + * decompositions). + * + * References: I. Gohberg, T. Kailath, V. Olshevsky, "Fast Gaussian elimination + * with partial pivoting for matrices with displacement structure," Math. Comp. + * 64 (1995). + * + * \tparam Scalar_ the scalar type, real or complex. + * + * \sa class Cauchy + */ +template <typename Scalar_> +class CauchyLU : public SolverBase<CauchyLU<Scalar_>> { + public: + using Base = SolverBase<CauchyLU>; + friend class SolverBase<CauchyLU>; + EIGEN_GENERIC_PUBLIC_INTERFACE(CauchyLU) + using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic>; + using DenseVector = Matrix<Scalar, Dynamic, 1>; + using IndexVector = Matrix<Index, Dynamic, 1>; + + /** Default constructor; call \ref compute before \ref solve. */ + CauchyLU() : m_isInitialized(false), m_info(InvalidInput) {} + + /** Constructs and factorizes from the square Cauchy matrix \a C. */ + template <int Rows_, int Cols_> + explicit CauchyLU(const Cauchy<Scalar, Rows_, Cols_>& C) : m_isInitialized(false), m_info(InvalidInput) { + compute(C); + } + + /** Factorizes the square Cauchy matrix \a C as \c P*C = L*U by the GKO + * recursion with partial pivoting. The Schur complement retains the form + * \f[ S_{ij}=\frac{a_i b_j}{x_i-y_j}, \qquad + * a_i\leftarrow a_i\frac{x_i-x_k}{x_i-y_k},\quad + * b_j\leftarrow b_j\frac{y_j-y_k}{y_j-x_k}. \f] + * This follows from the rank-one displacement + * \f$D_x C-C D_y=\mathbf{1}\mathbf{1}^T\f$. \sa solve */ + template <int Rows_, int Cols_> + CauchyLU& compute(const Cauchy<Scalar, Rows_, Cols_>& C) { + eigen_assert(C.rows() == C.cols() && "CauchyLU requires a square Cauchy matrix"); + const Index n = C.rows(); + DenseVector x = C.rowNodes(); + const DenseVector y = C.colNodes(); + DenseVector a = DenseVector::Ones(n); + DenseVector b = DenseVector::Ones(n); + IndexVector aExponent = IndexVector::Zero(n); + IndexVector bExponent = IndexVector::Zero(n); + m_lu.resize(n, n); + m_perm.resize(static_cast<std::size_t>(n)); + m_info = Success; + + for (Index k = 0; k < n; ++k) { + // Guarded differences preserve subnormal Schur-complement entries at the + // exponent boundary instead of creating zero pivots or NaN generators. + Index piv = k; + RealScalar best(-1); + for (Index i = k; i < n; ++i) { + m_lu(i, k) = internal::cauchy_scaled_entry(a[i], aExponent[i], b[k], bExponent[k], x[i], y[k]); + const RealScalar mag = numext::abs(m_lu(i, k)); + if ((numext::isnan)(mag) || mag > best) { + best = mag; + piv = i; + } + } + if (piv != k) { + m_lu.row(piv).head(k + 1).swap(m_lu.row(k).head(k + 1)); + std::swap(x[piv], x[k]); + std::swap(a[piv], a[k]); + std::swap(aExponent[piv], aExponent[k]); + } + m_perm[static_cast<std::size_t>(k)] = piv; + const Scalar pivot = m_lu(k, k); + if (pivot == Scalar(0) || !(numext::isfinite)(pivot)) { + m_info = NumericalIssue; + m_lu.row(k).tail(n - k - 1).setZero(); + m_lu.col(k).tail(n - k - 1).setZero(); + continue; + } + m_lu.col(k).tail(n - k - 1) /= pivot; + for (Index j = k + 1; j < n; ++j) { + m_lu(k, j) = internal::cauchy_scaled_entry(a[k], aExponent[k], b[j], bExponent[j], x[k], y[j]); + if (!(numext::isfinite)(m_lu(k, j))) m_info = NumericalIssue; + } + for (Index i = k + 1; i < n; ++i) { + Index ratioExponent; + const Scalar ratio = internal::cauchy_balanced_diff_ratio(x[i], x[k], y[k], ratioExponent); + aExponent[i] += ratioExponent; + a[i] = internal::structured_balance(Scalar(a[i] * ratio), aExponent[i]); + } + for (Index j = k + 1; j < n; ++j) { + Index ratioExponent; + const Scalar ratio = internal::cauchy_balanced_diff_ratio(y[j], y[k], x[k], ratioExponent); + bExponent[j] += ratioExponent; + b[j] = internal::structured_balance(Scalar(b[j] * ratio), bExponent[j]); + } + } + m_isInitialized = true; + return *this; + } + + Index rows() const noexcept { return m_lu.rows(); } + Index cols() const noexcept { return m_lu.cols(); } + + /** \returns \c Success, or \c NumericalIssue when the factorization encounters + * a zero or non-finite entry that prevents a usable factorization. */ + ComputationInfo info() const { + eigen_assert(m_isInitialized && "CauchyLU is not initialized."); + return m_info; + } + +#ifdef EIGEN_PARSED_BY_DOXYGEN + /** \returns the solution \c u of \c C*u = \a b, as a lazily evaluated + * expression. Supports multiple right-hand sides. The transposed and adjoint + * systems reuse the same factorization through \c transpose().solve(b) and + * \c adjoint().solve(b). \pre \ref compute has been called. */ + template <typename Rhs> + inline const Solve<CauchyLU, Rhs> solve(const MatrixBase<Rhs>& b) const; +#endif + +#ifndef EIGEN_PARSED_BY_DOXYGEN + /** \internal P*C = L*U, so C*u = b becomes L*U*u = P*b. */ + template <typename RhsType, typename DstType> + void _solve_impl(const RhsType& rhs, DstType& dst) const { + dst = rhs; + for (Index k = 0; k < rows(); ++k) { + const Index piv = m_perm[static_cast<std::size_t>(k)]; + if (piv != k) dst.row(k).swap(dst.row(piv)); + } + m_lu.template triangularView<UnitLower>().solveInPlace(dst); + m_lu.template triangularView<Upper>().solveInPlace(dst); + } + + /** \internal C^T = U^T L^T P, so C^T*w = b becomes U^T (L^T (P w)) = b: solve + * the transposed triangles, then undo the transpositions in reverse order; + * conjugated on the way in and out for the adjoint. */ + template <bool Conjugate, typename RhsType, typename DstType> + void _solve_impl_transposed(const RhsType& rhs, DstType& dst) const { + dst = rhs.template conjugateIf<Conjugate>(); + m_lu.template triangularView<Upper>().transpose().solveInPlace(dst); + m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst); + for (Index k = rows() - 1; k >= 0; --k) { + const Index piv = m_perm[static_cast<std::size_t>(k)]; + if (piv != k) dst.row(k).swap(dst.row(piv)); + } + if (Conjugate) dst = dst.conjugate().eval(); + } +#endif + + private: + DenseMatrix m_lu; + std::vector<Index> m_perm; // transposition applied at each elimination step + bool m_isInitialized; + ComputationInfo m_info; +}; + +namespace internal { + +// Single product specialization covering every product tag; see the note in +// Circulant.h. +template <typename Scalar_, int Rows_, int Cols_, typename Rhs, int ProductTag> +struct generic_product_impl<Cauchy<Scalar_, Rows_, Cols_>, Rhs, StructuredShape, DenseShape, ProductTag> + : structured_product_impl<Cauchy<Scalar_, Rows_, Cols_>, Rhs> {}; + +} // namespace internal + +} // namespace Eigen + +#endif // EIGEN_STRUCTURED_CAUCHY_H
diff --git a/contrib/benchmarks/StructuredMatrices/CMakeLists.txt b/contrib/benchmarks/StructuredMatrices/CMakeLists.txt index 26c2428..2a7e37e 100644 --- a/contrib/benchmarks/StructuredMatrices/CMakeLists.txt +++ b/contrib/benchmarks/StructuredMatrices/CMakeLists.txt
@@ -9,3 +9,4 @@ eigen_add_benchmark(bench_structured_kronecker bench_structured_kronecker.cpp) eigen_add_benchmark(bench_structured_dplr bench_structured_dplr.cpp) eigen_add_benchmark(bench_structured_vandermonde bench_structured_vandermonde.cpp) +eigen_add_benchmark(bench_structured_cauchy bench_structured_cauchy.cpp)
diff --git a/contrib/benchmarks/StructuredMatrices/bench_structured_cauchy.cpp b/contrib/benchmarks/StructuredMatrices/bench_structured_cauchy.cpp new file mode 100644 index 0000000..0bcfea9 --- /dev/null +++ b/contrib/benchmarks/StructuredMatrices/bench_structured_cauchy.cpp
@@ -0,0 +1,183 @@ +// Benchmarks for the Cauchy operator: the O(mn)-flop, O(m+n)-storage direct +// product against its dense equivalents (with and without materializing the +// matrix first), and the O(n^2) GKO pivoted solver (CauchyLU) against the +// O(n^3) dense PartialPivLU of the materialized matrix. +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Core> +#include <Eigen/LU> +#include <contrib/Eigen/StructuredMatrices> + +using namespace Eigen; + +typedef Matrix<double, Dynamic, 1> Vec; +typedef Matrix<double, Dynamic, Dynamic> Mat; + +// After reversing its rows, this Cauchy matrix has diagonal 4n and off-diagonal +// row sums below 2 H_n. It is well conditioned while still exercising row +// pivoting in CauchyLU. +static void wellConditionedNodes(Index n, Vec& x, Vec& y) { + x.resize(n); + y.resize(n); + const double delta = 1.0 / (4.0 * static_cast<double>(n)); + for (Index i = 0; i < n; ++i) { + x[i] = static_cast<double>(n - 1 - i); + y[i] = static_cast<double>(i) - delta; + } +} + +static Mat denseCauchy(const Vec& x, const Vec& y) { + Mat dense(x.size(), y.size()); + for (Index j = 0; j < y.size(); ++j) + for (Index i = 0; i < x.size(); ++i) dense(i, j) = 1.0 / (x[i] - y[j]); + return dense; +} + +static Vec benchmarkVector(Index n) { return Vec::LinSpaced(n, -1.0, 1.0); } + +static bool productIsAccurate(const Vec& product, const Vec& x, const Vec& y, const Vec& vector) { + if (!product.allFinite()) return false; + Vec reference = Vec::Zero(x.size()); + for (Index i = 0; i < x.size(); ++i) + for (Index j = 0; j < y.size(); ++j) reference[i] += vector[j] / (x[i] - y[j]); + const double tolerance = 64.0 * static_cast<double>(y.size()) * NumTraits<double>::epsilon(); + return (product - reference).norm() <= tolerance * (reference.norm() + 1.0); +} + +static bool solveIsAccurate(const Mat& matrix, const Vec& solution, const Vec& rhs) { + if (!solution.allFinite()) return false; + const Vec residual = matrix * solution - rhs; + if (!residual.allFinite()) return false; + const double scale = matrix.norm() * solution.norm() + rhs.norm(); + const double tolerance = 512.0 * static_cast<double>(matrix.rows()) * NumTraits<double>::epsilon(); + return residual.norm() <= tolerance * scale; +} + +// --- Matrix-vector product w = C * v --- +static void BM_CauchyProduct(benchmark::State& state) { + const Index n = state.range(0); + Vec x, y; + wellConditionedNodes(n, x, y); + Cauchy<double> C(x, y); // O(n) storage; entries generated on the fly + Vec v = benchmarkVector(n), w(n); + { + const Vec product = C * v; + if (!productIsAccurate(product, x, y, v)) { + state.SkipWithError("Cauchy product validation failed"); + return; + } + } + for (auto _ : state) { + w.noalias() = C * v; + benchmark::DoNotOptimize(w.data()); + benchmark::ClobberMemory(); + } +} +BENCHMARK(BM_CauchyProduct)->Arg(64)->Arg(256)->Arg(1024); + +static void BM_DenseCauchyProduct(benchmark::State& state) { + // The dense product once the matrix has been materialized (O(n^2) storage). + const Index n = state.range(0); + Vec x, y; + wellConditionedNodes(n, x, y); + Mat dense = denseCauchy(x, y); + Vec v = benchmarkVector(n), w(n); + { + const Vec product = dense * v; + if (!productIsAccurate(product, x, y, v)) { + state.SkipWithError("dense product validation failed"); + return; + } + } + for (auto _ : state) { + w.noalias() = dense * v; + benchmark::DoNotOptimize(w.data()); + benchmark::ClobberMemory(); + } +} +BENCHMARK(BM_DenseCauchyProduct)->Arg(64)->Arg(256)->Arg(1024); + +static void BM_DenseCauchyMaterializeProduct(benchmark::State& state) { + // The dense equivalent when the matrix is not kept around: materialize from + // the node vectors, then multiply. + const Index n = state.range(0); + Vec x, y; + wellConditionedNodes(n, x, y); + Cauchy<double> C(x, y); + Vec v = benchmarkVector(n), w(n); + { + const Mat dense = C; + const Vec product = dense * v; + if (!productIsAccurate(product, x, y, v)) { + state.SkipWithError("materialized product validation failed"); + return; + } + } + for (auto _ : state) { + Mat dense = C; + w.noalias() = dense * v; + benchmark::DoNotOptimize(w.data()); + benchmark::ClobberMemory(); + } +} +BENCHMARK(BM_DenseCauchyMaterializeProduct)->Arg(64)->Arg(256)->Arg(1024); + +// --- Square solve C * u = b from the node vectors --- +static void BM_CauchyLUSolve(benchmark::State& state) { + const Index n = state.range(0); + Vec x, y; + wellConditionedNodes(n, x, y); + Cauchy<double> C(x, y); + Vec b(n); + { + const Vec expected = benchmarkVector(n); + const Mat dense = C; + b = dense * expected; + CauchyLU<double> check(C); + if (check.info() != Success) { + state.SkipWithError("CauchyLU factorization failed"); + return; + } + const Vec solution = check.solve(b); + if (!solveIsAccurate(dense, solution, b)) { + state.SkipWithError("CauchyLU validation failed"); + return; + } + } + Vec u(n); + for (auto _ : state) { + CauchyLU<double> lu(C); // O(n^2) pivoted factorization (GKO) + u = lu.solve(b); + benchmark::DoNotOptimize(u.data()); + benchmark::ClobberMemory(); + } +} +BENCHMARK(BM_CauchyLUSolve)->Arg(64)->Arg(256)->Arg(1024); + +static void BM_DensePartialPivLUSolve(benchmark::State& state) { + const Index n = state.range(0); + Vec x, y; + wellConditionedNodes(n, x, y); + Mat dense = denseCauchy(x, y); + Vec b(n); + { + const Vec expected = benchmarkVector(n); + b = dense * expected; + PartialPivLU<Mat> check(dense); + const Vec solution = check.solve(b); + if (!solveIsAccurate(dense, solution, b)) { + state.SkipWithError("PartialPivLU validation failed"); + return; + } + } + Vec u(n); + for (auto _ : state) { + PartialPivLU<Mat> lu(dense); // O(n^3) pivoted factorization + u = lu.solve(b); + benchmark::DoNotOptimize(u.data()); + benchmark::ClobberMemory(); + } +} +BENCHMARK(BM_DensePartialPivLUSolve)->Arg(64)->Arg(256)->Arg(1024);
diff --git a/contrib/test/CMakeLists.txt b/contrib/test/CMakeLists.txt index ded2f7b..5ebd65a 100644 --- a/contrib/test/CMakeLists.txt +++ b/contrib/test/CMakeLists.txt
@@ -63,6 +63,7 @@ endif() ei_add_test(NNLS) +ei_add_test(structured_cauchy) ei_add_test(sparse_extra "" "")
diff --git a/contrib/test/structured_cauchy.cpp b/contrib/test/structured_cauchy.cpp new file mode 100644 index 0000000..68f7685 --- /dev/null +++ b/contrib/test/structured_cauchy.cpp
@@ -0,0 +1,816 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// 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-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include "main.h" + +#include <contrib/Eigen/StructuredMatrices> + +using namespace Eigen; + +// Reference dense Cauchy built entry-wise from the node vectors. +template <typename Scalar> +Matrix<Scalar, Dynamic, Dynamic> reference_cauchy(const Matrix<Scalar, Dynamic, 1>& x, + const Matrix<Scalar, Dynamic, 1>& y) { + Matrix<Scalar, Dynamic, Dynamic> dense(x.size(), y.size()); + for (Index j = 0; j < y.size(); ++j) + for (Index i = 0; i < x.size(); ++i) dense(i, j) = Scalar(1) / (x[i] - y[j]); + return dense; +} + +// Reference dense Cauchy built entry-wise through the guarded reciprocal, the +// single helper every coefficient evaluation must agree with at the overflow +// boundary. +template <typename Scalar> +Matrix<Scalar, Dynamic, Dynamic> reference_cauchy_guarded(const Matrix<Scalar, Dynamic, 1>& x, + const Matrix<Scalar, Dynamic, 1>& y) { + Matrix<Scalar, Dynamic, Dynamic> dense(x.size(), y.size()); + for (Index j = 0; j < y.size(); ++j) + for (Index i = 0; i < x.size(); ++i) dense(i, j) = internal::cauchy_reciprocal_diff(x[i], y[j]); + return dense; +} + +// Separated node sets: x in [2,3], y in [0,1], so all denominators are in [1,3]. +template <typename Scalar> +void separated_nodes(Index m, Index n, Matrix<Scalar, Dynamic, 1>& x, Matrix<Scalar, Dynamic, 1>& y) { + typedef typename NumTraits<Scalar>::Real RealScalar; + x = Matrix<Scalar, Dynamic, 1>::Random(m); + y = Matrix<Scalar, Dynamic, 1>::Random(n); + x = (x * Scalar(RealScalar(0.5))).array() + Scalar(RealScalar(2.5)); + y = (y * Scalar(RealScalar(0.5))).array() + Scalar(RealScalar(0.5)); +} + +template <typename Scalar> +void test_cauchy_product(Index m, Index n) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec x, y; + separated_nodes<Scalar>(m, n, x, y); + Cauchy<Scalar> C(x, y); + Mat dense = reference_cauchy<Scalar>(x, y); + + Mat Cd = C; + VERIFY_IS_APPROX(Cd, dense); + Mat accumulated = Mat::Random(m, n); + const Mat initial = accumulated; + accumulated += C; + VERIFY_IS_APPROX(accumulated, initial + dense); + accumulated -= C; + VERIFY_IS_APPROX(accumulated, initial); + for (Index t = 0; t < 5; ++t) { + Index i = internal::random<Index>(0, m - 1), j = internal::random<Index>(0, n - 1); + VERIFY_IS_APPROX(C.coeff(i, j), dense(i, j)); + } + + Vec v = Vec::Random(n); + VERIFY_IS_APPROX((C * v).eval(), (dense * v).eval()); + + Mat V = Mat::Random(n, 3); + VERIFY_IS_APPROX((C * V).eval(), (dense * V).eval()); + + // Accumulation form exercised by the iterative solvers. + Vec w = Vec::Random(m); + Vec w0 = w; + w.noalias() += C * v; + VERIFY_IS_APPROX(w, (w0 + dense * v).eval()); +} + +template <typename Scalar> +void test_cauchy_transpose(Index m, Index n) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec x, y; + separated_nodes<Scalar>(m, n, x, y); + Cauchy<Scalar> C(x, y); + Mat dense = reference_cauchy<Scalar>(x, y); + + Mat Td = C.transpose(); + VERIFY_IS_APPROX(Td, Mat(dense.transpose())); + Mat Ad = C.adjoint(); + VERIFY_IS_APPROX(Ad, Mat(dense.adjoint())); + Mat Kd = C.conjugate(); + VERIFY_IS_APPROX(Kd, Mat(dense.conjugate())); + + Vec w = Vec::Random(m); + VERIFY_IS_APPROX((C.transpose() * w).eval(), (dense.transpose() * w).eval()); + VERIFY_IS_APPROX((C.adjoint() * w).eval(), (dense.adjoint() * w).eval()); +} + +// GKO solve on separated random nodes, verified through residuals (Cauchy +// matrices are exponentially ill-conditioned, so forward-error comparisons +// between different algorithms would not be meaningful). +template <typename Scalar> +void test_cauchy_lu(Index n) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec x, y; + separated_nodes<Scalar>(n, n, x, y); + Cauchy<Scalar> C(x, y); + Mat dense = reference_cauchy<Scalar>(x, y); + + CauchyLU<Scalar> lu(C); + VERIFY(lu.info() == Success); + + // Backward-stability-style residual bound, scaled by the scalar's epsilon (the + // separated node sets bound the GKO generators by 1, so no growth term). + const RealScalar tol = RealScalar(1000) * RealScalar(n) * NumTraits<RealScalar>::epsilon(); + + Vec b = Vec::Random(n); + Vec u = lu.solve(b); + VERIFY((dense * u - b).norm() <= tol * (dense.norm() * u.norm() + b.norm())); + + // Multiple right-hand sides. + Mat B = Mat::Random(n, 3); + Mat U = lu.solve(B); + VERIFY((dense * U - B).norm() <= tol * (dense.norm() * U.norm() + B.norm())); + + // Transposed and adjoint systems reuse the same factorization. + Vec vt = lu.transpose().solve(b); + VERIFY((dense.transpose() * vt - b).norm() <= tol * (dense.norm() * vt.norm() + b.norm())); + Vec va = lu.adjoint().solve(b); + VERIFY((dense.adjoint() * va - b).norm() <= tol * (dense.norm() * va.norm() + b.norm())); +} + +// The Hilbert matrix as a Cauchy matrix (x_i = i+1, y_j = -j): deterministic +// residual bounds far below what the astronomical conditioning would allow a +// forward-error test, plus the closed-form determinant against the dense LU one. +void test_cauchy_hilbert() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + + auto residual = [](Index n) { + Vec x(n), y(n); + for (Index i = 0; i < n; ++i) { + x[i] = double(i + 1); + y[i] = -double(i); + } + Cauchy<double> C(x, y); + Mat dense = reference_cauchy<double>(x, y); + Vec b = dense * Vec::Ones(n); + CauchyLU<double> lu(C); + VERIFY(lu.info() == Success); + Vec u = lu.solve(b); + return (dense * u - b).norm() / b.norm(); + }; + // GKO partial pivoting is backward stable, so the residual stays a small + // multiple of epsilon even as the Hilbert conditioning explodes with n. + const double eps = NumTraits<double>::epsilon(); + VERIFY(residual(8) <= 5e5 * eps); // ~1e-10 + VERIFY(residual(12) <= 5e7 * eps); // ~1e-8 + + Vec x(5), y(5); + for (Index i = 0; i < 5; ++i) { + x[i] = double(i + 1); + y[i] = -double(i); + } + Cauchy<double> H5(x, y); + Mat dense = reference_cauchy<double>(x, y); + VERIFY_IS_APPROX(H5.determinant(), dense.determinant()); +} + +// Clustered row nodes make leading minors nearly singular: partial pivoting must +// keep the factorization backward stable (residual check). +void test_cauchy_lu_pivoting(Index n) { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + + Vec x(n), y(n); + for (Index i = 0; i < n; ++i) { + x[i] = 2.0 + 1e-13 * double(i * i); // tight cluster + y[i] = double(i) / double(n); // spread out + } + Cauchy<double> C(x, y); + Mat dense = reference_cauchy<double>(x, y); + CauchyLU<double> lu(C); + VERIFY(lu.info() == Success); + Vec b = Vec::Random(n); + Vec u = lu.solve(b); + const double tol = 5e7 * NumTraits<double>::epsilon(); // ~1e-8 + VERIFY((dense * u - b).norm() <= tol * (dense.norm() * u.norm() + b.norm())); +} + +// A duplicated row node makes two rows identical, hence the matrix exactly +// singular; a zero pivot must survive partial pivoting and be reported. +void test_cauchy_lu_singular() { + typedef Matrix<double, Dynamic, 1> Vec; + Vec x(5), y(5); + x << 2.0, 2.5, 2.0, 2.75, 2.25; // x[2] duplicates x[0] + y << 0.1, 0.3, 0.5, 0.7, 0.9; + Cauchy<double> C(x, y); + CauchyLU<double> lu(C); + VERIFY(lu.info() == NumericalIssue); +} + +// Generator mantissas can overflow even when every matrix entry and LU factor +// is finite. Exponent-tracked generators must preserve the cancellation against +// the node reciprocal, and still expose an exact duplicate-row singularity. +template <typename Scalar> +void test_cauchy_lu_scaled_generators() { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec x(2), y(2), b(2); + const RealScalar P = std::ldexp(RealScalar(1), 1000); + const RealScalar t = std::ldexp(RealScalar(1), -1000); + x << Scalar(0), Scalar(-P); + y << Scalar(P), Scalar(-t); + b << Scalar(1), Scalar(0.25); + Cauchy<Scalar> C(x, y); + Mat dense = C; + VERIFY(dense.allFinite()); + CauchyLU<Scalar> lu(C); + VERIFY(lu.info() == Success); + Vec u = lu.solve(b); + VERIFY(u.allFinite()); + const RealScalar tol = RealScalar(64) * NumTraits<RealScalar>::epsilon(); + VERIFY((dense * u - b).cwiseAbs().maxCoeff() <= tol * b.cwiseAbs().maxCoeff()); + + Vec xs(2); + xs << Scalar(0), Scalar(0); + Cauchy<Scalar> Cs(xs, y); + Mat denseSingular = Cs; + VERIFY(denseSingular.allFinite()); + CauchyLU<Scalar> singularLu(Cs); + VERIFY(singularLu.info() == NumericalIssue); + + Vec xi(1), yi(1); + xi << Scalar(0); + yi << Scalar(-std::numeric_limits<RealScalar>::denorm_min()); + Cauchy<Scalar> Ci(xi, yi); + CauchyLU<Scalar> nonfiniteLu(Ci); + VERIFY(nonfiniteLu.info() == NumericalIssue); +} + +// y = -x gives the symmetric generalized Hilbert matrix 1/(x_i + x_j): verify +// the symmetry closure of the transpose through the node identity C^T = C(-y,-x). +template <typename Scalar> +void test_cauchy_symmetric(Index n) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec x(n); + for (Index i = 0; i < n; ++i) x[i] = Scalar(1) + Scalar(i); + Vec negx = -x; + Cauchy<Scalar> C(x, negx); + Mat dense = reference_cauchy<Scalar>(x, negx); + VERIFY_IS_APPROX(dense, Mat(dense.transpose())); + Cauchy<Scalar> Ct = C.transpose(); + VERIFY_IS_EQUAL(Ct.rowNodes(), Vec(x)); // -(-x) round-trips exactly + Mat Ctd = Ct; + VERIFY_IS_APPROX(Ctd, dense); +} + +// The closed form is a product of exact node differences (relative error +// O(n^2 eps)); the dense LU determinant it is compared against carries a +// cond(C)*eps error, so keep n small and the tolerance loose: the reference is +// the less accurate side of this comparison. +template <typename Scalar> +void test_cauchy_determinant(Index n) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec x, y; + separated_nodes<Scalar>(n, n, x, y); + Cauchy<Scalar> C(x, y); + Mat dense = reference_cauchy<Scalar>(x, y); + // The closed form is accurate to O(n^2 eps); the dense LU reference is the less + // accurate side, its determinant carrying a relative error that grows like + // cond(C)*eps (empirically ~0.1*cond*eps). Scale the bound by the SVD condition + // number so the test is robust across random node draws. + JacobiSVD<Mat> svd(dense); + const RealScalar cond = svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size() - 1); + const RealScalar tol = RealScalar(100) * cond * NumTraits<RealScalar>::epsilon(); + VERIFY(numext::abs(C.determinant() - dense.determinant()) <= tol * numext::abs(dense.determinant())); +} + +// The products carry the default product tag, so assignment materializes a +// temporary exactly like a dense product: x = C * x and x += C * x must see the +// pre-assignment right-hand side. Without the temporary, x = C * x would read a +// zeroed right-hand side and x += C * x would interleave destination writes +// with right-hand-side reads. +template <typename Scalar> +void test_cauchy_aliased_product(Index n) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec x, y; + separated_nodes<Scalar>(n, n, x, y); + Cauchy<Scalar> C(x, y); + Mat dense = reference_cauchy<Scalar>(x, y); + + Vec v = Vec::Random(n); + Vec w = v; + w = C * w; + VERIFY_IS_APPROX(w, (dense * v).eval()); + + w = v; + w += C * w; + VERIFY_IS_APPROX(w, (v + dense * v).eval()); + + w = v; + w -= C * w; + VERIFY_IS_APPROX(w, (v - dense * v).eval()); + + Mat V = Mat::Random(n, 3); + Mat W = V; + W = C * W; + VERIFY_IS_APPROX(W, (dense * V).eval()); +} + +// Aliasing beyond the same-object case: the default-product temporary must also +// resolve right-hand-side expressions that reference the destination, +// overlapping views of one buffer, and rectangular self-assignment where the +// destination is resized by the assignment. +template <typename Scalar> +void test_cauchy_aliased_expression(Index n) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec xn, yn; + separated_nodes<Scalar>(n, n, xn, yn); + Cauchy<Scalar> C(xn, yn); + Mat dense = reference_cauchy<Scalar>(xn, yn); + + // Right-hand-side expression referencing the destination. + Vec v = Vec::Random(n), v0 = v; + v = C * (v + Vec::Ones(n)); + VERIFY_IS_APPROX(v, (dense * (v0 + Vec::Ones(n))).eval()); + + // Overlapping (shifted) segments of one buffer. + Vec buf = Vec::Random(n + 1); + Vec expected = dense * buf.tail(n); + buf.head(n) = C * buf.tail(n); + VERIFY_IS_APPROX(buf.head(n).eval(), expected); + + // Rectangular self-assignment: z = R * z resizes the destination from n to m, + // so the product must be captured before the destination storage is touched. + const Index m = n + 3; + Vec xr, yr; + separated_nodes<Scalar>(m, n, xr, yr); + Cauchy<Scalar> R(xr, yr); + Mat denseR = reference_cauchy<Scalar>(xr, yr); + Vec z = Vec::Random(n), z0 = z; + z = R * z; + VERIFY_IS_EQUAL(z.size(), m); + VERIFY_IS_APPROX(z, (denseR * z0).eval()); +} + +// transpose()/conjugate()/adjoint() return owning temporaries, so the product +// expression must nest the structured operand by value: a delayed-evaluated +// expression has to outlive the temporary operator it was built from. The static +// check pins the value nesting; the behavioral check would read freed memory if +// the product held a reference instead. +template <typename Scalar> +void test_cauchy_delayed_product(Index m, Index n) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + STATIC_CHECK(!std::is_reference<typename internal::ref_selector<Cauchy<Scalar>>::type>::value); + + Vec x, y; + separated_nodes<Scalar>(m, n, x, y); + Cauchy<Scalar> C(x, y); + Mat dense = reference_cauchy<Scalar>(x, y); + + Vec w = Vec::Random(m); + auto expr = C.adjoint() * w; // the adjoint temporary dies with the full expression + Vec scribble = Vec::Random(m + n); // reuses the temporary's freed heap storage + Vec u = expr; + VERIFY_IS_APPROX(u, (dense.adjoint() * w).eval()); + VERIFY_IS_EQUAL(scribble.size(), m + n); // keep the scribble alive across the evaluation +} + +// Mixed-scalar products: a real operator applied to a complex right-hand side +// (and a complex operator applied to a real one) promotes to the complex product +// scalar, so alpha and the accumulation must run in the promoted type rather than +// the operator scalar. +template <typename RealScalar> +void test_cauchy_mixed_scalar(Index m, Index n) { + typedef std::complex<RealScalar> Complex; + typedef Matrix<RealScalar, Dynamic, 1> RVec; + typedef Matrix<Complex, Dynamic, 1> CVec; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + RVec x, y; + separated_nodes<RealScalar>(m, n, x, y); + Cauchy<RealScalar> C(x, y); + CMat dense = reference_cauchy<RealScalar>(x, y).template cast<Complex>(); + + CVec v = CVec::Random(n); + CVec w = C * v; + VERIFY_IS_APPROX(w, (dense * v).eval()); + + CVec w0 = CVec::Random(m); + w = w0; + w.noalias() += C * v; + VERIFY_IS_APPROX(w, (w0 + dense * v).eval()); + + CVec xc, yc; + separated_nodes<Complex>(m, n, xc, yc); + Cauchy<Complex> Cc(xc, yc); + CMat denseC = reference_cauchy<Complex>(xc, yc); + RVec vr = RVec::Random(n); + CVec z = Cc * vr; + VERIFY_IS_APPROX(z, (denseC * vr).eval()); +} + +// Wide-dynamic-range determinants: the closed-form factors overflow or underflow +// individually while the determinant itself is representable, so the balanced +// m * 2^e accumulation must carry the exponent past the intermediate extremes. +void test_cauchy_determinant_range() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + + // Reviewer repro: x = [2s, 3s], y = [0, s] with s = 1e160. The determinant is + // -1/(12 s^2) ~ -8.35e-322, a subnormal, while the naive numerator product + // -s^2 = -1e320 already overflows. Both the closed form and the dense LU + // reference round once into the subnormal range (spacing denorm_min), so the + // comparison needs an absolute term of a few subnormal spacings on top of a + // relative term; the relative term alone would be far below denorm_min. + { + const double s = 1e160; + Vec x(2), y(2); + x << 2 * s, 3 * s; + y << 0.0, s; + Cauchy<double> C(x, y); + Mat dense = reference_cauchy<double>(x, y); + const double det = C.determinant(); + const double ref = dense.partialPivLu().determinant(); + const double kRelTol = 16 * NumTraits<double>::epsilon(); + const double kAbsTol = 4 * std::numeric_limits<double>::denorm_min(); + VERIFY((numext::isfinite)(det)); + VERIFY(numext::abs(det - ref) <= kRelTol * numext::abs(ref) + kAbsTol); + } + + // Overflow-side analogue: x = [2s, 3s, 4s], y = [0, s/2, s] with s = 1e-60. + // The determinant is -1/(3780 s^3) ~ -2.65e176, huge but representable, while + // the naive numerator product -s^6/2 = -5e-361 underflows to zero. The dense + // LU reference is the less accurate side (its error carries the conditioning + // of the matrix, measured ~64 eps here); the deterministic nodes make 1000 eps + // comfortable headroom. + { + const double s = 1e-60; + Vec x(3), y(3); + x << 2 * s, 3 * s, 4 * s; + y << 0.0, s / 2, s; + Cauchy<double> C(x, y); + Mat dense = reference_cauchy<double>(x, y); + const double det = C.determinant(); + const double ref = dense.partialPivLu().determinant(); + const double kRelTol = 1000 * NumTraits<double>::epsilon(); + VERIFY((numext::isfinite)(det)); + VERIFY(numext::abs(det - ref) <= kRelTol * numext::abs(ref)); + } + + // Genuinely infinite determinant: x[0] == y[0] makes entry (0,0) infinite (a + // zero denominator factor with a non-zero numerator); the signed limit as + // x[0] -> y[0] from above is -inf, and division by the exact zero factor must + // produce it rather than a balanced-away finite value. + { + Vec x(2), y(2); + x << 1.0, 2.0; + y << 1.0, 3.0; + Cauchy<double> C(x, y); + const double det = C.determinant(); + VERIFY((numext::isinf)(det)); + VERIFY(det < 0.0); + } + + // Coincident row nodes: two identical rows, so the matrix is exactly singular + // and the zero numerator factor must propagate to an exact zero. + { + Vec x(3), y(3); + x << 2.0, 3.0, 2.0; // x[2] duplicates x[0] + y << 0.0, 0.5, 1.0; + Cauchy<double> C(x, y); + VERIFY_IS_EQUAL(C.determinant(), 0.0); + } +} + +// Node differences at the overflow boundary: forming x_j - x_i or x_i - y_j can +// overflow to Inf even though every matrix entry is finite and unexceptional. +// Such factors must enter the balanced accumulation through the exact +// halved-operand recomputation, so only the determinant's own overflow or +// underflow is visible in the result -- saturated to a zero or infinity of the +// mathematically correct sign. +void test_cauchy_determinant_overflow_boundary() { + typedef Matrix<double, Dynamic, 1> Vec; + const double M = 0.6 * (std::numeric_limits<double>::max)(); + + // Reviewer reproducer: every coefficient 1/(x_i - y_j) is finite and the exact + // determinant (2M)(-0.8M) / (0.7056 M^4) underflows, but the numerator + // difference x_1 - x_0 = 1.2 * DBL_MAX overflows if formed naively (the old + // code returned -Inf). The result must be a zero of the correct sign: one + // negative numerator factor against two negative denominator factors. + { + Vec x(2), y(2); + x << -M, M; + y << -0.4 * M, 0.4 * M; + const double det = Cauchy<double>(x, y).determinant(); + VERIFY(det == 0.0 && std::signbit(det)); + } + + // Sign flip of the same configuration: swapping the y nodes negates the + // determinant, so the underflow must land on +0. + { + Vec x(2), y(2); + x << -M, M; + y << 0.4 * M, -0.4 * M; + const double det = Cauchy<double>(x, y).determinant(); + VERIFY(det == 0.0 && !std::signbit(det)); + } + + // Purely imaginary nodes of the same magnitudes: the halved-operand + // recomputation applies componentwise to complex nodes. Scaling every node by + // i multiplies the 2x2 determinant by i^2 / i^4 = -1, so this underflow lands + // on a real part of +0 (the value is real: the four denominator divisions + // rotate the accumulation back onto the real axis). + { + typedef std::complex<double> Cplx; + Matrix<Cplx, Dynamic, 1> x(2), y(2); + x << Cplx(0.0, -M), Cplx(0.0, M); + y << Cplx(0.0, -0.4 * M), Cplx(0.0, 0.4 * M); + const Cplx det = Cauchy<Cplx>(x, y).determinant(); + VERIFY(numext::real(det) == 0.0 && numext::imag(det) == 0.0); + VERIFY(!std::signbit(numext::real(det))); + } + + // A representable determinant whose evaluation crosses the boundary, guard on + // a denominator factor (accumulated exponent decremented): with P = 2^1023, + // x = [-P, 0], y = [P, c], the difference x_0 - y_0 = -2^1024 overflows, yet + // det = (P - c) / (2 P c (P + c)), which for c = 2^-60 is 2^-964 up to a + // relative correction c/P ~ 2^-1083, far below roundoff. Every rounded factor + // is a power of two, so the balanced accumulation is exact here. + { + const double P = std::ldexp(1.0, 1023); + const double c = std::ldexp(1.0, -60); + Vec x(2), y(2); + x << -P, 0.0; + y << P, c; + const double det = Cauchy<double>(x, y).determinant(); + VERIFY_IS_APPROX(det, std::ldexp(1.0, -964)); + } + + // Guard on a numerator factor (accumulated exponent incremented): + // x = [-P, P, 0], y = [c, -c, d] with c = 2^1000 and d = 2^-1050. The + // numerator difference x_1 - x_0 = 2^1024 overflows; the determinant + // -4 P^3 (c^2 - d^2) / ((P^2 - c^2)^2 (P^2 - d^2) c d) equals + // -2^-1017 / (1 - 2^-46)^2 up to relative corrections of order 2^-2050. + { + const double P = std::ldexp(1.0, 1023); + const double c = std::ldexp(1.0, 1000); + const double d = std::ldexp(1.0, -1050); + Vec x(3), y(3); + x << -P, P, 0.0; + y << c, -c, d; + const double det = Cauchy<double>(x, y).determinant(); + const double r = 1.0 - std::ldexp(1.0, -46); // 1 - (c/P)^2 + VERIFY_IS_APPROX(det, -std::ldexp(1.0, -1017) / (r * r)); + } + + // Genuine overflow: clustered tiny nodes push the determinant past the + // representable range while every factor stays finite; the accumulated + // exponent must saturate to a correctly signed infinity, in both signs. + { + const double t = std::ldexp(1.0, -537); + Vec x(2), y(2); + x << 0.0, 3.0 * t; + y << t, 2.0 * t; // det = -3 / (4 t^2) ~ -1.6e323 + const double det = Cauchy<double>(x, y).determinant(); + VERIFY((numext::isinf)(det) && std::signbit(det)); + } + { + const double t = std::ldexp(1.0, -537); + Vec x(2), y(2); + x << 0.0, 3.0 * t; + y << 2.0 * t, t; // swapped y nodes: det = +3 / (4 t^2) + const double det = Cauchy<double>(x, y).determinant(); + VERIFY((numext::isinf)(det) && !std::signbit(det)); + } +} + +// Boundary nodes in every coefficient evaluation: when x_i - y_j overflows for +// finite nodes, a naively formed coefficient collapses to 1/Inf = 0 while the +// true value is a representable subnormal. coeff(), the dense materialization, +// the products and the CauchyLU factorization must all produce the guarded +// value -- and agree with determinant(), which derives the same quantity +// through its balanced accumulation. +void test_cauchy_boundary_coefficients() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + const double X = (std::numeric_limits<double>::max)(); + const double e = std::ldexp(1.0, -1025); // 1/(2*DBL_MAX), correctly rounded (subnormal) + + // 1x1 reviewer case: the single coefficient is 1/(2*DBL_MAX) = 2^-1025 + // (2.7813423231340017e-309). Every step is correctly rounded IEEE arithmetic + // with a power-of-two result, so the checks are exact equalities. + { + Vec x(1), y(1); + x << X; + y << -X; + Cauchy<double> C(x, y); + VERIFY_IS_EQUAL(C.coeff(0, 0), e); + Mat d = C; + VERIFY_IS_EQUAL(d(0, 0), e); + Mat acc = Mat::Zero(1, 1); + acc += C; // addTo + VERIFY_IS_EQUAL(acc(0, 0), e); + acc -= C; // subTo + VERIFY_IS_EQUAL(acc(0, 0), 0.0); + Vec p = C * Vec::Ones(1); + VERIFY_IS_EQUAL(p[0], e); + // determinant() reaches 2^-1025 through the balanced accumulation; the + // guarded coefficient matches it exactly. + VERIFY_IS_EQUAL(C.determinant(), e); + // The old code materialized a zero coefficient and reported the (regular) + // matrix as singular; the guarded pivot is subnormal but non-zero. + CauchyLU<double> lu(C); + VERIFY(lu.info() == Success); + Vec u = lu.solve(p); // p = C * [1]; recovered exactly (e / e = 1) + VERIFY_IS_EQUAL(u[0], 1.0); + } + + // A non-power-of-two boundary difference catches double rounding: computing + // (1/t)*0.5 is one subnormal ULP above the correctly rounded 0.5/t. MPFR gives + // the expected value below for x=1.1777865466541931e308 and + // y=-1.7435684109415042e308, represented by their exact double bit patterns. + { + Vec x(1), y(1); + x << numext::bit_cast<double>(numext::uint64_t(0x7fe4f71daafe5a86ull)); + y << numext::bit_cast<double>(numext::uint64_t(0xffef095b3493b755ull)); + const double expected = numext::bit_cast<double>(numext::uint64_t(0x00027621a9baa547ull)); + Cauchy<double> C(x, y); + VERIFY_IS_EQUAL(C.coeff(0, 0), expected); + Mat dense = C; + VERIFY_IS_EQUAL(dense(0, 0), expected); + Vec product = C * Vec::Ones(1); + VERIFY_IS_EQUAL(product[0], expected); + VERIFY_IS_EQUAL(C.determinant(), expected); + CauchyLU<double> lu(C); + VERIFY(lu.info() == Success); + Vec rhs(1); + rhs << expected; + Vec solution = lu.solve(rhs); + VERIFY_IS_EQUAL(solution[0], 1.0); + } + + // 2x2 with boundary pairs among moderate ones: three of the four differences + // overflow (2X and two 1.5X), one stays finite (X). All APIs agree with the + // reference built from the guarded reciprocal; scaled by its magnitude the + // matrix is well conditioned (cond ~ 38), so the GKO solve -- whose column + // generation and generator updates cross the boundary too -- recovers the + // solution accurately. + { + Vec x(2), y(2); + x << X, X / 2; + y << -X, -X / 2; + Cauchy<double> C(x, y); + Mat ref = reference_cauchy_guarded<double>(x, y); + VERIFY_IS_EQUAL(ref(0, 0), e); // 1/(2*DBL_MAX) + VERIFY_IS_EQUAL(ref(1, 1), std::ldexp(1.0, -1024)); // 1/DBL_MAX + VERIFY(ref.allFinite()); + VERIFY((ref.array() != 0.0).all()); + for (Index j = 0; j < 2; ++j) + for (Index i = 0; i < 2; ++i) VERIFY_IS_EQUAL(C.coeff(i, j), ref(i, j)); + Mat d = C; + for (Index j = 0; j < 2; ++j) + for (Index i = 0; i < 2; ++i) VERIFY_IS_EQUAL(d(i, j), ref(i, j)); + Vec v(2); + v << 0.75, -0.5; + VERIFY_IS_APPROX((C * v).eval(), (ref * v).eval()); + CauchyLU<double> lu(C); + VERIFY(lu.info() == Success); + Vec b = C * Vec::Ones(2); + Vec u = lu.solve(b); + VERIFY_IS_APPROX(u, Vec::Ones(2).eval()); + } + + // Complex nodes: purely imaginary boundary nodes overflow in the imaginary + // component of the difference; the guard applies componentwise. The paths + // sharing the helper agree exactly; the value is -i/(2*DBL_MAX) and the + // determinant derives it independently, both to within a few subnormal + // spacings (complex division rounds per component). + { + typedef std::complex<double> Cplx; + const double tiny = 4.0 * std::numeric_limits<double>::denorm_min(); + Matrix<Cplx, Dynamic, 1> xc(1), yc(1); + xc << Cplx(0.0, X); + yc << Cplx(0.0, -X); + Cauchy<Cplx> C(xc, yc); + const Cplx cval = C.coeff(0, 0); + Matrix<Cplx, Dynamic, Dynamic> dc = C; + VERIFY_IS_EQUAL(dc(0, 0), cval); + Matrix<Cplx, Dynamic, 1> pc = C * Matrix<Cplx, Dynamic, 1>::Ones(1); + VERIFY_IS_EQUAL(pc[0], cval); + VERIFY(numext::abs(numext::real(cval)) <= tiny); + VERIFY(numext::abs(numext::imag(cval) + e) <= tiny); + const Cplx det = C.determinant(); + VERIFY(numext::abs(numext::real(det) - numext::real(cval)) <= tiny); + VERIFY(numext::abs(numext::imag(det) - numext::imag(cval)) <= tiny); + } +} + +template <typename Scalar, int M, int N> +void test_cauchy_fixed() { + typedef Matrix<Scalar, M, 1> XVec; + typedef Matrix<Scalar, N, 1> YVec; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, M, N> MatMN; + + Vec xd, yd; + separated_nodes<Scalar>(M, N, xd, yd); + XVec x = xd; + YVec y = yd; + typedef Cauchy<Scalar, M, N> CauchyType; + CauchyType C(x, y); + STATIC_CHECK((Cauchy<Scalar, M, N>::RowsAtCompileTime == M)); + STATIC_CHECK((Cauchy<Scalar, M, N>::ColsAtCompileTime == N)); + STATIC_CHECK((internal::remove_all_t<decltype(makeCauchy(x, y))>::ColsAtCompileTime == N)); + STATIC_CHECK((internal::remove_all_t<decltype(C.transpose())>::RowsAtCompileTime == N)); + STATIC_CHECK((internal::remove_all_t<decltype(C.transpose())>::ColsAtCompileTime == M)); + + MatMN dense = C; + VERIFY_IS_APPROX(dense, MatMN(reference_cauchy<Scalar>(xd, yd))); + + CauchyType* heapC = new CauchyType(x, y); + VERIFY(std::uintptr_t(heapC) % std::alignment_of<CauchyType>::value == 0); + MatMN heapDense = *heapC; + delete heapC; + VERIFY_IS_APPROX(heapDense, dense); + +#if EIGEN_MAX_ALIGN_BYTES > 0 && !EIGEN_HAS_CXX17_OVERALIGN + void* raw = (CauchyType::operator new)(sizeof(CauchyType)); + VERIFY(std::uintptr_t(raw) % std::alignment_of<CauchyType>::value == 0); + (CauchyType::operator delete)(raw); +#endif + + YVec v = YVec::Random(); + Matrix<Scalar, M, 1> w = C * v; + VERIFY_IS_APPROX(w, (dense * v).eval()); + + // .noalias() keeps the direct (temporary-free) path of the default product tag; + // the matching fixed dimensions also pin the compile-time product check. + Matrix<Scalar, M, 1> w2; + w2.noalias() = C * v; + VERIFY_IS_APPROX(w2, (dense * v).eval()); +} + +EIGEN_DECLARE_TEST(structured_cauchy) { + for (int i = 0; i < g_repeat; ++i) { + // Products, dense assignment, coefficient access. + CALL_SUBTEST_1((test_cauchy_product<double>(1, 1))); + CALL_SUBTEST_1((test_cauchy_product<double>(8, 8))); + CALL_SUBTEST_1((test_cauchy_product<double>(20, 12))); // tall + CALL_SUBTEST_1((test_cauchy_product<double>(12, 20))); // wide + CALL_SUBTEST_1((test_cauchy_product<float>(10, 10))); + CALL_SUBTEST_1((test_cauchy_product<std::complex<double>>(9, 7))); + CALL_SUBTEST_1((test_cauchy_product<std::complex<float>>(7, 9))); + CALL_SUBTEST_1((test_cauchy_transpose<double>(10, 14))); + CALL_SUBTEST_1((test_cauchy_transpose<std::complex<double>>(8, 6))); + + // GKO pivoted LU solves. + CALL_SUBTEST_2((test_cauchy_lu<double>(1))); + CALL_SUBTEST_2((test_cauchy_lu<double>(2))); + CALL_SUBTEST_2((test_cauchy_lu<double>(12))); + CALL_SUBTEST_2((test_cauchy_lu<double>(30))); + CALL_SUBTEST_2((test_cauchy_lu<std::complex<double>>(16))); + CALL_SUBTEST_2((test_cauchy_lu<float>(10))); + CALL_SUBTEST_2(test_cauchy_hilbert()); + CALL_SUBTEST_2(test_cauchy_lu_pivoting(20)); + CALL_SUBTEST_2(test_cauchy_lu_singular()); + CALL_SUBTEST_2((test_cauchy_lu_scaled_generators<double>())); + CALL_SUBTEST_2((test_cauchy_lu_scaled_generators<std::complex<double>>())); + + // Closed-form determinant, symmetric generalized Hilbert, fixed sizes. + CALL_SUBTEST_3((test_cauchy_determinant<double>(4))); + CALL_SUBTEST_3((test_cauchy_determinant<std::complex<double>>(4))); + CALL_SUBTEST_3((test_cauchy_symmetric<double>(9))); + CALL_SUBTEST_3((test_cauchy_fixed<double, 6, 4>())); + CALL_SUBTEST_3((test_cauchy_fixed<std::complex<float>, 4, 5>())); + + // Numerical and lifetime boundaries: aliased and value-nested (owning) + // delayed products, mixed-scalar products, wide-dynamic-range determinants. + CALL_SUBTEST_4((test_cauchy_aliased_product<double>(11))); + CALL_SUBTEST_4((test_cauchy_aliased_product<std::complex<double>>(8))); + CALL_SUBTEST_4((test_cauchy_aliased_expression<double>(11))); + CALL_SUBTEST_4((test_cauchy_aliased_expression<std::complex<double>>(8))); + CALL_SUBTEST_4((test_cauchy_delayed_product<double>(12, 9))); + CALL_SUBTEST_4((test_cauchy_delayed_product<std::complex<double>>(7, 10))); + CALL_SUBTEST_4((test_cauchy_mixed_scalar<double>(10, 13))); + CALL_SUBTEST_4((test_cauchy_mixed_scalar<float>(9, 6))); + CALL_SUBTEST_4(test_cauchy_determinant_range()); + CALL_SUBTEST_4(test_cauchy_determinant_overflow_boundary()); + CALL_SUBTEST_4(test_cauchy_boundary_coefficients()); + } +}