StructuredMatrices: Add Bccb (block circulant with circulant blocks) libeigen/eigen!2690 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com> Co-authored-by: Rasmus Munk Larsen <rlarsen@nvidia.com>
diff --git a/contrib/Eigen/StructuredMatrices b/contrib/Eigen/StructuredMatrices index 8d3b18b..469f50c 100644 --- a/contrib/Eigen/StructuredMatrices +++ b/contrib/Eigen/StructuredMatrices
@@ -45,11 +45,14 @@ * (\c CauchyLU, Gohberg-Kailath-Olshevsky); * - \c DPR1EigenSolver : the direct O(n^2) secular-equation eigensolver for * real symmetric diagonal-plus-rank-one matrices D + rho*z*z^T, with - * LAPACK-style deflation and Gu-Eisenstat orthogonal eigenvectors. + * LAPACK-style deflation and Gu-Eisenstat orthogonal eigenvectors; + * - \c Bccb : a block circulant matrix with circulant blocks (the matrix + * of a 2-D circular convolution), diagonalized by the 2-D DFT, with + * closed-form solves, eigendecomposition and SVD. * * 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. + * or factors. The FFT-backed operators (\c Circulant, \c Toeplitz, \c Hankel and + * \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 @@ -69,9 +72,10 @@ * \c conjugate(), \c adjoint() return operators of the same kind, reusing the * cached symbol), which in particular feeds the least-squares solvers \c LSMR * and \c LeastSquaresConjugateGradient (again with \c IdentityPreconditioner). - * \c Circulant additionally exposes its closed-form eigendecomposition and SVD - * in the Fourier basis, a pseudo-inverse (minimum-norm least-squares) solve, - * \c rank(), \c inverse() and \c determinant(). + * \c Circulant and \c Bccb additionally expose their closed-form + * eigendecomposition and SVD in the Fourier basis, a pseudo-inverse + * (minimum-norm least-squares) solve, \c rank(), \c inverse() and + * \c determinant(). * * \code * #include <contrib/Eigen/StructuredMatrices> @@ -93,6 +97,7 @@ #include "src/StructuredMatrices/Vandermonde.h" #include "src/StructuredMatrices/Cauchy.h" #include "src/StructuredMatrices/DPR1EigenSolver.h" +#include "src/StructuredMatrices/Bccb.h" // IWYU pragma: end_exports #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
diff --git a/contrib/Eigen/src/StructuredMatrices/Bccb.h b/contrib/Eigen/src/StructuredMatrices/Bccb.h new file mode 100644 index 0000000..84c73f9 --- /dev/null +++ b/contrib/Eigen/src/StructuredMatrices/Bccb.h
@@ -0,0 +1,803 @@ +// 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] P. J. Davis, "Circulant Matrices", Wiley, 1979. Diagonalization of +// circulant and block circulant matrices by the DFT; the closed-form +// eigenstructure used by eigenvalues()/eigenvectors(), and the SVD and +// pseudo-inverse below, follow from it by taking moduli/phases of the +// eigenvalues. +// [2] R. H. Chan and X.-Q. Jin, "An Introduction to Iterative Toeplitz +// Solvers", SIAM, 2007. BCCB matrices are diagonalized by the 2-D DFT +// F_{n1} (x) F_{n2}; the FFT-based products and solves below, and the use +// of BCCB operators as preconditioners for two-level Toeplitz (BTTB) +// systems, follow this reference. +// [3] G. H. Golub and C. F. Van Loan, "Matrix Computations", 4th ed., Johns +// Hopkins University Press, 2013, chapter 5.4 (numerical rank conventions). +// [4] 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. +// [5] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974. +// Scaling by a power of two is exact, the property the balanced +// accumulation, the rescaled rank threshold and the scaled 2-D FFT +// products rely on. + +#ifndef EIGEN_STRUCTURED_BCCB_H +#define EIGEN_STRUCTURED_BCCB_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +template <typename Scalar_, int BlockSize_ = Dynamic, int NumBlocks_ = Dynamic> +class Bccb; + +namespace internal { + +/** \internal Compile-time product of the two circulant levels, Dynamic-aware. */ +constexpr int bccb_dim(int a, int b) { return (a == Dynamic || b == Dynamic) ? Dynamic : a * b; } + +template <typename Scalar_, int BlockSize_, int NumBlocks_> +struct traits<Bccb<Scalar_, BlockSize_, NumBlocks_>> { + using Scalar = Scalar_; + using StorageKind = Dense; + using XprKind = MatrixXpr; + using StorageIndex = int; + static constexpr int RowsAtCompileTime = bccb_dim(BlockSize_, NumBlocks_); + static constexpr int ColsAtCompileTime = RowsAtCompileTime; + static constexpr int MaxRowsAtCompileTime = RowsAtCompileTime; + static constexpr int MaxColsAtCompileTime = RowsAtCompileTime; + // 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(N), negligible against the O(N log N) product evaluation. + static constexpr unsigned int Flags = 0; +}; + +template <typename Scalar_, int BlockSize_, int NumBlocks_> +struct evaluator_traits<Bccb<Scalar_, BlockSize_, NumBlocks_>> { + using Kind = IndexBased; + using Shape = StructuredShape; +}; + +} // namespace internal + +/** \ingroup StructuredMatrices_Module + * \class Bccb + * \brief A block circulant matrix with circulant blocks (BCCB), the matrix of a + * two-dimensional circular convolution, represented by its n2 x n1 generating + * array. + * + * A BCCB matrix is the \c N x \c N matrix, \c N = n1*n2, that is circulant at two + * levels: it is an n1 x n1 block circulant whose n2 x n2 blocks are themselves + * circulant. With the generating array \c G (column \c k holds the first column + * of the k-th block), entry \c (i,j) with \c i = b1*n2 + i2, \c j = c1*n2 + j2 + * equals + * \f[ C_{i,j}=G_{(i_2-j_2)\bmod n_2,\,(b_1-c_1)\bmod n_1}. \f] + * On a column-major reshaped vector, \f$C\,\operatorname{vec}(X) + * =\operatorname{vec}(G\mathbin{\circledast}X)\f$, where \f$\circledast\f$ + * denotes 2-D circular convolution. + * + * BCCB matrices are diagonalized by the 2-D discrete Fourier transform + * \f$ F_{n_1} \otimes F_{n_2} \f$ ([1], [2]): the operator's \em symbol -- the 2-D + * DFT of \c G -- holds the eigenvalues. Products reuse that symbol or, for an + * awkward transform size, an equivalent cached padded-embedding symbol. This + * yields O(N log N) products (\c operator*), an O(N log N) direct + * (pseudo-inverse) solve (\ref solve), and closed-form factorizations: the + * eigendecomposition (\ref eigenvalues, \ref eigenvectors) and the SVD + * (\ref singularValues, \ref matrixU, \ref matrixV) in the 2-D Fourier basis, + * plus \ref rank, \ref inverse and \ref determinant. The class is closed under + * \ref transpose, \ref conjugate and \ref adjoint, which reuse the cached symbols. + * BCCB operators are the workhorse of image deblurring with periodic boundary + * conditions, and the natural preconditioners for two-level Toeplitz (BTTB) + * systems [2]. + * + * The operator stores its own copy of the generating array and derives from + * \c EigenBase. Because \c operator* returns an Eigen product expression, a + * \c Bccb 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 ConjugateGradient<Bccb<double>,Lower|Upper,IdentityPreconditioner>): + * the default preconditioners read individual coefficients through \c col() or + * \c InnerIterator, which the structured operators do not expose. + * + * Spectral operations use FFTs of the exact sizes \c n1 and \c n2. Products use + * those sizes when both are 5-smooth; otherwise an equivalent per-axis + * circulant embedding pads each awkward dimension to a 5-smooth size, avoiding + * the default kissfft backend's slow generic butterfly for large prime factors. + * + * \tparam Scalar_ the scalar type, real or complex. + * \tparam BlockSize_ the circulant block dimension \c n2 at compile time, or + * \c Dynamic (the default). + * \tparam NumBlocks_ the number of blocks \c n1 at compile time, or \c Dynamic + * (the default). + * + * \sa class Circulant, makeBccb() + */ +template <typename Scalar_, int BlockSize_, int NumBlocks_> +class Bccb : public EigenBase<Bccb<Scalar_, BlockSize_, NumBlocks_>> { + enum DenseAssignment { SetAssignment, AddAssignment, SubAssignment }; + + public: + using Scalar = Scalar_; + using RealScalar = typename NumTraits<Scalar>::Real; + using StorageIndex = int; + using Complex = std::complex<RealScalar>; + // The two-level structure keys on the *column-major* flattening of the n2 x n1 + // generating array and symbol (entry (k2, k1) is flat index k1*n2 + k2), so + // every internal 2-D type is pinned to ColMajor explicitly: the semantics must + // not change under EIGEN_DEFAULT_TO_ROW_MAJOR. The single-row special case + // only satisfies Eigen's rule that 1 x n matrices be row-major; with one row + // the two orders coincide. + using GeneratorType = + Matrix<Scalar, BlockSize_, NumBlocks_, (BlockSize_ == 1 && NumBlocks_ != 1) ? int(RowMajor) : int(ColMajor)>; + using ComplexArray = Matrix<Complex, Dynamic, Dynamic, ColMajor>; + using ComplexVector = Matrix<Complex, Dynamic, 1>; + using RealVector = Matrix<RealScalar, Dynamic, 1>; + using RealArray = Matrix<RealScalar, Dynamic, Dynamic, ColMajor>; + using ComplexMatrix = Matrix<Complex, Dynamic, Dynamic, ColMajor>; + + static constexpr int RowsAtCompileTime = internal::bccb_dim(BlockSize_, NumBlocks_); + static constexpr int ColsAtCompileTime = RowsAtCompileTime; + static constexpr int MaxRowsAtCompileTime = RowsAtCompileTime; + static constexpr int MaxColsAtCompileTime = RowsAtCompileTime; + static constexpr int SizeAtCompileTime = internal::size_at_compile_time(RowsAtCompileTime, ColsAtCompileTime); + static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime; + static constexpr bool IsRowMajor = false; + // Deliberately no IsVectorAtCompileTime: Ref<const Bccb>'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 BCCB matrix from its generating array \a generator: column \c k is + * the first column of the k-th circulant block. + * + * When the matrix is large enough for products to take the FFT path, the 2-D + * DFT of the array -- the eigenvalues of the matrix -- and any padded product + * symbol are computed here. Subsequent products and solves reuse the applicable + * cached transform. */ + template <typename Derived> + explicit Bccb(const MatrixBase<Derived>& generator) : m_g(generator) { + eigen_assert(m_g.rows() > 0 && m_g.cols() > 0 && "Bccb generator must be non-empty"); + if (m_g.size() > internal::structured_direct_threshold()) { + m_symbol = computeSymbol(); + m_prodSymbol = computeProdSymbol(m_g); + } + m_fftUsable = computeFftUsable(); + } + + EIGEN_DEVICE_FUNC Index rows() const { return m_g.size(); } + EIGEN_DEVICE_FUNC Index cols() const { return m_g.size(); } + + /** \returns the circulant block dimension \c n2. */ + Index blockSize() const { return m_g.rows(); } + /** \returns the number of blocks \c n1 in each block row. */ + Index numBlocks() const { return m_g.cols(); } + /** \returns the generating array. */ + const GeneratorType& generator() const { return m_g; } + + /** \returns the symbol of the operator: the 2-D DFT of the generating array, + * an n2 x n1 complex array whose entries are the eigenvalues of the matrix + * (see \ref eigenvalues for the ordering). Cached when the operator is large + * enough for products to take the FFT path, computed on the fly for small + * operators. */ + ComplexArray symbol() const { return m_symbol.size() > 0 ? m_symbol : computeSymbol(); } + + /** \returns the coefficient at row \a row and column \a col. */ + Scalar coeff(Index row, Index col) const { + const Index n2 = blockSize(); + Index k2 = row % n2 - col % n2; + if (k2 < 0) k2 += n2; + Index k1 = row / n2 - col / n2; + if (k1 < 0) k1 += numBlocks(); + return m_g.coeff(k2, k1); + } + + /** \returns the transpose of \c *this, itself a \c Bccb operator: the one + * generated by the array index-reversed in both dimensions. The cached + * symbols, when present, are reused -- the symbols of the transpose are their + * two-dimensional index reversals (embedding a generator commutes with + * index-reversing it, per axis, so the padded product symbol follows the same + * rule at its own grid) -- so no FFT is recomputed. */ + Bccb transpose() const { return Bccb(reverse2(m_g), reverse2(m_symbol), reverse2(m_prodSymbol)); } + + /** \returns the complex conjugate of \c *this, itself a \c Bccb operator. The + * cached symbols, when present, are reused: the symbols of the conjugate are + * the conjugated two-dimensional index reversals of the symbols. */ + Bccb conjugate() const { + return Bccb(m_g.conjugate(), reverse2(m_symbol).conjugate(), reverse2(m_prodSymbol).conjugate()); + } + + /** \returns the adjoint of \c *this, itself a \c Bccb operator. The cached + * symbols, when present, are reused: the symbols of the adjoint are the + * elementwise conjugates of the symbols (the eigenvalues conjugate while the + * 2-D Fourier eigenbasis stays fixed). */ + Bccb adjoint() const { return Bccb(reverse2(m_g).conjugate(), m_symbol.conjugate(), m_prodSymbol.conjugate()); } + + /** \returns the minimum-norm least-squares solution of \c (*this) * x = b, + * computed directly in the 2-D Fourier domain. Symbol entries whose modulus + * reaches the rank threshold (see \ref rank) are inverted; the remaining ones + * are treated as exact zeros, so the result is the pseudo-inverse applied to + * \a b. For a non-singular operator this is the exact solution. Supports + * multiple right-hand sides. */ + template <typename Rhs> + Matrix<Scalar, RowsAtCompileTime, Rhs::ColsAtCompileTime> solve(const MatrixBase<Rhs>& b) const { + EIGEN_STATIC_ASSERT(RowsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic || + int(RowsAtCompileTime) == int(Rhs::RowsAtCompileTime), + YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES) + const Index N = rows(); + eigen_assert(b.rows() == N && "right-hand side has the wrong number of rows"); + const ComplexArray s = symbol(); + RealArray mods; + RealScalar tol; + scaledModuli(s, mods, tol); + ComplexArray sinv(s.rows(), s.cols()); + // Strictly-below-threshold entries are zeroed, matching SVDBase::rank(), so a + // smallest-normal 1x1 operator stays invertible. A comparison with NaN is + // false, so NaN symbol entries stay in the inverted set and propagate to the + // output instead of being silently zeroed. + // Keep the reciprocal behind a scalar branch: select() evaluates both arms, + // so it would divide by thresholded zeros and potentially raise floating-point + // exceptions even though those coefficients are discarded. + for (Index k = 0; k < s.size(); ++k) sinv(k) = mods(k) < tol ? Complex(0) : Complex(1) / s(k); + Matrix<Scalar, RowsAtCompileTime, Rhs::ColsAtCompileTime> x(N, b.cols()); + if (!b.allFinite()) { + // A non-finite right-hand side cannot go through the transforms (see + // addProduct): apply the pseudo-inverse -- itself a BCCB operator, + // generated by the 2-D inverse DFT of the thresholded reciprocal symbol -- + // through the direct kernel so Inf/NaN propagate entrywise. + const GeneratorType pg = generatorFromSymbol(sinv); + x.setZero(); + Bccb(pg, ComplexArray(), ComplexArray()).directProduct(x, b.derived(), Scalar(1)); + return x; + } + // The right-hand side is finite here (checked above), so no per-column + // direct kernel can be needed. + applySymbol(x, sinv, b.derived(), Scalar(1), /*accumulate=*/false, + [](Index) { eigen_assert(false && "non-finite column requires a direct kernel"); }); + return x; + } + + /** \returns the numerical rank: the number of symbol entries whose modulus is + * no smaller than the threshold \c N * epsilon * max|symbol|, both evaluated + * in an exactly rescaled frame so the moduli cannot overflow (see + * scaledModuli()) and clamped from below like \c SVDBase::rank(). This is the + * same threshold \ref solve uses to decide which Fourier components to invert, + * and the comparison is strict like SVDBase's, so an entry sitting exactly on + * the threshold still counts as non-zero. */ + Index rank() const { + const ComplexArray s = symbol(); + RealArray mods; + RealScalar tol; + scaledModuli(s, mods, tol); + return (!(mods.array() < tol)).count(); // negated so NaN entries count as non-zero + } + + /** \returns the inverse of \c *this, itself a \c Bccb operator: the one + * generated by the 2-D inverse DFT of the entrywise-inverted symbol. + * \warning The operator must be non-singular; use \ref solve for a + * pseudo-inverse solve of a rank-deficient operator. */ + Bccb inverse() const { + ComplexArray sinv = symbol().cwiseInverse(); + const GeneratorType ginv = generatorFromSymbol(sinv); + const bool cache = rows() > internal::structured_direct_threshold(); + return Bccb(ginv, cache ? sinv : ComplexArray(), cache ? computeProdSymbol(ginv) : ComplexArray()); + } + + /** \returns the determinant, i.e. the product of the eigenvalues (the symbol + * entries). The product is accumulated in the balanced form \c m * 2^e (the + * split fraction/exponent determinant convention of LINPACK's xGEDI [4]) -- + * every factor and the running product are renormalized to unit magnitude with + * the power of two tracked separately -- so the partial products can neither + * overflow nor underflow when the determinant itself is representable, whatever + * the ordering of large and small eigenvalues. For a real operator the product + * is real up to roundoff, and its real part is returned. */ + Scalar determinant() const { + const ComplexArray s = symbol(); + Complex det(1); + Index exponent = 0; + for (Index k1 = 0; k1 < s.cols(); ++k1) + for (Index k2 = 0; k2 < s.rows(); ++k2) + det = internal::structured_balance(Complex(det * internal::structured_balance(s(k2, k1), exponent)), exponent); + return toScalar(internal::structured_ldexp_clamped(det, exponent), std::is_same<RealScalar, Scalar>()); + } + + /** \returns the eigenvalues as the column-major flattening of the symbol: + * eigenvalue \c f1*n2 + f2 is \c symbol()(f2, f1), and its (unit-norm) + * eigenvector is the Kronecker product of the 1-D Fourier vectors of + * frequencies \c f1 and \c f2, i.e. column \c f1*n2 + f2 of \ref eigenvectors. + * Every BCCB matrix is diagonalized by this same 2-D Fourier basis. */ + ComplexVector eigenvalues() const { return symbol().reshaped(); } + + /** \returns the unitary matrix of eigenvectors: column \c f1*n2 + f2 is the + * 2-D Fourier vector matching \c eigenvalues()[f1*n2 + f2]. + * \note The eigenvector matrix is materialized as a dense \c N x \c N matrix; + * unlike the other methods of this class this costs O(N^2) storage. */ + ComplexMatrix eigenvectors() const { + const Index n2 = blockSize(), n1 = numBlocks(), N = rows(); + ComplexMatrix F(N, N); + for (Index f1 = 0; f1 < n1; ++f1) + for (Index f2 = 0; f2 < n2; ++f2) fourierColumn(F, f1, f2, f1 * n2 + f2); + return F; + } + + /** \returns the singular values, sorted in decreasing order: the moduli of the + * symbol entries. The ordering is shared with \ref matrixU and \ref matrixV, + * so together they form the SVD \c *this = U * singularValues().asDiagonal() * V^H. */ + RealVector singularValues() const { + const ComplexVector s = eigenvalues(); + const RealVector mods = s.cwiseAbs(); + const std::vector<Index> perm = internal::structured_svd_permutation(mods); + RealVector sv(s.size()); + for (Index t = 0; t < s.size(); ++t) sv[t] = mods[perm[t]]; + return sv; + } + + /** \returns the matrix of left singular vectors \c U: column \c t is the 2-D + * Fourier vector of the t-th largest symbol entry, scaled by its phase (phase 1 + * for a zero entry). Dense \c N x \c N, see the note in \ref eigenvectors. */ + ComplexMatrix matrixU() const { + const Index n2 = blockSize(), N = rows(); + const ComplexVector s = eigenvalues(); + const RealVector mods = s.cwiseAbs(); + const std::vector<Index> perm = internal::structured_svd_permutation(mods); + ComplexMatrix U(N, N); + for (Index t = 0; t < N; ++t) { + fourierColumn(U, perm[t] / n2, perm[t] % n2, t); + const RealScalar a = mods[perm[t]]; + if (a > RealScalar(0)) U.col(t) *= s[perm[t]] / a; + } + return U; + } + + /** \returns the matrix of right singular vectors \c V: column \c t is the 2-D + * Fourier vector of the t-th largest symbol entry. Dense \c N x \c N, see the + * note in \ref eigenvectors. */ + ComplexMatrix matrixV() const { + const Index n2 = blockSize(), N = rows(); + const ComplexVector s = eigenvalues(); + const std::vector<Index> perm = internal::structured_svd_permutation(RealVector(s.cwiseAbs())); + ComplexMatrix V(N, N); + for (Index t = 0; t < N; ++t) fourierColumn(V, perm[t] / n2, perm[t] % n2, t); + return V; + } + + /** \internal Writes the dense representation into \a dst: the (b1,c1) block is + * the generator column (b1-c1) mod n1 rotated downwards by the within-block + * column index, so only contiguous segment copies are involved. Invoked + * through \c dense = bccb; */ + template <typename Dest> + void evalTo(Dest& dst) const { + assignTo<SetAssignment>(dst); + } + + /** \internal Computes \c dst += (*this), see evalTo(). */ + template <typename Dest> + void addTo(Dest& dst) const { + assignTo<AddAssignment>(dst); + } + + /** \internal Computes \c dst -= (*this), see evalTo(). */ + template <typename Dest> + void subTo(Dest& dst) const { + assignTo<SubAssignment>(dst); + } + + /** \returns the product expression \c (*this) * \a x, evaluated through a fast + * 2-D-FFT-based matrix-vector product. The expression carries the default + * product tag, so assigning it behaves like any dense product: a temporary + * resolves aliasing between the destination and \a x, and \c .noalias() skips + * it. */ + template <typename Rhs> + Product<Bccb, Rhs> operator*(const MatrixBase<Rhs>& x) const { + EIGEN_STATIC_ASSERT(ColsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic || + int(ColsAtCompileTime) == int(Rhs::RowsAtCompileTime), + INVALID_MATRIX_PRODUCT) + eigen_assert(x.rows() == cols() && "invalid product: dimensions do not match"); + return Product<Bccb, Rhs>(*this, x.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. + * + * Non-finite data takes the direct O(N^2) kernel: the transforms would smear a + * single Inf/NaN into NaNs across the whole output, where the dense product + * only propagates it through the dot products that touch it. A non-finite + * generating array or cached symbol (which can overflow even for a finite + * generator: the 2-D DFT accumulates up to N addends) routes the whole + * product; a non-finite right-hand-side column is detected inside the FFT + * loop -- in the same pass that derives its scaling exponent, so finite data + * pays no extra scan -- and falls back per column. */ + template <typename Dest, typename Rhs, typename ProductScalar> + void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const { + const Index N = rows(); + eigen_assert(rhs.rows() == N && "invalid product: dimensions do not match"); + if (N <= internal::structured_direct_threshold() || !m_fftUsable) { + directProduct(dst, rhs, alpha); + return; + } + // Products use the padded embedding symbol when a block dimension is not + // 5-smooth (see computeProdSymbol()); the exact-size transform of such a + // dimension runs through kissfft's quadratic generic butterfly. + const ComplexArray& s = m_prodSymbol.size() > 0 ? m_prodSymbol : m_symbol; + applySymbol(dst, s, rhs, alpha, /*accumulate=*/true, [&](Index k) { directProductColumn(dst, rhs, k, alpha); }); + } + + private: + /** \internal Materializes the dense BCCB representation with the selected + * assignment mode; shared by evalTo(), addTo() and subTo(). */ + template <DenseAssignment Assignment, typename Dest> + void assignTo(Dest& dst) const { + const Index n2 = blockSize(), n1 = numBlocks(); + for (Index c1 = 0; c1 < n1; ++c1) + for (Index j2 = 0; j2 < n2; ++j2) { + auto col = dst.col(c1 * n2 + j2); + for (Index b1 = 0; b1 < n1; ++b1) { + Index k1 = b1 - c1; + if (k1 < 0) k1 += n1; + EIGEN_IF_CONSTEXPR (Assignment == SetAssignment) { + col.segment(b1 * n2, j2) = m_g.col(k1).tail(j2); + col.segment(b1 * n2 + j2, n2 - j2) = m_g.col(k1).head(n2 - j2); + } else EIGEN_IF_CONSTEXPR (Assignment == AddAssignment) { + col.segment(b1 * n2, j2) += m_g.col(k1).tail(j2); + col.segment(b1 * n2 + j2, n2 - j2) += m_g.col(k1).head(n2 - j2); + } else { + col.segment(b1 * n2, j2) -= m_g.col(k1).tail(j2); + col.segment(b1 * n2 + j2, n2 - j2) -= m_g.col(k1).head(n2 - j2); + } + } + } + } + + /** \internal Builds an operator from a generating array and already-known + * symbols (empty for small operators; \a prodSymbol also empty when both block + * dimensions are 5-smooth, which needs no padded product embedding), skipping + * the FFTs of the public constructor. Used by transpose(), conjugate(), + * adjoint() and inverse(), whose symbols are cheap transformations of the + * existing ones. */ + Bccb(const GeneratorType& g, const ComplexArray& symbol, const ComplexArray& prodSymbol) + : m_g(g), m_symbol(symbol), m_prodSymbol(prodSymbol) { + m_fftUsable = computeFftUsable(); + } + + /** \internal Whether products may take the FFT path: the generating array and + * the cached symbols must be finite. The symbols accumulate up to N + * (respectively p2*p1) addends, so they can overflow to Inf even for a finite + * generator; such operators fall back to the direct kernel, which stays + * exact. */ + bool computeFftUsable() const { + return m_g.allFinite() && (m_symbol.size() == 0 || m_symbol.allFinite()) && + (m_prodSymbol.size() == 0 || m_prodSymbol.allFinite()); + } + + /** \internal \returns the symbol products should use for a generating array + * \a g: empty when both block dimensions are 5-smooth (the exact-size + * transforms are already fast), otherwise the 2-D DFT of the per-axis + * circulant embedding of \a g. kissfft falls back to a quadratic generic + * butterfly for prime factors other than 2, 3 and 5, and the row-column + * algorithm pays that cost in every transform of the awkward axis, so each + * such axis is padded to \c p_i = fft_next_good_size(2*n_i - 1) and the cyclic + * convolution along it is evaluated as a padded linear convolution -- exactly + * the Toeplitz embedding of a circulant, applied per dimension. Along a + * padded axis the generator is laid out as [g; 0...0; g(1..)] (the wrapped + * band g(1..) sits at the top indices), so that for reachable offsets + * \c d = i - j in [-(n-1), n-1] the embedded array satisfies + * \c Ge(d mod p) = g(d mod n); offsets outside that range only touch the zero + * band. The two axes decouple, so the identity holds per entry of the 2-D + * grid and the leading (n2, n1) block of the padded cyclic convolution equals + * the exact one. A 5-smooth axis is kept at its exact size with no embedding. + * The spectral operations (eigenvalues, solve, rank, determinant, SVD) keep + * the exact-size symbol, whose entries are the eigenvalues. */ + ComplexArray computeProdSymbol(const GeneratorType& g) const { + const Index n2 = g.rows(), n1 = g.cols(); + const Index p2 = internal::fft_next_good_size(n2) == n2 ? n2 : internal::fft_next_good_size(2 * n2 - 1); + const Index p1 = internal::fft_next_good_size(n1) == n1 ? n1 : internal::fft_next_good_size(2 * n1 - 1); + if (p2 == n2 && p1 == n1) return ComplexArray(); + ComplexArray Ge = ComplexArray::Zero(p2, p1); + Ge.topLeftCorner(n2, n1) = g.template cast<Complex>(); + if (p2 != n2) Ge.bottomLeftCorner(n2 - 1, n1) = g.bottomRows(n2 - 1).template cast<Complex>(); + if (p1 != n1) Ge.topRightCorner(n2, n1 - 1) = g.rightCols(n1 - 1).template cast<Complex>(); + if (p2 != n2 && p1 != n1) + Ge.bottomRightCorner(n2 - 1, n1 - 1) = g.bottomRightCorner(n2 - 1, n1 - 1).template cast<Complex>(); + fft2(Ge); + return Ge; + } + + /** \internal Direct O(N^2) kernel for column \a k of the right-hand side: + * computes \c dst.col(k) += alpha * (*this) * rhs.col(k) without transforms, + * as a plain scalar loop (the two-level segment-based middle tier of the 1-D + * operators does not pay off here: the within-block segments are too short at + * these sizes). Serves operators below the FFT threshold and any column + * involving non-finite data, whose entrywise IEEE semantics the transforms + * cannot preserve. */ + template <typename Dest, typename Rhs, typename ProductScalar> + void directProductColumn(Dest& dst, const Rhs& rhs, Index k, const ProductScalar& alpha) const { + const Index N = rows(); + // A unit alpha must not multiply: even the identity complex scalar (1,0) + // pollutes an (Inf,0) value with NaN through the 0*Inf cross term. + const bool unitAlpha = alpha == ProductScalar(1); + for (Index i = 0; i < N; ++i) { + ProductScalar acc(0); + for (Index j = 0; j < N; ++j) acc += coeff(i, j) * rhs.coeff(j, k); + dst.coeffRef(i, k) += unitAlpha ? acc : ProductScalar(alpha * acc); + } + } + + /** \internal Direct O(N^2) product kernel over every column, see + * directProductColumn(). */ + template <typename Dest, typename Rhs, typename ProductScalar> + void directProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const { + for (Index k = 0; k < rhs.cols(); ++k) directProductColumn(dst, rhs, k, alpha); + } + + /** \internal Applies the operator whose 2-D symbol is \a s to every column of + * \a rhs: reshape to n2 x n1 (column-major), 2-D FFT, multiply by \a s, + * back-transform, take the \c ProductScalar part. Adds into \a dst when + * \a accumulate, overwrites otherwise. The transforms run on the grid of + * \a s itself: for the padded product symbol (see computeProdSymbol()) the + * column is zero-padded into the (p2, p1) grid and the leading (n2, n1) block + * of the back-transform is read out; for an exact-size symbol the grid is + * (n2, n1) and no padding is involved. + * + * The transforms sum up to \c N inputs, so intermediates can overflow even + * when every entry of the true result is representable. To keep the whole + * pipeline overflow-free for finite inputs anywhere in the floating-point + * range, the symbol and each right-hand-side column are rescaled by an exact + * power of two that brings their maximum modulus below one -- huge generators + * make huge symbols, hence both sides -- and the removed exponent is folded + * back into the output through per-entry ldexp, which saturates cleanly to + * zero or infinity if the true result itself leaves the representable range. + * The exponents come from the component-wise magnitudes: the modulus of a + * finite complex value near the overflow threshold is not representable, + * which would silently disable the scaling exactly where it is needed. + * + * A single plain (fast-max) reduction pass per column yields the finiteness + * routing, the zero-column shortcut and the scaling exponent (see + * internal::structured_exponent_bound_finite() for why missing a NaN in the + * fast max cannot change a routing decision). Non-finite columns are handed + * to \a directColumn, the caller's per-column direct kernel, so the remaining + * columns keep the fast path; a non-finite symbol is the caller's + * responsibility (addProduct routes it to the direct kernel up front, and + * solve() takes its dedicated pseudo-inverse fallback for non-finite + * right-hand sides). Zero columns short-circuit to zero, but only under a + * finite symbol (Inf*0 and NaN*0 are NaN entrywise) and after an exact + * recheck: the fast max can miss a NaN hiding among zeros (an Inf always + * surfaces), and such a column must fall back, not shortcut. */ + template <typename Dest, typename Rhs, typename ProductScalar, typename DirectColumn> + void applySymbol(Dest& dst, const ComplexArray& s, const Rhs& rhs, const ProductScalar& alpha, bool accumulate, + DirectColumn&& directColumn) const { + const Index n2 = blockSize(), n1 = numBlocks(), N = rows(); + const Index p2 = s.rows(), p1 = s.cols(); // the symbol's transform grid: (n2, n1) or the padded product grid + const bool padded = p2 != n2 || p1 != n1; + const bool sFinite = s.allFinite(); + // max|s| < 2^es; component bounds avoid overflow in complex moduli. + const int es = internal::structured_exponent_bound(s); + ComplexArray sScaled; + if (es != 0) { + sScaled = s; + ldexpInPlace(sScaled, -es); + } + const ComplexArray& sUse = es != 0 ? sScaled : s; + Matrix<ProductScalar, Dynamic, 1> xc(N); + ComplexArray X(p2, p1), Xn; // Xn: the leading-block extraction of a padded grid + for (Index k = 0; k < rhs.cols(); ++k) { + xc = rhs.col(k).template cast<ProductScalar>(); + const RealScalar m = xc.realView().cwiseAbs().maxCoeff(); + if (!(numext::isfinite)(m)) { + directColumn(k); + continue; + } + if (m == RealScalar(0)) { + // The fast max cannot hide an Inf (those comparisons are ordered), but + // it can miss a NaN among zeros: recheck exactly before shortcutting. + if ((xc.array() == ProductScalar(0)).all()) { + if (sFinite) { + // An exactly zero column maps to an exactly zero column -- unless + // the symbol holds Inf or NaN, whose products with zero are NaN; + // falling through to the transforms produces exactly that. + if (!accumulate) dst.col(k).setZero(); + continue; + } + } else { + directColumn(k); + continue; + } + } + int ex = 0; // stays 0 for an all-zero column: no scaling + if (m > RealScalar(0)) { + EIGEN_USING_STD(frexp); + frexp(m, &ex); + EIGEN_IF_CONSTEXPR (NumTraits<ProductScalar>::IsComplex) ++ex; + } + // reshaped() defaults to column-major traversal, the flattening the + // two-level structure keys on, regardless of EIGEN_DEFAULT_TO_ROW_MAJOR. + if (padded) { + X.setZero(); + X.topLeftCorner(n2, n1) = xc.reshaped(n2, n1).template cast<Complex>(); + } else { + X = xc.reshaped(n2, n1).template cast<Complex>(); + } + ldexpInPlace(X, -ex); + fft2(X); + X.array() *= sUse.array(); + ifft2(X); + if (padded) Xn = X.topLeftCorner(n2, n1); + ComplexArray& out = padded ? Xn : X; + ldexpInPlace(out, ex + es); + if (accumulate) + dst.col(k) += alpha * internal::structured_scalar_part_impl<ProductScalar>::run(out.reshaped()); + else + dst.col(k) = alpha * internal::structured_scalar_part_impl<ProductScalar>::run(out.reshaped()); + } + } + + /** \internal Multiplies every entry of \a X by 2^e, exactly. The per-entry + * ldexp saturates to zero / infinity component-wise without ever forming the + * (possibly unrepresentable) scale factor 2^e itself. */ + static void ldexpInPlace(ComplexArray& X, int e) { + if (e == 0) return; + X.realView() = X.realView().array().ldexp(e).matrix(); + } + + /** \internal In-place forward or inverse 2-D FFT by the row-column algorithm: + * transform every column (length n2), then every row (length n1). Length-1 + * transforms are the identity (and unsupported by the kissfft backend), hence + * the guards. */ + template <bool Inverse> + void transform2(ComplexArray& X) const { + auto&& fft = internal::structured_fft_engine<RealScalar>(); + const Index n2 = X.rows(), n1 = X.cols(); + // Both passes share one output buffer (resized by the engine) and one packed + // input buffer, allocated at most twice for the whole transform. A column of + // the column-major grid is already packed and goes to the engine as is; a row + // is strided, and the engine would pack it into a fresh temporary per call. + ComplexVector tmp, rowv; + if (n2 > 1) { + for (Index k1 = 0; k1 < n1; ++k1) { + EIGEN_IF_CONSTEXPR (Inverse) { + fft.inv(tmp, X.col(k1), n2); + } else { + fft.fwd(tmp, X.col(k1), n2); + } + X.col(k1) = tmp; + } + } + if (n1 > 1) { + for (Index k2 = 0; k2 < n2; ++k2) { + rowv = X.row(k2).transpose(); + EIGEN_IF_CONSTEXPR (Inverse) { + fft.inv(tmp, rowv, n1); + } else { + fft.fwd(tmp, rowv, n1); + } + X.row(k2) = tmp.transpose(); + } + } + } + + /** \internal In-place 2-D forward FFT. */ + void fft2(ComplexArray& X) const { transform2<false>(X); } + + /** \internal In-place 2-D inverse FFT. */ + void ifft2(ComplexArray& X) const { transform2<true>(X); } + + /** \internal Reconstructs the generating array from an exact-size symbol. */ + GeneratorType generatorFromSymbol(ComplexArray symbol) const { + ifft2(symbol); + GeneratorType generator = internal::structured_scalar_part_impl<Scalar>::run(symbol); + return generator; + } + + /** \internal \returns the 2-D DFT of the generating array. */ + ComplexArray computeSymbol() const { + ComplexArray s = m_g.template cast<Complex>(); + fft2(s); + return s; + } + + /** \internal \returns \a M index-reversed in both dimensions: + * result(k2, k1) = M((-k2) mod n2, (-k1) mod n1). Row 0 and column 0 stay in + * place; the rest is a two-dimensional reversal. Empty input stays empty. */ + template <typename MatType> + static MatType reverse2(const MatType& M) { + const Index r = M.rows(), c = M.cols(); + MatType R(r, c); + if (M.size() == 0) return R; + R(0, 0) = M(0, 0); + if (c > 1) R.row(0).tail(c - 1) = M.row(0).tail(c - 1).reverse(); + if (r > 1) R.col(0).tail(r - 1) = M.col(0).tail(r - 1).reverse(); + if (r > 1 && c > 1) R.bottomRightCorner(r - 1, c - 1) = M.bottomRightCorner(r - 1, c - 1).reverse(); + return R; + } + + /** \internal Computes the moduli of the symbol entries and the matching + * rank/pseudo-inversion threshold, both evaluated in an exactly rescaled frame + * [5]: the entries are pre-scaled by a power of two chosen so no modulus can + * overflow. A finite complex entry near the overflow threshold has a + * non-representable modulus, which would otherwise turn the threshold into + * infinity and misclassify every other entry (rank under-reported, solve() + * zeroing valid Fourier modes). The rescaling is exact, so comparing scaled + * moduli against the scaled threshold is equivalent to the unscaled + * comparison. The threshold keeps the N * epsilon * max|s| convention of [3], + * chapter 5.4, and the smallest-normal clamp -- carried into the scaled frame, + * where its underflowing to zero for a huge frame is correct: no entry of such + * a symbol can sit below the smallest normal number. Entries at or above the + * threshold, in particular a smallest-normal entry of a moderate symbol, are + * inverted (their reciprocals are finite). */ + static void scaledModuli(const ComplexArray& s, RealArray& mods, RealScalar& tol) { + const int e = numext::maxi(internal::structured_exponent_bound(s), 0); + // Two exact factors, as in Circulant::scaledModuli(): a single 2^-e is itself + // subnormal once the frame exceeds the exponent range, and reads as zero under + // flush-to-zero -- collapsing every scaled modulus, and the threshold with + // them. Neither factor exceeds one, so no intermediate underflows on its own. + const RealScalar down1 = numext::ldexp(RealScalar(1), -(e / 2)), down2 = numext::ldexp(RealScalar(1), -(e - e / 2)); + mods = ((s * down1) * down2).cwiseAbs(); + tol = numext::maxi(RealScalar(s.size()) * NumTraits<RealScalar>::epsilon() * mods.maxCoeff(), + ((std::numeric_limits<RealScalar>::min)() * down1) * down2); + } + + /** \internal Writes the unit-norm 2-D Fourier eigenvector of frequencies + * \c (f1, f2) into column \a dstCol of \a F: entry \c b1*n2 + i2 is + * exp(2 pi i (b1 f1 / n1 + i2 f2 / n2)) / sqrt(N). All frequency products are + * accumulated incrementally modulo their length, so the angles stay O(2 pi) at + * full accuracy and no Index overflow can occur. */ + void fourierColumn(ComplexMatrix& F, Index f1, Index f2, Index dstCol) const { + const Index n2 = blockSize(), n1 = numBlocks(); + const RealScalar scale = RealScalar(1) / numext::sqrt(RealScalar(rows())); + ComplexVector w2(n2); + Index jf = 0; // i2 * f2 mod n2 + for (Index i2 = 0; i2 < n2; ++i2) { + w2[i2] = std::polar(scale, RealScalar(2 * EIGEN_PI) * RealScalar(jf) / RealScalar(n2)); + jf += f2; + if (jf >= n2) jf -= n2; + } + Index bf = 0; // b1 * f1 mod n1 + for (Index b1 = 0; b1 < n1; ++b1) { + const Complex w1 = std::polar(RealScalar(1), RealScalar(2 * EIGEN_PI) * RealScalar(bf) / RealScalar(n1)); + F.col(dstCol).segment(b1 * n2, n2) = w1 * w2; + bf += f1; + if (bf >= n1) bf -= n1; + } + } + + /** \internal Projects the complex determinant onto \c Scalar. */ + static Scalar toScalar(const Complex& z, std::true_type /*scalar_is_real*/) { return numext::real(z); } + static Scalar toScalar(const Complex& z, std::false_type /*scalar_is_real*/) { return z; } + + GeneratorType m_g; + ComplexArray m_symbol; + // The padded embedding symbol products use when a block dimension is not + // 5-smooth; empty otherwise. See computeProdSymbol(). + ComplexArray m_prodSymbol; + bool m_fftUsable; +}; + +/** \ingroup StructuredMatrices_Module + * \returns a \ref Bccb operator with generating array \a generator. The + * compile-time dimensions of the operator are deduced from the array. */ +template <typename Derived> +Bccb<typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime> makeBccb( + const MatrixBase<Derived>& generator) { + return Bccb<typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime>(generator); +} + +namespace internal { + +template <typename Scalar_, int BlockSize_, int NumBlocks_, typename Rhs, int ProductTag> +struct generic_product_impl<Bccb<Scalar_, BlockSize_, NumBlocks_>, Rhs, StructuredShape, DenseShape, ProductTag> + : structured_product_impl<Bccb<Scalar_, BlockSize_, NumBlocks_>, Rhs> {}; + +} // namespace internal + +} // namespace Eigen + +#endif // EIGEN_STRUCTURED_BCCB_H
diff --git a/contrib/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h b/contrib/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h index 510ba6a..a2d368c 100644 --- a/contrib/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h +++ b/contrib/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h
@@ -11,7 +11,7 @@ // [1] N. J. Higham, "Accuracy and Stability of Numerical Algorithms", 2nd ed., // SIAM, 2002, chapter 27. Avoiding spurious overflow by rescaling with // powers of two, the technique behind structured_exponent_bound() and the -// column scaling in structured_fft_apply(). +// column scaling in structured_fft_apply() and Bccb::applySymbol(). // [2] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974. // Scaling by a power of two is exact, so the scaled transforms introduce // no roundoff beyond the transforms themselves.
diff --git a/contrib/benchmarks/StructuredMatrices/CMakeLists.txt b/contrib/benchmarks/StructuredMatrices/CMakeLists.txt index 5eb6524..57d8ce0 100644 --- a/contrib/benchmarks/StructuredMatrices/CMakeLists.txt +++ b/contrib/benchmarks/StructuredMatrices/CMakeLists.txt
@@ -11,3 +11,4 @@ eigen_add_benchmark(bench_structured_vandermonde bench_structured_vandermonde.cpp) eigen_add_benchmark(bench_structured_cauchy bench_structured_cauchy.cpp) eigen_add_benchmark(bench_structured_dpr1 bench_structured_dpr1.cpp) +eigen_add_benchmark(bench_structured_bccb bench_structured_bccb.cpp)
diff --git a/contrib/benchmarks/StructuredMatrices/bench_structured_bccb.cpp b/contrib/benchmarks/StructuredMatrices/bench_structured_bccb.cpp new file mode 100644 index 0000000..30d71f8 --- /dev/null +++ b/contrib/benchmarks/StructuredMatrices/bench_structured_bccb.cpp
@@ -0,0 +1,96 @@ +// Benchmarks for the Bccb operator: the O(N log N) 2-D-FFT-based product and +// direct solve against their dense counterparts (GEMV, and a solve with a +// precomputed PartialPivLU factorization -- the analogue of Bccb's symbol, +// which is precomputed at construction). The generating array is n x n, so the +// operator is N x N with N = n^2; the dense variants stop at n = 64 (N = 4096) +// to keep the dense matrix and its factorization affordable. +// 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; + +// Diagonally dominant generating array: keeps the symbol away from zero, so the +// solves are well conditioned. +static Mat generatingArray(Index n) { + Mat G = Mat::Random(n, n); + G(0, 0) += double(2 * n * n); + return G; +} + +// The dense representation, built entry-wise: entry (i,j) with i = b1*n2 + i2, +// j = c1*n2 + j2 is G((i2-j2) mod n2, (b1-c1) mod n1). +static Mat denseBccb(const Mat& G) { + const Index n2 = G.rows(), n1 = G.cols(), N = n1 * n2; + Mat dense(N, N); + for (Index j = 0; j < N; ++j) + for (Index i = 0; i < N; ++i) { + Index k2 = i % n2 - j % n2; + if (k2 < 0) k2 += n2; + Index k1 = i / n2 - j / n2; + if (k1 < 0) k1 += n1; + dense(i, j) = G(k2, k1); + } + return dense; +} + +// --- product: y = C * x --- +static void BM_BccbProduct(benchmark::State& state) { + const Index n = state.range(0), N = n * n; + Mat G = generatingArray(n); + Bccb<double> C(G); + Vec x = Vec::Random(N), y(N); + for (auto _ : state) { + y.noalias() = C * x; + benchmark::DoNotOptimize(y.data()); + } +} +// 97 is prime in both block dimensions: the product pins the padded embedding +// (the exact-size transform would run kissfft's quadratic generic butterfly). +BENCHMARK(BM_BccbProduct)->Arg(8)->Arg(16)->Arg(32)->Arg(64)->Arg(97)->Arg(128); + +static void BM_DenseProduct(benchmark::State& state) { + const Index n = state.range(0), N = n * n; + Mat dense = denseBccb(generatingArray(n)); + Vec x = Vec::Random(N), y(N); + for (auto _ : state) { + y.noalias() = dense * x; + benchmark::DoNotOptimize(y.data()); + } +} +BENCHMARK(BM_DenseProduct)->Arg(8)->Arg(16)->Arg(32)->Arg(64); + +// --- solve: x = C^{-1} * b, both sides reusing their precomputed factorization +// (the 2-D DFT symbol for Bccb, the LU factors for the dense matrix) --- +static void BM_BccbSolve(benchmark::State& state) { + const Index n = state.range(0), N = n * n; + Mat G = generatingArray(n); + Bccb<double> C(G); + Vec b = Vec::Random(N), x(N); + for (auto _ : state) { + x = C.solve(b); + benchmark::DoNotOptimize(x.data()); + } +} +// All 5-smooth: solve transforms at the exact size, so a prime dimension would +// only time kissfft's generic butterfly instead of the operator. +BENCHMARK(BM_BccbSolve)->Arg(8)->Arg(16)->Arg(32)->Arg(64)->Arg(96)->Arg(128); + +static void BM_DenseSolve(benchmark::State& state) { + const Index n = state.range(0), N = n * n; + Mat dense = denseBccb(generatingArray(n)); + PartialPivLU<Mat> lu(dense); + Vec b = Vec::Random(N), x(N); + for (auto _ : state) { + x = lu.solve(b); + benchmark::DoNotOptimize(x.data()); + } +} +BENCHMARK(BM_DenseSolve)->Arg(8)->Arg(16)->Arg(32)->Arg(64);
diff --git a/contrib/test/CMakeLists.txt b/contrib/test/CMakeLists.txt index 0f9879c..d18e077 100644 --- a/contrib/test/CMakeLists.txt +++ b/contrib/test/CMakeLists.txt
@@ -51,6 +51,7 @@ ei_add_test(alignedvector3) ei_add_test(FFT) +ei_add_test(structured_bccb) ei_add_test(EulerAngles)
diff --git a/contrib/test/structured_bccb.cpp b/contrib/test/structured_bccb.cpp new file mode 100644 index 0000000..d5aa38e --- /dev/null +++ b/contrib/test/structured_bccb.cpp
@@ -0,0 +1,938 @@ +// 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 "fp_control.h" + +#include <contrib/Eigen/StructuredMatrices> + +using namespace Eigen; + +// Dense inverse DFT matrix, used to synthesize generators independently of the +// FFT implementation under test. +template <typename RealScalar> +Matrix<std::complex<RealScalar>, Dynamic, Dynamic> inverse_dft_matrix(Index n) { + typedef std::complex<RealScalar> Complex; + Matrix<Complex, Dynamic, Dynamic> F(n, n); + for (Index a = 0; a < n; ++a) + for (Index c = 0; c < n; ++c) + F(a, c) = + std::polar(RealScalar(1) / RealScalar(n), RealScalar(2 * EIGEN_PI) * RealScalar((a * c) % n) / RealScalar(n)); + return F; +} + +// Reference dense BCCB built entry-wise from the generating array, independently +// of the operator under test: entry (i,j) with i = b1*n2+i2, j = c1*n2+j2 is +// G((i2-j2) mod n2, (b1-c1) mod n1). +template <typename Scalar> +Matrix<Scalar, Dynamic, Dynamic> reference_bccb(const Matrix<Scalar, Dynamic, Dynamic>& G) { + const Index n2 = G.rows(), n1 = G.cols(), N = n1 * n2; + Matrix<Scalar, Dynamic, Dynamic> dense(N, N); + for (Index j = 0; j < N; ++j) + for (Index i = 0; i < N; ++i) { + Index k2 = i % n2 - j % n2; + if (k2 < 0) k2 += n2; + Index k1 = i / n2 - j / n2; + if (k1 < 0) k1 += n1; + dense(i, j) = G(k2, k1); + } + return dense; +} + +template <typename Scalar> +void test_bccb_product(Index n2, Index n1) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + VERIFY_IS_EQUAL(C.rows(), N); + + Mat Cd = C; + VERIFY_IS_APPROX(Cd, dense); + Mat accumulated = Mat::Random(N, N); + const Mat initial = accumulated; + accumulated += C; + VERIFY_IS_APPROX(accumulated, (initial + dense).eval()); + accumulated = initial; + accumulated -= C; + VERIFY_IS_APPROX(accumulated, (initial - dense).eval()); + for (Index t = 0; t < 5; ++t) { + Index i = internal::random<Index>(0, N - 1), j = internal::random<Index>(0, N - 1); + VERIFY_IS_APPROX(C.coeff(i, j), dense(i, j)); + } + + Vec x = Vec::Random(N); + VERIFY_IS_APPROX((C * x).eval(), (dense * x).eval()); + + Mat X = Mat::Random(N, 3); + VERIFY_IS_APPROX((C * X).eval(), (dense * X).eval()); + + // Accumulation form exercised by the iterative solvers. + Vec y = Vec::Random(N); + Vec y0 = y; + y.noalias() += C * x; + VERIFY_IS_APPROX(y, (y0 + dense * x).eval()); +} + +template <typename Scalar> +void test_bccb_transpose(Index n2, Index n1) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + + 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 x = Vec::Random(N); + VERIFY_IS_APPROX((C.transpose() * x).eval(), (dense.transpose() * x).eval()); + VERIFY_IS_APPROX((C.adjoint() * x).eval(), (dense.adjoint() * x).eval()); + + // Exact round trips: generators and symbols are pure permutations/conjugations. + Bccb<Scalar> Ctt = C.transpose().transpose(); + VERIFY_IS_EQUAL(Ctt.generator(), G); + VERIFY_IS_EQUAL(Ctt.symbol(), C.symbol()); + Bccb<Scalar> Caa = C.adjoint().adjoint(); + VERIFY_IS_EQUAL(Caa.generator(), G); + VERIFY_IS_EQUAL(Caa.symbol(), C.symbol()); +} + +// 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_bccb_delayed_product(Index n2, Index n1) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + STATIC_CHECK(!std::is_reference<typename internal::ref_selector<Bccb<Scalar>>::type>::value); + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + Vec x = Vec::Random(N); + + auto expr = C.adjoint() * x; // the adjoint temporary dies with the full expression + Vec scribble = Vec::Random(2 * N); // reuses the temporary's freed heap storage + Vec y = expr; + VERIFY_IS_APPROX(y, (dense.adjoint() * x).eval()); + VERIFY_IS_EQUAL(scribble.size(), 2 * N); // keep the scribble alive across the evaluation +} + +// The products carry the default product tag, so plain assignment materializes +// a temporary exactly like a dense product: x = C * x and x += C * x get the +// ordinary dense-product aliasing semantics. +template <typename Scalar> +void test_bccb_aliased_product(Index n2, Index n1) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + + Vec x = Vec::Random(N); + Vec y = x; + y = C * y; + VERIFY_IS_APPROX(y, (dense * x).eval()); + + y = x; + y += C * y; + VERIFY_IS_APPROX(y, (x + dense * x).eval()); + + y = x; + y -= C * y; + VERIFY_IS_APPROX(y, (x - dense * x).eval()); + + Mat X = Mat::Random(N, 3); + Mat Y = X; + Y = C * Y; + VERIFY_IS_APPROX(Y, (dense * X).eval()); +} + +// Aliasing beyond the same-object case: the default-product temporary must also +// resolve right-hand-side expressions that reference the destination and +// overlapping views of one buffer, neither of which is_same_dense can see. +template <typename Scalar> +void test_bccb_aliased_expression(Index n2, Index n1) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + + // Right-hand-side expression referencing the destination. + Vec x = Vec::Random(N), x0 = x; + x = C * (x + Vec::Ones(N)); + VERIFY_IS_APPROX(x, (dense * (x0 + 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); +} + +// 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_bccb_mixed_scalar(Index n2, Index n1) { + typedef std::complex<RealScalar> Complex; + typedef Matrix<RealScalar, Dynamic, 1> RVec; + typedef Matrix<RealScalar, Dynamic, Dynamic> RMat; + typedef Matrix<Complex, Dynamic, 1> CVec; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + RMat G = RMat::Random(n2, n1); + Bccb<RealScalar> C(G); + CMat dense = reference_bccb<RealScalar>(G).template cast<Complex>(); + const Index N = n1 * n2; + + CVec x = CVec::Random(N); + CVec y = C * x; + VERIFY_IS_APPROX(y, (dense * x).eval()); + + CVec y0 = CVec::Random(N); + y = y0; + y.noalias() += C * x; + VERIFY_IS_APPROX(y, (y0 + dense * x).eval()); + + CMat Gc = CMat::Random(n2, n1); + Bccb<Complex> Cc(Gc); + CMat denseC = reference_bccb<Complex>(Gc); + RVec xr = RVec::Random(N); + CVec z = Cc * xr; + VERIFY_IS_APPROX(z, (denseC * xr).eval()); +} + +// Non-5-smooth block dimensions: products transform on the padded per-axis +// circulant-embedding grid (kissfft's generic butterfly is quadratic in prime +// factors), while the spectral operations keep the exact-size symbol. Cover +// each padded combination -- n2 awkward, n1 awkward, both -- against the dense +// reference, including the transposition family (which reuses the padded +// symbol through the reversal/conjugation rules), solve, and inverse (which +// rebuilds the padded symbol for the inverse operator). +template <typename Scalar> +void test_bccb_prime_dimensions(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + const Index N = n1 * n2; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + + Vec x = Vec::Random(N); + VERIFY_IS_APPROX((C * x).eval(), (dense * x).eval()); + Mat X = Mat::Random(N, 3); + VERIFY_IS_APPROX((C * X).eval(), (dense * X).eval()); + + VERIFY_IS_APPROX((C.transpose() * x).eval(), (dense.transpose() * x).eval()); + VERIFY_IS_APPROX((C.adjoint() * x).eval(), (dense.adjoint() * x).eval()); + VERIFY_IS_APPROX((C.conjugate() * x).eval(), (dense.conjugate() * x).eval()); + + // Exact symbol round trip through the transposition family: both cached + // symbols are pure permutations/conjugations of the originals. + Bccb<Scalar> Ctt = C.transpose().transpose(); + VERIFY_IS_EQUAL(Ctt.generator(), G); + VERIFY_IS_EQUAL(Ctt.symbol(), C.symbol()); + + // Solve (exact-size spectral path) and inverse (whose product rebuilds the + // padded symbol): diagonal dominance keeps the symbol away from zero. + Mat Gd = Mat::Random(n2, n1); + Gd(0, 0) += Scalar(RealScalar(2 * N)); + Bccb<Scalar> Cd(Gd); + Mat densed = reference_bccb<Scalar>(Gd); + Vec b = Vec::Random(N); + Vec xs = Cd.solve(b); + VERIFY_IS_APPROX((densed * xs).eval(), b); + VERIFY_IS_APPROX((Cd.inverse() * b).eval(), xs); +} + +// Dynamic dimension mismatches must trip the runtime assertion when the product +// or solve expression is built. (Incompatible *fixed* sizes are rejected at +// compile time by a static assertion in operator* / solve, which a runtime test +// cannot exercise.) +void test_bccb_dimension_asserts() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(3, 4); + Bccb<double> C(G); + Vec bad = Vec::Random(11); + VERIFY_RAISES_ASSERT(C * bad); + VERIFY_RAISES_ASSERT(C.solve(bad)); +} + +template <typename Scalar> +void test_bccb_solve(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + // Diagonal dominance (through the (0,0) generator entry) keeps the symbol away + // from zero, so the direct 2-D FFT solve is exact up to roundoff. + Mat G = Mat::Random(n2, n1); + G(0, 0) += Scalar(RealScalar(2 * n1 * n2)); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + + Vec b = Vec::Random(N); + Vec x = C.solve(b); + VERIFY_IS_APPROX((dense * x).eval(), b); + + Mat B = Mat::Random(N, 3); + Mat Xs = C.solve(B); + VERIFY_IS_APPROX((dense * Xs).eval(), B); +} + +// Rank-deficient BCCB synthesized by zeroing 2-D symbol entries: the rank counts +// the surviving entries and solve() matches the SVD pseudo-inverse. +template <typename Scalar> +void test_bccb_rank_deficient(Index n2, Index n1, Index defect) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> Complex; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + const Index N = n1 * n2; + CMat S = CMat::Random(n2, n1); + S.array() += Complex(2); // keep the surviving moduli away from the threshold + for (Index k = 0; k < defect; ++k) S((2 * k + 1) % n2, (3 * k) % n1) = Complex(0); + + const CMat F2 = inverse_dft_matrix<RealScalar>(n2), F1 = inverse_dft_matrix<RealScalar>(n1); + Mat G = (F2 * S * F1.transpose()).eval(); // Scalar is complex here + Bccb<Scalar> C(G); + VERIFY_IS_EQUAL(C.rank(), N - defect); + + Mat dense = reference_bccb<Scalar>(G); + JacobiSVD<Mat> svd(dense, ComputeThinU | ComputeThinV); + Vec b = Vec::Random(N); + VERIFY_IS_APPROX(C.solve(b), svd.solve(b).eval()); +} + +template <typename Scalar> +void test_bccb_zero(Index n2, Index n1) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Bccb<Scalar> C(Mat(Mat::Zero(n2, n1))); + VERIFY_IS_EQUAL(C.rank(), 0); + VERIFY(C.singularValues().isZero()); + Vec b = Vec::Random(n1 * n2); + VERIFY(C.solve(b).isZero()); +} + +void test_bccb_nan_propagation(Index n2, Index n1) { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(n2, n1); + G(n2 / 2, n1 / 2) = std::numeric_limits<double>::quiet_NaN(); + Bccb<double> C(G); + VERIFY_IS_EQUAL(C.rank(), n1 * n2); + Vec b = Vec::Random(n1 * n2); + Vec x = C.solve(b); + VERIFY(!(x.array() == x.array()).all()); +} + +// Inputs at the very top (and bottom) of the exponent range: the transforms sum +// up to N terms, so without scaling the FFT intermediates overflow even though +// every entry of the true result is representable. The apply path rescales the +// symbol and each right-hand-side column by exact powers of two and folds the +// exponent back at the end. +void test_bccb_finite_overflow() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + + const double huge = (std::numeric_limits<double>::max)() / 16; + // The 2-D FFT round trip costs a few ulps per entry; a norm-based comparison + // would itself overflow at these magnitudes, hence the entrywise bound. + const double kFftRoundTripTol = 100 * NumTraits<double>::epsilon(); + + // Reviewer repro: a 36x36 identity BCCB applied to a vector of DBL_MAX/16 + // must return the vector, not 36 NaNs. + Mat G = Mat::Zero(6, 6); + G(0, 0) = 1.0; + Bccb<double> C(G); + Vec x = Vec::Constant(36, huge); + Vec y = C * x; + VERIFY(y.allFinite()); + VERIFY(((y - x).array().abs() <= kFftRoundTripTol * x.array().abs()).all()); + + // A huge generator makes a huge symbol: the scaling must come from the symbol + // side as well. + Mat Gh = Mat::Zero(6, 6); + Gh(0, 0) = huge; // C == huge * Identity + Bccb<double> Ch(Gh); + Vec ones = Vec::Ones(36); + Vec z = Ch * ones; + VERIFY(z.allFinite()); + VERIFY(((z.array() - huge).abs() <= kFftRoundTripTol * huge).all()); + + // The inverse symbol of a tiny operator is huge; solve() shares the scaled + // apply path. + Mat Gt = Mat::Zero(6, 6); + Gt(0, 0) = 1.0 / huge; + Bccb<double> Ct(Gt); + Vec w = Ct.solve(ones); + VERIFY(w.allFinite()); + VERIFY(((w.array() - huge).abs() <= kFftRoundTripTol * huge).all()); + + // An exactly zero column maps to an exactly zero column (short-circuited + // before any scaling), while nonzero columns are still transformed. + Mat B(36, 2); + B.col(0).setConstant(huge); + B.col(1).setZero(); + Mat Y = C * B; + VERIFY(Y.col(1).isZero()); + VERIFY(((Y.col(0) - x).array().abs() <= kFftRoundTripTol * x.array().abs()).all()); + + // Genuine NaN and infinity inputs keep propagating; such right-hand sides + // take the direct kernel (see test_bccb_nonfinite_product for the entrywise + // semantics). + Vec xn = x; + xn[7] = std::numeric_limits<double>::quiet_NaN(); + Vec yn = C * xn; + VERIFY(!(yn.array() == yn.array()).all()); + Vec xi = Vec::Ones(36); + xi[3] = std::numeric_limits<double>::infinity(); + Vec yi = C * xi; + VERIFY(!yi.allFinite()); +} + +// The scaling exponents are derived from component-wise magnitudes: a finite +// complex value near the overflow threshold has a non-representable modulus, +// which would otherwise disable the scaling and turn an exactly representable +// product into NaN. +template <typename RealScalar> +void test_bccb_fft_complex_boundary(Index n2, Index n1) { + typedef std::complex<RealScalar> Complex; + typedef Matrix<RealScalar, Dynamic, Dynamic> RMat; + typedef Matrix<Complex, Dynamic, 1> CVec; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + const RealScalar kFftRoundTripTol = RealScalar(100) * NumTraits<RealScalar>::epsilon(); + const RealScalar big = RealScalar(0.75) * (std::numeric_limits<RealScalar>::max)(); + const Index N = n1 * n2; + + // Identity BCCB with a complex generator: the product returns the right-hand + // side unchanged even though |x_k| overflows. + CMat G = CMat::Zero(n2, n1); + G(0, 0) = Complex(1); + Bccb<Complex> C(G); + CVec x = CVec::Constant(N, Complex(big, big)); + CVec y = C * x; + VERIFY(y.allFinite()); + VERIFY(((y - x).cwiseAbs() / big).maxCoeff() <= kFftRoundTripTol); + + // A real identity operator applied to the same complex right-hand side takes + // the mixed-scalar product path. + RMat Gr = RMat::Zero(n2, n1); + Gr(0, 0) = RealScalar(1); + Bccb<RealScalar> Cr(Gr); + y = Cr * x; + VERIFY(y.allFinite()); + VERIFY(((y - x).cwiseAbs() / big).maxCoeff() <= kFftRoundTripTol); +} + +// Entrywise IEEE comparison for the non-finite tests: NaNs match NaNs, +// infinities match by value (sign included), finite entries match to roundoff. +// VERIFY_IS_APPROX would reject any output containing NaN. +template <typename D1, typename D2> +bool ieee_entrywise_match(const D1& a, const D2& b) { + if (a.rows() != b.rows() || a.cols() != b.cols()) return false; + for (Index j = 0; j < a.cols(); ++j) + for (Index i = 0; i < a.rows(); ++i) { + const typename D1::Scalar x = a(i, j), y = b(i, j); + if (x == y) continue; // finite match or same-signed infinities + if ((numext::isnan)(x) && (numext::isnan)(y)) continue; // both NaN + if (!test_isApprox(x, y)) return false; // finite roundoff + } + return true; +} + +// Scalar-loop product: the mathematically transparent IEEE reference for the +// non-finite tests. Eigen's own vectorized complex kernels can smear a single +// infinity into NaN (Inf - Inf across the split real/imaginary accumulators), so +// the dense product is not a faithful entrywise reference for non-finite data. +template <typename Scalar> +Matrix<Scalar, Dynamic, 1> reference_product_ieee(const Matrix<Scalar, Dynamic, Dynamic>& A, + const Matrix<Scalar, Dynamic, 1>& x) { + Matrix<Scalar, Dynamic, 1> y(A.rows()); + for (Index i = 0; i < A.rows(); ++i) { + Scalar acc(0); + for (Index j = 0; j < A.cols(); ++j) acc += A(i, j) * x[j]; + y[i] = acc; + } + return y; +} + +// A single Inf or NaN in the data must propagate like the reference product -- +// through the dot products that touch it -- instead of being smeared into NaNs +// across the whole output by the transforms. +template <typename Scalar> +void test_bccb_nonfinite_product(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + const RealScalar inf = std::numeric_limits<RealScalar>::infinity(); + const RealScalar nan = std::numeric_limits<RealScalar>::quiet_NaN(); + const Index N = n1 * n2; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + + // Inf in the right-hand side. + Vec x = Vec::Random(N); + x[N / 2] = Scalar(inf); + VERIFY(ieee_entrywise_match((C * x).eval(), reference_product_ieee(dense, x))); + + // NaN in the right-hand side. + Vec xn = Vec::Random(N); + xn[N - 1] = Scalar(nan); + VERIFY(ieee_entrywise_match((C * xn).eval(), reference_product_ieee(dense, xn))); + + // Mixed multi-column right-hand side: the non-finite column falls back to the + // direct kernel individually while the finite column keeps the FFT path. + Mat Xm(N, 2); + Xm.col(0) = Vec::Random(N); + Xm.col(1) = x; + Mat Ym = C * Xm; + VERIFY_IS_APPROX(Ym.col(0).eval(), (dense * Xm.col(0)).eval()); + VERIFY(ieee_entrywise_match(Ym.col(1).eval(), reference_product_ieee(dense, Vec(Xm.col(1))))); + + // A zero column carrying a single NaN must not take the zero-column shortcut: + // the fast-max routing scan can miss a NaN among zeros (an Inf always + // surfaces), so the shortcut rechecks exactly and such a column falls back. + Vec xz = Vec::Zero(N); + xz[0] = Scalar(nan); + VERIFY(ieee_entrywise_match((C * xz).eval(), reference_product_ieee(dense, xz))); + + // Inf in the generating array: the operator itself is non-finite, whatever the + // right-hand side. + Mat G2 = Mat::Random(n2, n1); + G2(n2 / 2, n1 / 2) = Scalar(-inf); + Bccb<Scalar> C2(G2); + Mat dense2 = reference_bccb<Scalar>(G2); + Vec x2 = Vec::Random(N); + VERIFY(ieee_entrywise_match((C2 * x2).eval(), reference_product_ieee(dense2, x2))); + + // Non-finite right-hand sides of solve() apply the pseudo-inverse -- itself a + // BCCB operator -- through the direct kernel, so the Inf propagates entrywise + // instead of NaN-ing the whole output through the transforms. The operator is + // diagonally dominant, so no symbol entry is thresholded and the pseudo-inverse + // coincides with inverse(), whose dense form serves as the reference matrix. + Mat Gd = Mat::Random(n2, n1); + Gd(0, 0) += Scalar(RealScalar(2 * N)); + Bccb<Scalar> Cd(Gd); + Mat pinv = Mat(Cd.inverse()); + Vec binf = Vec::Random(N); + binf[N / 3] = Scalar(inf); + VERIFY(ieee_entrywise_match(Cd.solve(binf), reference_product_ieee(pinv, binf))); +} + +// An all-zero right-hand side must not short-circuit to an exact zero when the +// operator holds non-finite data: every row of a BCCB touches every generator +// entry, so each result entry is a dot product with an Inf coefficient times +// zero -- NaN under IEEE, not zero. +void test_bccb_nonfinite_zero_rhs(Index n2, Index n1) { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + const Index N = n1 * n2; + + Mat G = Mat::Random(n2, n1); + G(n2 / 2, n1 / 2) = std::numeric_limits<double>::infinity(); + Bccb<double> C(G); + Mat dense = reference_bccb<double>(G); + + const Vec z = Vec::Zero(N); + Vec y = C * z; + VERIFY(y.hasNaN()); + VERIFY(ieee_entrywise_match(y, reference_product_ieee(dense, z))); + + // The same holds for solve(): the symbol of a non-finite operator holds NaN, + // so the pseudo-inverse applied to a zero right-hand side is NaN, not zero. + Vec xs = C.solve(z); + VERIFY(xs.hasNaN()); +} + +// Spectra with a wide dynamic range: the determinant is representable, but a +// plain product of the eigenvalues in FFT order overflows to infinity (or +// underflows to an exact zero) partway through. Pins the balanced accumulation +// in determinant(). +void test_bccb_determinant_scaled() { + typedef std::complex<double> Complex; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + const Index n2 = 25, n1 = 40; // N = 1000 + const Index nLead = 653; + + // Symbol magnitudes `lead` on the first nLead column-major-flattened indices + // (the accumulation order of determinant()) and `rest` elsewhere, so the + // naive running product leaves the representable range partway through while + // the true determinant lead^nLead * rest^(N - nLead) is representable. The + // generator is recovered through the dense inverse 2-D DFT matrices, + // independently of the implementation under test. + auto makeOperator = [n2, n1](double lead, double rest) { + CMat S(n2, n1); + for (Index k1 = 0; k1 < n1; ++k1) + for (Index k2 = 0; k2 < n2; ++k2) S(k2, k1) = Complex(k1 * n2 + k2 < nLead ? lead : rest); + const CMat F2 = inverse_dft_matrix<double>(n2), F1 = inverse_dft_matrix<double>(n1); + return Bccb<Complex>((F2 * S * F1.transpose()).eval()); + }; + + // The spectrum is only reproduced up to the FFT round trip's forward error, + // and the determinant multiplies ~1e3 such factors. + const double kSpectrumRoundTripTol = 1e8 * NumTraits<double>::epsilon(); + { + // det = 10^653 * 10^-347 = 1e306; the naive partial product reaches 1e327. + Bccb<Complex> C = makeOperator(10.0, 0.1); + const Complex det = C.determinant(); + VERIFY((numext::isfinite)(numext::abs(det))); + VERIFY(numext::abs(det / 1e306 - Complex(1)) <= kSpectrumRoundTripTol); + } + { + // det = 10^-653 * 10^347 = 1e-306; the naive partial product reaches + // 1e-327, well below the smallest subnormal, and flushes to an exact zero. + Bccb<Complex> C = makeOperator(0.1, 10.0); + const Complex det = C.determinant(); + VERIFY(det != Complex(0)); + VERIFY(numext::abs(det / 1e-306 - Complex(1)) <= kSpectrumRoundTripTol); + } + { + // A genuinely overflowing determinant must still saturate to infinity. + typedef Matrix<double, Dynamic, Dynamic> Mat; + Mat G = Mat::Zero(6, 6); + G(0, 0) = 1e308; // the symbol is 1e308 everywhere: det = 10^11088 + Bccb<double> C(G); + VERIFY((numext::isinf)(C.determinant())); + } +} + +// The rank decision at the clamped threshold (the smallest normal number): a +// smallest-normal symbol entry is still inverted -- the comparison is strict, +// matching SVDBase::rank(), which likewise reports rank one -- a subnormal entry +// is treated as an exact zero, and non-finite entries count as non-zero. +void test_bccb_rank_boundaries() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + + const double mn = (std::numeric_limits<double>::min)(); + const Vec b = Vec::Ones(1); + { + Bccb<double> C(Mat(Mat::Constant(1, 1, mn))); + VERIFY_IS_EQUAL(C.rank(), 1); + Vec x = C.solve(b); + VERIFY((numext::isfinite)(x[0])); + VERIFY_IS_APPROX(x[0], 1.0 / mn); + } + { + Bccb<double> C(Mat(Mat::Constant(1, 1, mn / 2))); // subnormal + VERIFY_IS_EQUAL(C.rank(), 0); + VERIFY(C.solve(b).isZero()); + } + { + Bccb<double> C(Mat(Mat::Constant(1, 1, std::numeric_limits<double>::infinity()))); + VERIFY_IS_EQUAL(C.rank(), 1); + VERIFY((numext::isinf)(C.determinant())); + } +} + +// A finite complex symbol entry near the overflow threshold has a +// non-representable modulus. The rank threshold used to be computed from the raw +// moduli, turning it into infinity: the rank was under-reported and solve() +// zeroed valid Fourier modes. Both are now evaluated in an exactly rescaled +// frame. +void test_bccb_rank_complex_boundary() { + typedef std::complex<double> Complex; + typedef Matrix<Complex, Dynamic, 1> CVec; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + const double mx = (std::numeric_limits<double>::max)(); + + // n2 = 2, n1 = 1: the symbol is exactly [g0 + g1, g0 - g1], so pick the + // generator from the desired spectrum. |s0| overflows while both of its + // components are finite. + const Complex s0(0.75 * mx, 0.75 * mx), s1(1e300, 0.0); + CMat G(2, 1); + G(0, 0) = (s0 + s1) * 0.5; + G(1, 0) = (s0 - s1) * 0.5; + Bccb<Complex> C(G); + VERIFY_IS_EQUAL(C.rank(), 2); + + // The second Fourier mode must be inverted, not zeroed: a product of a small + // vector solves back to that vector (the accuracy is limited by the condition + // number |s0| / |s1| ~ 2.5e8). + CVec x0(2); + x0[0] = Complex(1e-10, -2e-10); + x0[1] = Complex(-3e-10, 1e-10); + CVec b = C * x0; + VERIFY(b.allFinite()); + CVec x = C.solve(b); + VERIFY(((x - x0).cwiseAbs().maxCoeff() / x0.cwiseAbs().maxCoeff()) <= 1e-6); + + // A genuinely negligible second entry still truncates in the scaled frame. + CMat G2(2, 1); + G2(0, 0) = (s0 + Complex(1)) * 0.5; + G2(1, 0) = (s0 - Complex(1)) * 0.5; + VERIFY_IS_EQUAL(Bccb<Complex>(G2).rank(), 1); +} + +// A single 2^-e frame factor is itself subnormal once the frame exceeds the +// exponent range, and reads as zero under flush-to-zero: every scaled modulus +// and the threshold collapse together, the rank is over-reported, and solve() +// inverts a mode it should have truncated. Hence the two exact factors. +void test_bccb_rank_flush_to_zero() { + ScopedFlushToZero flush_to_zero; + if (!flush_to_zero.isSupported()) return; + + using Complex = std::complex<double>; + using CMat = Matrix<Complex, Dynamic, Dynamic>; + const double mx = (std::numeric_limits<double>::max)(); + + // Spectrum [s0, 1] with |s0| at the overflow boundary: the second mode is + // negligible against it, so the operator is rank one. + const Complex s0(0.75 * mx, 0.75 * mx); + CMat G(2, 1); + G(0, 0) = (s0 + Complex(1)) * 0.5; + G(1, 0) = (s0 - Complex(1)) * 0.5; + VERIFY_IS_EQUAL(Bccb<Complex>(G).rank(), 1); +} + +template <typename Scalar> +void test_bccb_eigen(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> Complex; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + CMat denseC = reference_bccb<Scalar>(G).template cast<Complex>(); + const Index N = n1 * n2; + + Matrix<Complex, Dynamic, 1> lam = C.eigenvalues(); + CMat V = C.eigenvectors(); + VERIFY_IS_APPROX((denseC * V).eval(), (V * lam.asDiagonal()).eval()); + VERIFY_IS_APPROX((V.adjoint() * V).eval(), CMat(CMat::Identity(N, N))); +} + +template <typename Scalar> +void test_bccb_svd(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> Complex; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + Mat G = Mat::Random(n2, n1); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + + Matrix<RealScalar, Dynamic, 1> sv = C.singularValues(); + JacobiSVD<Mat> svd(dense); + VERIFY_IS_APPROX(sv, svd.singularValues()); + + CMat U = C.matrixU(), V = C.matrixV(); + VERIFY_IS_APPROX((U * sv.template cast<Complex>().asDiagonal() * V.adjoint()).eval(), + CMat(dense.template cast<Complex>())); + VERIFY_IS_APPROX((U.adjoint() * U).eval(), CMat(CMat::Identity(N, N))); + VERIFY_IS_APPROX((V.adjoint() * V).eval(), CMat(CMat::Identity(N, N))); +} + +template <typename Scalar> +void test_bccb_inverse(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(n2, n1); + G(0, 0) += Scalar(RealScalar(2 * n1 * n2)); // safely invertible, tame determinant scale + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + const Index N = n1 * n2; + + Mat inv = C.inverse(); + VERIFY_IS_APPROX((inv * dense).eval(), Mat(Mat::Identity(N, N))); + + Vec b = Vec::Random(N); + VERIFY_IS_APPROX((C.inverse() * b).eval(), C.solve(b)); +} + +template <typename Scalar> +void test_bccb_determinant(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Mat G = Mat::Random(n2, n1); + G(0, 0) += Scalar(RealScalar(2 * n1 * n2)); + Bccb<Scalar> C(G); + Mat dense = reference_bccb<Scalar>(G); + VERIFY_IS_APPROX(C.determinant(), dense.determinant()); +} + +template <typename Scalar> +void test_bccb_matrix_free_cg(Index n2, Index n1) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + // Symmetrizing the generator under the two-dimensional index reversal makes + // the BCCB matrix symmetric; diagonal dominance then makes it positive + // definite. + Mat G = Mat::Random(n2, n1); + Mat Gs = G; + for (Index k1 = 0; k1 < n1; ++k1) + for (Index k2 = 0; k2 < n2; ++k2) Gs(k2, k1) = (G(k2, k1) + G((n2 - k2) % n2, (n1 - k1) % n1)) / Scalar(2); + Gs(0, 0) += Scalar(RealScalar(2 * n1 * n2)); + Bccb<Scalar> C(Gs); + Mat dense = reference_bccb<Scalar>(Gs); + VERIFY_IS_APPROX(dense, Mat(dense.transpose())); + + const Index N = n1 * n2; + Vec b = Vec::Random(N); + ConjugateGradient<Bccb<Scalar>, Lower | Upper, IdentityPreconditioner> cg; + cg.compute(C); + Vec x = cg.solve(b); + VERIFY(cg.info() == Success); + VERIFY_IS_APPROX((dense * x).eval(), b); +} + +template <typename Scalar, int N2, int N1> +void test_bccb_fixed() { + typedef Matrix<Scalar, N2, N1> GenMat; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + GenMat G = GenMat::Random(); + Bccb<Scalar, N2, N1> C(G); + STATIC_CHECK((Bccb<Scalar, N2, N1>::RowsAtCompileTime == N1 * N2)); + STATIC_CHECK((internal::remove_all_t<decltype(makeBccb(G))>::RowsAtCompileTime == N1 * N2)); + + Mat dense = reference_bccb<Scalar>(Mat(G)); + Matrix<Scalar, N1 * N2, N1 * N2> Cd = C; + VERIFY_IS_APPROX(Mat(Cd), dense); + + Matrix<Scalar, N1 * N2, 1> x = Matrix<Scalar, N1 * N2, 1>::Random(); + Matrix<Scalar, N1 * N2, 1> y = C * x; + VERIFY_IS_APPROX(y, (dense * x).eval()); +} + +EIGEN_DECLARE_TEST(structured_bccb) { + for (int i = 0; i < g_repeat; ++i) { + // Products, dense assignment, coefficient access: scalar tier (N <= 32) and + // 2-D FFT tier, including single-row/column-of-blocks degenerate shapes. + CALL_SUBTEST_1((test_bccb_product<double>(1, 1))); + CALL_SUBTEST_1((test_bccb_product<double>(3, 4))); + CALL_SUBTEST_1((test_bccb_product<double>(5, 5))); + CALL_SUBTEST_1((test_bccb_product<double>(8, 6))); + CALL_SUBTEST_1((test_bccb_product<double>(7, 5))); // odd prime dimensions + CALL_SUBTEST_1((test_bccb_product<double>(1, 40))); + CALL_SUBTEST_1((test_bccb_product<double>(40, 1))); + CALL_SUBTEST_1((test_bccb_product<float>(6, 8))); + CALL_SUBTEST_1((test_bccb_product<std::complex<double>>(4, 3))); + CALL_SUBTEST_1((test_bccb_product<std::complex<double>>(9, 6))); + CALL_SUBTEST_1((test_bccb_product<std::complex<float>>(6, 7))); + + // Transposition family with exact symbol round trips. + CALL_SUBTEST_2((test_bccb_transpose<double>(3, 4))); + CALL_SUBTEST_2((test_bccb_transpose<double>(8, 6))); + CALL_SUBTEST_2((test_bccb_transpose<double>(1, 12))); + CALL_SUBTEST_2((test_bccb_transpose<std::complex<double>>(6, 8))); + CALL_SUBTEST_2((test_bccb_transpose<std::complex<float>>(5, 7))); + + // Direct and pseudo-inverse solves, degenerate operators. + CALL_SUBTEST_3((test_bccb_solve<double>(1, 1))); + CALL_SUBTEST_3((test_bccb_solve<double>(6, 8))); + CALL_SUBTEST_3((test_bccb_solve<float>(5, 6))); + CALL_SUBTEST_3((test_bccb_solve<std::complex<double>>(6, 6))); + CALL_SUBTEST_3((test_bccb_rank_deficient<std::complex<double>>(6, 8, 4))); + CALL_SUBTEST_3((test_bccb_rank_deficient<std::complex<float>>(4, 5, 2))); + CALL_SUBTEST_3((test_bccb_zero<double>(4, 5))); + CALL_SUBTEST_3(test_bccb_nan_propagation(6, 7)); + + // Closed-form eigendecomposition, SVD, inverse, determinant. + CALL_SUBTEST_4((test_bccb_eigen<double>(4, 5))); + CALL_SUBTEST_4((test_bccb_eigen<double>(1, 8))); + CALL_SUBTEST_4((test_bccb_eigen<std::complex<double>>(3, 5))); + CALL_SUBTEST_4((test_bccb_svd<double>(4, 5))); + CALL_SUBTEST_4((test_bccb_svd<std::complex<double>>(4, 4))); + CALL_SUBTEST_4((test_bccb_svd<float>(3, 4))); + CALL_SUBTEST_4((test_bccb_inverse<double>(6, 8))); + CALL_SUBTEST_4((test_bccb_inverse<std::complex<double>>(5, 5))); + CALL_SUBTEST_4((test_bccb_determinant<double>(3, 4))); + CALL_SUBTEST_4((test_bccb_determinant<std::complex<double>>(3, 3))); + + // Matrix-free iterative solve and fixed-size generators. + CALL_SUBTEST_5((test_bccb_matrix_free_cg<double>(8, 10))); + CALL_SUBTEST_5((test_bccb_fixed<double, 3, 4>())); + CALL_SUBTEST_5((test_bccb_fixed<std::complex<float>, 4, 3>())); + + // Non-5-smooth (prime) block dimensions: the padded product-embedding grid, + // in each axis combination (n2 awkward, n1 awkward, both). + CALL_SUBTEST_5((test_bccb_prime_dimensions<double>(7, 6))); + CALL_SUBTEST_5((test_bccb_prime_dimensions<double>(6, 7))); + CALL_SUBTEST_5((test_bccb_prime_dimensions<double>(7, 7))); + CALL_SUBTEST_5((test_bccb_prime_dimensions<std::complex<double>>(11, 7))); + CALL_SUBTEST_5((test_bccb_prime_dimensions<float>(7, 6))); + + // Product lifetime, aliasing, mixed-scalar promotion, dimension mismatches, + // across the scalar-loop and FFT dispatch tiers. + CALL_SUBTEST_6((test_bccb_delayed_product<double>(4, 5))); + CALL_SUBTEST_6((test_bccb_delayed_product<std::complex<double>>(6, 8))); + CALL_SUBTEST_6((test_bccb_aliased_product<double>(3, 4))); + CALL_SUBTEST_6((test_bccb_aliased_product<double>(8, 6))); + CALL_SUBTEST_6((test_bccb_aliased_product<std::complex<double>>(6, 8))); + CALL_SUBTEST_6((test_bccb_aliased_expression<double>(3, 4))); + CALL_SUBTEST_6((test_bccb_aliased_expression<double>(6, 8))); + CALL_SUBTEST_6((test_bccb_aliased_expression<std::complex<double>>(6, 8))); + CALL_SUBTEST_6((test_bccb_mixed_scalar<double>(3, 4))); + CALL_SUBTEST_6((test_bccb_mixed_scalar<double>(6, 8))); + CALL_SUBTEST_6((test_bccb_mixed_scalar<float>(8, 8))); + CALL_SUBTEST_6(test_bccb_dimension_asserts()); + + // Finite-range robustness: scaled transforms, the complex overflow boundary, + // balanced determinant accumulation, and the rank threshold boundary. + CALL_SUBTEST_7(test_bccb_finite_overflow()); + CALL_SUBTEST_7((test_bccb_fft_complex_boundary<double>(6, 8))); + CALL_SUBTEST_7((test_bccb_fft_complex_boundary<float>(6, 8))); + CALL_SUBTEST_7(test_bccb_determinant_scaled()); + CALL_SUBTEST_7(test_bccb_rank_boundaries()); + CALL_SUBTEST_7(test_bccb_rank_complex_boundary()); + CALL_SUBTEST_7(test_bccb_rank_flush_to_zero()); + + // Entrywise Inf/NaN propagation: FFT-sized operators must fall back to the + // direct kernel; small ones are IEEE-exact already. A zero right-hand side + // under a non-finite operator yields NaN, not zero. + CALL_SUBTEST_7((test_bccb_nonfinite_product<double>(6, 8))); + CALL_SUBTEST_7((test_bccb_nonfinite_product<double>(3, 4))); + CALL_SUBTEST_7((test_bccb_nonfinite_product<std::complex<double>>(6, 8))); + CALL_SUBTEST_7(test_bccb_nonfinite_zero_rhs(6, 8)); + } +}