TRSM: Deepen the k-blocks and panel the right-hand side once the operand outgrows the cache libeigen/eigen!3053 Closes #3162 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 96efe53..f75365e 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h
@@ -385,6 +385,9 @@ } // WARNING Below, we assume that Traits::nr is a power of two. Index nc = numext::mini<Index>(rhs_panel_budget / (2 * k * sizeof(RhsScalar)), max_nc) & (~(Traits::nr - 1)); + // As in the threaded branch, a budget below one kernel-width panel still takes nr columns: a zero + // width would divide by zero below. + nc = numext::maxi<Index>(nc, Traits::nr); if (n > nc) { // We are really blocking over the columns: // -> reduce blocking size to make sure the last block is as large as possible
diff --git a/Eigen/src/Core/products/TriangularSolverMatrix.h b/Eigen/src/Core/products/TriangularSolverMatrix.h index f720f73..d939411 100644 --- a/Eigen/src/Core/products/TriangularSolverMatrix.h +++ b/Eigen/src/Core/products/TriangularSolverMatrix.h
@@ -186,6 +186,51 @@ } }; +/** \internal Entries of the solved operand a blocked triangular solve may keep in cache while its + * k-blocks sweep it: a quarter of the L3, which on a multi-die part is the package total of which one + * core reaches a fraction, and never less than the L2. */ +template <typename Scalar> +std::ptrdiff_t triangular_solve_budget(std::ptrdiff_t l2, std::ptrdiff_t l3) { + return (numext::maxi)(l3 / 4, l2) / std::ptrdiff_t(sizeof(Scalar)); +} + +/** \internal Columns of the right-hand side a solve on the left takes per panel: the widest multiple + * of nr whose size x nc entries fit the budget, or all cols. Every panel re-packs the triangle, + * size^2/2 copies against size^2*nc/2 multiply-adds, and panels narrower than 512 columns measured + * slower than none, so a right-hand side that would need them is solved whole. */ +template <typename Index> +Index triangular_solve_panel_columns(Index size, Index cols, std::ptrdiff_t budget, Index nr) { + eigen_internal_assert(size > 0); + const std::ptrdiff_t width = numext::round_down<std::ptrdiff_t>(budget / size, nr); + return width >= 512 && width < cols ? Index(width) : cols; +} + +/** \internal Depth of the k-blocks of a solve against otherSize columns (on the left) or rows (on the + * right), whose packed buffers hold kc x extent entries. The KcFactor 4 depth SolveTriangular.h blocks + * with keeps the share of flops in the kc x kc diagonal blocks, which run below gebp's rate, near + * kc/size. Every k-block also sweeps the operand and packs a slab of the triangle, costs that scale as + * 1/kc and outweigh that share once the operand exceeds the budget or, for a slab packed one column + * run at a time (slabRuns), once half the triangle does. The depth then grows toward the + * single-threaded GEMM depth, no further than size/8 and 160, where deeper diagonal blocks measured + * slower than the sweeps they save (by up to 10% at the GEMM depth in AVX2 builds), and, when the + * buffers fit on the stack at the blocking's depth, no further than the depth at which they still do. + * A caller that preallocated the buffers sized them for the blocking's depth. */ +template <typename Scalar, typename Index> +Index triangular_solve_kc(Index size, Index otherSize, Index extent, std::ptrdiff_t budget, bool slabRuns, + level3_blocking<Scalar, Scalar>& blocking) { + const bool deep = std::ptrdiff_t(size) * otherSize > budget || (slabRuns && std::ptrdiff_t(size) * size / 2 > budget); + if (!deep || blocking.blockA() != nullptr) return blocking.kc(); + Index kc = size, mc = size, nc = otherSize; + computeProductBlockingSizes<Scalar, Scalar>(kc, mc, nc); + kc = (numext::mini)(kc, numext::round_down((numext::mini)(size / 8, Index(160)), Index(8))); +#if defined(EIGEN_ALLOCA) && !defined(EIGEN_NO_ALLOCA) + const std::ptrdiff_t stackKc = + std::ptrdiff_t(EIGEN_STACK_ALLOCATION_LIMIT) / (std::ptrdiff_t(sizeof(Scalar)) * extent); + if (blocking.kc() <= stackKc) kc = Index((numext::mini)(std::ptrdiff_t(kc), stackKc)); +#endif + return (numext::maxi)(kc, blocking.kc()); +} + /* Optimized triangular solver with multiple right hand side and the triangular matrix on the left */ template <typename Scalar, typename Index, int Mode, bool Conjugate, int TriStorageOrder, int OtherInnerStride> @@ -200,8 +245,6 @@ Index triStride, Scalar* _other, Index otherIncr, Index otherStride, level3_blocking<Scalar, Scalar>& blocking) { - Index cols = otherSize; - std::ptrdiff_t l1, l2, l3; manage_caching_sizes(GetAction, &l1, &l2, &l3); @@ -213,9 +256,9 @@ // For small problem sizes trsmKernel compiled with clang is generally faster. // TODO: Investigate better heuristics for cutoffs. double L2Cap = 0.5; // 50% of L2 size - if (size < avx512_trsm_cutoff<Scalar>(l2, cols, L2Cap)) { + if (size < avx512_trsm_cutoff<Scalar>(l2, otherSize, L2Cap)) { trsmKernelL<Scalar, Index, Mode, Conjugate, TriStorageOrder, 1, /*Specialized=*/true>::kernel( - size, cols, _tri, triStride, _other, 1, otherStride); + size, otherSize, _tri, triStride, _other, 1, otherStride); return; } } @@ -224,17 +267,26 @@ using TriMapper = const_blas_data_mapper<Scalar, Index, TriStorageOrder>; using OtherMapper = blas_data_mapper<Scalar, Index, ColMajor, Unaligned, OtherInnerStride>; TriMapper tri(_tri, triStride); - OtherMapper other(_other, otherStride, otherIncr); using Traits = gebp_traits<Scalar, Scalar>; enum { SmallPanelWidth = plain_enum_max(Traits::mr, Traits::nr), IsLower = (Mode & Lower) == Lower }; - Index kc = blocking.kc(); // cache block size along the K direction - Index mc = (std::min)(size, blocking.mc()); // cache block size along the M direction + // Every k-block updates the rows of the right-hand side beyond it through gebp, so a right-hand side + // solved whole streams through the caches size/kc times (issue #3162); column panels that stay in + // cache keep those sweeps out of memory. This kernel packs the slabs of the triangle by rows, which + // costs more the deeper they are, so a large triangle alone does not deepen it. + const std::ptrdiff_t budget = triangular_solve_budget<Scalar>(l2, l3); + const Index nc = triangular_solve_panel_columns(size, otherSize, budget, Index(Traits::nr)); + const Index mc = (numext::mini)(size, blocking.mc()); // cache block size along the M direction + // The tr solve below packs up to SmallPanelWidth x kc entries of the triangle into blockA. + const Index blockARows = (numext::maxi)(mc, Index(SmallPanelWidth)); + // cache block size along the K direction + const Index kc = triangular_solve_kc<Scalar>(size, otherSize, (numext::maxi)(blockARows, nc), budget, + /*slabRuns=*/false, blocking); - std::size_t sizeA = kc * mc; - std::size_t sizeB = kc * cols; + std::size_t sizeA = kc * blockARows; + std::size_t sizeB = kc * nc; ei_declare_aligned_stack_constructed_variable(Scalar, blockA, sizeA, blocking.blockA()); ei_declare_aligned_stack_constructed_variable(Scalar, blockB, sizeB, blocking.blockB()); @@ -247,74 +299,82 @@ // the goal here is to subdivide the Rhs panels such that we keep some cache // coherence when accessing the rhs elements - Index subcols = cols > 0 ? l2 / (4 * sizeof(Scalar) * std::max<Index>(otherStride, size)) : 0; - subcols = std::max<Index>((subcols / Traits::nr) * Traits::nr, Traits::nr); + Index subcols = otherSize > 0 ? l2 / (4 * sizeof(Scalar) * numext::maxi<Index>(otherStride, size)) : 0; + subcols = numext::maxi<Index>((subcols / Traits::nr) * Traits::nr, Traits::nr); - for (Index k2 = IsLower ? 0 : size; IsLower ? k2 < size : k2 > 0; IsLower ? k2 += kc : k2 -= kc) { - const Index actual_kc = (std::min)(IsLower ? size - k2 : k2, kc); + for (Index j0 = 0; j0 < otherSize; j0 += nc) { + const Index cols = (numext::mini)(otherSize - j0, nc); + Scalar* _panel = _other + j0 * otherStride; + OtherMapper other(_panel, otherStride, otherIncr); - // We have selected and packed a big horizontal panel R1 of rhs. Let B be the packed copy of this panel, - // and R2 the remaining part of rhs. The corresponding vertical panel of lhs is split into - // A11 (the triangular part) and A21 the remaining rectangular part. - // Then the high level algorithm is: - // - B = R1 => general block copy (done during the next step) - // - R1 = A11^-1 B => tricky part - // - update B from the new R1 => actually this has to be performed continuously during the above step - // - R2 -= A21 * B => GEPP + for (Index k2 = IsLower ? 0 : size; IsLower ? k2 < size : k2 > 0; IsLower ? k2 += kc : k2 -= kc) { + const Index actual_kc = (numext::mini)(IsLower ? size - k2 : k2, kc); - // The tricky part: compute R1 = A11^-1 B while updating B from R1 - // The idea is to split A11 into multiple small vertical panels. - // Each panel can be split into a small triangular part T1k which is processed without optimization, - // and the remaining small part T2k which is processed using gebp with appropriate block strides - for (Index j2 = 0; j2 < cols; j2 += subcols) { - Index actual_cols = (std::min)(cols - j2, subcols); - // for each small vertical panels [T1k^T, T2k^T]^T of lhs - for (Index k1 = 0; k1 < actual_kc; k1 += SmallPanelWidth) { - Index actualPanelWidth = std::min<Index>(actual_kc - k1, SmallPanelWidth); - // tr solve - { - Index i = IsLower ? k2 + k1 : k2 - k1 - 1; + // We have selected and packed a big horizontal panel R1 of rhs. Let B be the packed copy of this panel, + // and R2 the remaining part of rhs. The corresponding vertical panel of lhs is split into + // A11 (the triangular part) and A21 the remaining rectangular part. + // Then the high level algorithm is: + // - B = R1 => general block copy (done during the next step) + // - R1 = A11^-1 B => tricky part + // - update B from the new R1 => actually this has to be performed continuously during the above step + // - R2 -= A21 * B => GEPP + + // The tricky part: compute R1 = A11^-1 B while updating B from R1 + // The idea is to split A11 into multiple small vertical panels. + // Each panel can be split into a small triangular part T1k which is processed without optimization, + // and the remaining small part T2k which is processed using gebp with appropriate block strides + for (Index j2 = 0; j2 < cols; j2 += subcols) { + Index actual_cols = (numext::mini)(cols - j2, subcols); + // for each small vertical panels [T1k^T, T2k^T]^T of lhs + for (Index k1 = 0; k1 < actual_kc; k1 += SmallPanelWidth) { + Index actualPanelWidth = numext::mini<Index>(actual_kc - k1, SmallPanelWidth); + // tr solve + { + Index i = IsLower ? k2 + k1 : k2 - k1 - 1; #if defined(EIGEN_VECTORIZE_AVX512) && defined(EIGEN_USE_AVX512_TRSM_L_KERNELS) && EIGEN_USE_AVX512_TRSM_L_KERNELS - EIGEN_IF_CONSTEXPR ((OtherInnerStride == 1 && - (std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value))) { - i = IsLower ? k2 + k1 : k2 - k1 - actualPanelWidth; - } + EIGEN_IF_CONSTEXPR ((OtherInnerStride == 1 && + (std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value))) { + i = IsLower ? k2 + k1 : k2 - k1 - actualPanelWidth; + } #endif - trsmKernelL<Scalar, Index, Mode, Conjugate, TriStorageOrder, OtherInnerStride, /*Specialized=*/true>::kernel( - actualPanelWidth, actual_cols, _tri + i + (i)*triStride, triStride, - _other + i * otherIncr + j2 * otherStride, otherIncr, otherStride); - } + trsmKernelL<Scalar, Index, Mode, Conjugate, TriStorageOrder, OtherInnerStride, + /*Specialized=*/true>::kernel(actualPanelWidth, actual_cols, _tri + i + (i)*triStride, + triStride, _panel + i * otherIncr + j2 * otherStride, otherIncr, + otherStride); + } - Index lengthTarget = actual_kc - k1 - actualPanelWidth; - Index startBlock = IsLower ? k2 + k1 : k2 - k1 - actualPanelWidth; - Index blockBOffset = IsLower ? k1 : lengthTarget; + Index lengthTarget = actual_kc - k1 - actualPanelWidth; + Index startBlock = IsLower ? k2 + k1 : k2 - k1 - actualPanelWidth; + Index blockBOffset = IsLower ? k1 : lengthTarget; - // update the respective rows of B from other - pack_rhs(blockB + actual_kc * j2, other.getSubMapper(startBlock, j2), actualPanelWidth, actual_cols, actual_kc, - blockBOffset); + // update the respective rows of B from other + pack_rhs(blockB + actual_kc * j2, other.getSubMapper(startBlock, j2), actualPanelWidth, actual_cols, + actual_kc, blockBOffset); - // GEBP - if (lengthTarget > 0) { - Index startTarget = IsLower ? k2 + k1 + actualPanelWidth : k2 - actual_kc; + // GEBP + if (lengthTarget > 0) { + Index startTarget = IsLower ? k2 + k1 + actualPanelWidth : k2 - actual_kc; - pack_lhs(blockA, tri.getSubMapper(startTarget, startBlock), actualPanelWidth, lengthTarget); + pack_lhs(blockA, tri.getSubMapper(startTarget, startBlock), actualPanelWidth, lengthTarget); - gebp_kernel(other.getSubMapper(startTarget, j2), blockA, blockB + actual_kc * j2, lengthTarget, - actualPanelWidth, actual_cols, Scalar(-1), actualPanelWidth, actual_kc, 0, blockBOffset); + gebp_kernel(other.getSubMapper(startTarget, j2), blockA, blockB + actual_kc * j2, lengthTarget, + actualPanelWidth, actual_cols, Scalar(-1), actualPanelWidth, actual_kc, 0, blockBOffset); + } } } - } - // R2 -= A21 * B => GEPP - { - Index start = IsLower ? k2 + kc : 0; - Index end = IsLower ? size : k2 - kc; - for (Index i2 = start; i2 < end; i2 += mc) { - const Index actual_mc = (std::min)(mc, end - i2); - if (actual_mc > 0) { - pack_lhs(blockA, tri.getSubMapper(i2, IsLower ? k2 : k2 - kc), actual_kc, actual_mc); + // R2 -= A21 * B => GEPP + { + Index start = IsLower ? k2 + kc : 0; + Index end = IsLower ? size : k2 - kc; + for (Index i2 = start; i2 < end; i2 += mc) { + const Index actual_mc = (numext::mini)(mc, end - i2); + if (actual_mc > 0) { + pack_lhs(blockA, tri.getSubMapper(i2, IsLower ? k2 : k2 - kc), actual_kc, actual_mc); - gebp_kernel(other.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc, cols, Scalar(-1), -1, -1, 0, 0); + gebp_kernel(other.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc, cols, Scalar(-1), -1, -1, 0, + 0); + } } } } @@ -338,13 +398,14 @@ level3_blocking<Scalar, Scalar>& blocking) { Index rows = otherSize; + std::ptrdiff_t l1, l2, l3; + manage_caching_sizes(GetAction, &l1, &l2, &l3); + #if defined(EIGEN_VECTORIZE_AVX512) && defined(EIGEN_USE_AVX512_TRSM_R_KERNELS) && EIGEN_USE_AVX512_TRSM_R_KERNELS && \ EIGEN_ENABLE_AVX512_NOCOPY_TRSM_R_CUTOFFS EIGEN_IF_CONSTEXPR ((OtherInnerStride == 1 && (std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value))) { // TODO: Investigate better heuristics for cutoffs. - std::ptrdiff_t l1, l2, l3; - manage_caching_sizes(GetAction, &l1, &l2, &l3); double L2Cap = 0.5; // 50% of L2 size if (size < avx512_trsm_cutoff<Scalar>(l2, rows, L2Cap)) { trsmKernelR<Scalar, Index, Mode, Conjugate, TriStorageOrder, OtherInnerStride, /*Specialized=*/true>::kernel( @@ -366,8 +427,21 @@ IsLower = (Mode & Lower) == Lower }; - Index kc = blocking.kc(); // cache block size along the K direction - Index mc = (std::min)(rows, blocking.mc()); // cache block size along the M direction + // Every k-block sweeps all rows of the left-hand side through gebp and packs its slab of a + // column-major triangle one column run at a time, so a large triangle alone deepens this kernel + // (issue #3162); a row-major triangle is packed by rows, as on the left. + const std::ptrdiff_t budget = triangular_solve_budget<Scalar>(l2, l3); + Index mc = (numext::mini)(rows, blocking.mc()); // cache block size along the M direction + // cache block size along the K direction + const Index kc = triangular_solve_kc<Scalar>(size, rows, (numext::maxi)(mc, size), budget, + /*slabRuns=*/TriStorageOrder == ColMajor, blocking); + // blockA packs kc x mc entries of the left-hand side, and rows can far exceed size. Past half the + // budget a deeper kc takes proportionally fewer rows per pass, so blockA does not outgrow the buffer + // the blocking chose. + if (kc > blocking.kc()) { + const std::ptrdiff_t maxA = (numext::maxi)(std::ptrdiff_t(blocking.kc()) * mc, budget / 2); + mc = (numext::mini)(mc, (numext::maxi)(Index(Traits::mr), numext::round_down(Index(maxA / kc), Index(Traits::mr)))); + } std::size_t sizeA = kc * mc; std::size_t sizeB = kc * size; @@ -383,7 +457,7 @@ pack_lhs_panel; for (Index k2 = IsLower ? size : 0; IsLower ? k2 > 0 : k2 < size; IsLower ? k2 -= kc : k2 += kc) { - const Index actual_kc = (std::min)(IsLower ? k2 : size - k2, kc); + const Index actual_kc = (numext::mini)(IsLower ? k2 : size - k2, kc); Index actual_k2 = IsLower ? k2 - actual_kc : k2; Index startPanel = IsLower ? 0 : k2 + actual_kc; @@ -396,7 +470,7 @@ // neglecting the blocks overlapping the diagonal { for (Index j2 = 0; j2 < actual_kc; j2 += SmallPanelWidth) { - Index actualPanelWidth = std::min<Index>(actual_kc - j2, SmallPanelWidth); + Index actualPanelWidth = numext::mini<Index>(actual_kc - j2, SmallPanelWidth); Index actual_j2 = actual_k2 + j2; Index panelOffset = IsLower ? j2 + actualPanelWidth : 0; Index panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2; @@ -408,7 +482,7 @@ } for (Index i2 = 0; i2 < rows; i2 += mc) { - const Index actual_mc = (std::min)(mc, rows - i2); + const Index actual_mc = (numext::mini)(mc, rows - i2); // triangular solver kernel { @@ -417,7 +491,7 @@ : Index(SmallPanelWidth))) : 0; IsLower ? j2 >= 0 : j2 < actual_kc; IsLower ? j2 -= SmallPanelWidth : j2 += SmallPanelWidth) { - Index actualPanelWidth = std::min<Index>(actual_kc - j2, SmallPanelWidth); + Index actualPanelWidth = numext::mini<Index>(actual_kc - j2, SmallPanelWidth); Index absolute_j2 = actual_k2 + j2; Index panelOffset = IsLower ? j2 + actualPanelWidth : 0; Index panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2;
diff --git a/benchmarks/Core/bench_trsm.cpp b/benchmarks/Core/bench_trsm.cpp index f7f7582..4ab6b81 100644 --- a/benchmarks/Core/bench_trsm.cpp +++ b/benchmarks/Core/bench_trsm.cpp
@@ -75,9 +75,15 @@ // clang-format off BENCHMARK(BM_TRSM_Left<float, Lower>)->ArgsProduct({{64, 256, 512}, {1, 16, 64}})->Name("TRSM_Left_float_Lower"); BENCHMARK(BM_TRSM_Left<double, Lower>)->ArgsProduct({{64, 256, 512}, {1, 16, 64}})->Name("TRSM_Left_double_Lower"); +// Operands from below to well beyond the last-level cache (issue #3162). Clang AVX-512 builds on hosts +// with an L2 of 2 MiB or more solve the float ones of size 1024 with the unblocked kernel. +BENCHMARK(BM_TRSM_Left<float, Lower>)->ArgsProduct({{1024, 2048, 4096}, {1024, 4096}})->Name("TRSM_Left_float_Lower_Large"); +BENCHMARK(BM_TRSM_Left<double, Lower>)->ArgsProduct({{1024, 2048, 4096}, {1024, 4096}})->Name("TRSM_Left_double_Lower_Large"); // ---------- TRSM Right benchmarks ---------- BENCHMARK(BM_TRSM_Right<float, Lower>)->ArgsProduct({{64, 256, 512}, {1, 16, 64}})->Name("TRSM_Right_float_Lower"); BENCHMARK(BM_TRSM_Right<double, Lower>)->ArgsProduct({{64, 256, 512}, {1, 16, 64}})->Name("TRSM_Right_double_Lower"); +BENCHMARK(BM_TRSM_Right<float, Lower>)->ArgsProduct({{1024, 2048, 4096}, {1024, 4096}})->Name("TRSM_Right_float_Lower_Large"); +BENCHMARK(BM_TRSM_Right<double, Lower>)->ArgsProduct({{1024, 2048, 4096}, {1024, 4096}})->Name("TRSM_Right_double_Lower_Large"); // clang-format on
diff --git a/test/product_trsolve.cpp b/test/product_trsolve.cpp index 7f95e7f..ea7b7a1 100644 --- a/test/product_trsolve.cpp +++ b/test/product_trsolve.cpp
@@ -360,6 +360,95 @@ } } +// The blocked solvers take their k-block depth, and a solve on the left its column panels, from the +// cache sizes (issue #3162). Both take effect only for operands larger than the cache, which test-sized +// solves are not on most hosts, so the caches are set here for a budget of 520 columns of the operand: +// 1100 right-hand sides then span two full panels and a partial one, and on x86 the depth rises above +// the blocking's on either side. Both storage orders of the right-hand side reach both kernels from +// either side, since the dispatcher transposes a row-major operand. +template <typename Scalar> +void trsolve_panels(int size, int cols, bool allCases) { + using RealScalar = typename NumTraits<Scalar>::Real; + using MatrixX = Matrix<Scalar, Dynamic, Dynamic>; + using RowMatrixX = Matrix<Scalar, Dynamic, Dynamic, RowMajor>; + + const std::ptrdiff_t l1 = 32768, l2 = 65536; + const std::ptrdiff_t panelL3 = 4 * 520 * std::ptrdiff_t(size) * std::ptrdiff_t(sizeof(Scalar)); + setCpuCacheSizes(l1, l2, panelL3); + const std::ptrdiff_t budget = internal::triangular_solve_budget<Scalar>(l2, panelL3); + // Reach: the budget holds fewer columns than the right-hand side has, and the operand exceeds it, + // which deepens the k-blocks wherever the blocking leaves room. + VERIFY(internal::triangular_solve_panel_columns(Index(size), Index(cols), budget, + Index(internal::gebp_traits<Scalar, Scalar>::nr)) < cols); + VERIFY(std::ptrdiff_t(size) * cols > budget); + + // Scaling the off-diagonal part by 1/size keeps the triangles well conditioned whether the diagonal + // is used or, for the unit views, taken as 1. + MatrixX cmLhs = MatrixX::Random(size, size) / RealScalar(size); + cmLhs.diagonal().array() += RealScalar(1); + RowMatrixX rmLhs = cmLhs; + MatrixX cmRhs(size, cols); + RowMatrixX rmRhs(size, cols); + MatrixX ref(size, cols); + + VERIFY_TRSM(cmLhs.template triangularView<Lower>(), cmRhs); + VERIFY_TRSM(cmLhs.template triangularView<Upper>(), cmRhs); + VERIFY_TRSM(cmLhs.template triangularView<Lower>(), rmRhs); + VERIFY_TRSM(cmLhs.template triangularView<Upper>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.template triangularView<Lower>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.template triangularView<Upper>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.template triangularView<Lower>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.template triangularView<Upper>(), rmRhs); + if (!allCases) return; + + VERIFY_TRSM(cmLhs.template triangularView<UnitLower>(), cmRhs); + VERIFY_TRSM(cmLhs.adjoint().template triangularView<Upper>(), rmRhs); + VERIFY_TRSM(rmLhs.template triangularView<Lower>(), cmRhs); + VERIFY_TRSM(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs); + + VERIFY_TRSM_ONTHERIGHT(cmLhs.template triangularView<UnitUpper>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(cmLhs.adjoint().template triangularView<Lower>(), cmRhs); + VERIFY_TRSM_ONTHERIGHT(rmLhs.template triangularView<Lower>(), rmRhs); + VERIFY_TRSM_ONTHERIGHT(rmLhs.conjugate().template triangularView<UnitUpper>(), cmRhs); + + // A runtime inner stride, by which the panel origins must not be scaled. + MatrixX buffer(2 * size, 2 * cols); + Map<MatrixX, 0, Stride<Dynamic, Dynamic> > cmMap(buffer.data(), size, cols, Stride<Dynamic, Dynamic>(2 * size, 2)); + Map<RowMatrixX, 0, Stride<Dynamic, Dynamic> > rmMap(buffer.data(), size, cols, Stride<Dynamic, Dynamic>(2 * cols, 2)); + buffer.setZero(); + VERIFY_TRSM(cmLhs.template triangularView<Lower>(), cmMap); + buffer.setZero(); + VERIFY_TRSM(cmLhs.template triangularView<Upper>(), rmMap); + buffer.setZero(); + VERIFY_TRSM_ONTHERIGHT(cmLhs.template triangularView<Upper>(), cmMap); + buffer.setZero(); + VERIFY_TRSM_ONTHERIGHT(cmLhs.template triangularView<Lower>(), rmMap); + + // A few rows beside a column-major triangle half of which exceeds the budget: only the triangle makes + // this solve on the right deep. + const std::ptrdiff_t triangleL3 = std::ptrdiff_t(size) * size * std::ptrdiff_t(sizeof(Scalar)); + setCpuCacheSizes(l1, l2, triangleL3); + const std::ptrdiff_t triangleBudget = internal::triangular_solve_budget<Scalar>(l2, triangleL3); + VERIFY(std::ptrdiff_t(size) * size / 2 > triangleBudget && std::ptrdiff_t(size) * 7 <= triangleBudget); + RowMatrixX fewRows(size, 7); + VERIFY_TRSM_ONTHERIGHT(rmLhs.template triangularView<Lower>(), fewRows); + VERIFY_TRSM_ONTHERIGHT(rmLhs.template triangularView<Upper>(), fewRows); +} + +template <int> +void trsolve_panels_all() { + std::ptrdiff_t l1, l2, l3, l3_per_cpu; + internal::manage_caching_sizes(GetAction, &l1, &l2, &l3, &l3_per_cpu); + // 193 raises the depth of every x86 build except SSE2 float and double, which 391 raises; both leave a + // partial last k-block at every depth they reach. + trsolve_panels<float>(193, 1100, true); + trsolve_panels<double>(193, 1100, true); + trsolve_panels<std::complex<double> >(193, 1100, true); + trsolve_panels<float>(391, 1100, false); + trsolve_panels<double>(391, 1100, false); + internal::manage_caching_sizes(SetAction, &l1, &l2, &l3, &l3_per_cpu); +} + EIGEN_DECLARE_TEST(product_trsolve) { for (int i = 0; i < g_repeat; i++) { // matrices @@ -391,4 +480,5 @@ CALL_SUBTEST_15(trsolve_strided_boundary<0>()); CALL_SUBTEST_16(trsolve_indexed_view()); CALL_SUBTEST_17(trsolve_no_malloc_all<0>()); + CALL_SUBTEST_18(trsolve_panels_all<0>()); }