Align temporary array in TensorSelectOp packet evaluator.

Fix MKL enum conversion warning.
### Description  Fix MKL enum conversion warning. ### Reference issue  Fixes #2999. Closes #2999 See merge request libeigen/eigen!2061
2026-04-10 11:34:33 +08:00 · 2025-11-05 19:44:47 +00:00 · 2025-11-05 18:03:22 +00:00 · 2025-11-05 17:12:33 +00:00 · 2025-11-04 22:26:27 +00:00 · 2025-11-04 21:42:32 +00:00
52 changed files with 1013 additions and 601 deletions
--- a/.gitlab/issue_templates/Bug
+++ b/.gitlab/issue_templates/Bug
@@ -1,42 +1,37 @@
 <!--
-Please read this!
+Thank you for submitting an issue!

-Before opening a new issue, make sure to search for keywords in the issues
-filtered by "bug::confirmed" or "bug::unconfirmed" and "bugzilla" label:
-
- https://gitlab.com/libeigen/eigen/-/issues?scope=all&utf8=%E2%9C%93&state=opened&label_name[]=bug%3A%3Aconfirmed
- https://gitlab.com/libeigen/eigen/-/issues?scope=all&utf8=%E2%9C%93&state=opened&label_name[]=bug%3A%3Aunconfirmed
- https://gitlab.com/libeigen/eigen/-/issues?scope=all&utf8=%E2%9C%93&state=opened&label_name[]=bugzilla
-
-and verify the issue you're about to submit isn't a duplicate. -->
+Before opening a new issue, please search for keywords in the existing [list of issues](https://gitlab.com/libeigen/eigen/-/issues?state=opened) to verify it isn't a duplicate.
+-->

 ### Summary
 <!-- Summarize the bug encountered concisely. -->

 ### Environment
-<!-- Please provide your development environment here -->
+<!-- Please provide your development environment. -->
 - **Operating System** : Windows/Linux
 - **Architecture** : x64/Arm64/PowerPC ...
- **Eigen Version** : 3.3.9
- **Compiler Version** : Gcc7.0
+- **Eigen Version** : 5.0.0
+- **Compiler Version** : gcc-12.0
 - **Compile Flags** : -O3 -march=native
 - **Vector Extension** : SSE/AVX/NEON ...

 ### Minimal Example
-<!-- If possible, please create a minimal example here that exhibits the problematic behavior.
-You can also link to [godbolt](https://godbolt.org). But please note that you need to click 
-the "Share" button in the top right-hand corner of the godbolt page where you reproduce the sample 
-code to get the share link instead of in your browser address bar. 
+<!--
+Please create a minimal reproducing example here that exhibits the problematic behavior.
+The example should be complete, in that it can fully build and run.  See the [the guidelines on stackoverflow](https://stackoverflow.com/help/minimal-reproducible-example) for how to create a good minimal example.

-You can read [the guidelines on stackoverflow](https://stackoverflow.com/help/minimal-reproducible-example)
-on how to create a good minimal example. -->
+You can also link to [godbolt](https://godbolt.org). Note that you need to click 
+the "Share" button in the top right-hand corner of the godbolt page to get the share link
+instead of the URL in your browser address bar. 
+-->

 ```cpp
-//show your code here
+// Insert  your code here.
 ```

-### Steps to reproduce
-<!-- Describe how one can reproduce the issue - this is very important. Please use an ordered list. -->
+### Steps to reproduce the issue
+<!-- Describe the necessary steps to reproduce the issue. -->

 1. first step
 2. second step
@@ -49,21 +44,16 @@ on how to create a good minimal example. -->
 <!-- Describe what you should see instead. -->

 ### Relevant logs
-<!-- Add relevant code snippets or program output within blocks marked by " ``` " -->
+<!-- Add relevant build logs or program output within blocks marked by " ``` " -->

-<!-- OPTIONAL: remove this section if you are not reporting a compilation warning issue.-->
-### Warning Messages
-<!-- Show us the warning messages you got! -->
-
-<!-- OPTIONAL: remove this section if you are not reporting a performance issue. -->
-### Benchmark scripts and results
+### [Optional] Benchmark scripts and results
 <!-- Please share any benchmark scripts - either standalone, or using [Google Benchmark](https://github.com/google/benchmark). -->

 ### Anything else that might help
-<!-- It will be better to provide us more information to help narrow down the cause. 
+<!--
+It will be better to provide us more information to help narrow down the cause. 
 Including but not limited to the following: 
 - lines of code that might help us diagnose the problem. 
 - potential ways to address the issue.
- last known working/first broken version (release number or commit hash). --> 
-
- [ ] Have a plan to fix this issue.
+- last known working/first broken version (release number or commit hash).
+--> 
--- a/.gitlab/issue_templates/Feature
+++ b/.gitlab/issue_templates/Feature
@@ -1,6 +1,13 @@
+<!--
+Thank you for submitting a Feature Request!
+
+If you want to run ideas by the maintainers and the Eigen community first, 
+you can chat about them on the [Eigen Discord server](https://discord.gg/2SkEJGqZjR).
+-->
+
 ### Describe the feature you would like to be implemented.

-### Would such a feature be useful for other users? Why?
+### Why Would such a feature be useful for other users?

 ### Any hints on how to implement the requested feature?

--- a/.gitlab/merge_request_templates/Default.md
+++ b/.gitlab/merge_request_templates/Default.md
@@ -0,0 +1,30 @@
+<!--
+Thanks for contributing a merge request!
+
+We recommend that first-time contributors read our [contribution guidelines](https://eigen.tuxfamily.org/index.php?title=Contributing_to_Eigen).
+
+Before submitting the MR, please complete the following checks:
+- Create one PR per feature or bugfix,
+- Run the test suite to verify your changes.
+  See our [test guidelines](https://eigen.tuxfamily.org/index.php?title=Tests).
+- Add tests to cover the bug addressed or any new feature.
+- Document new features.  If it is a substantial change, add it to the [Changelog](https://gitlab.com/libeigen/eigen/-/blob/master/CHANGELOG.md).
+- Leave the following box checked when submitting: `Allow commits from members who can merge to the target branch`.
+  This allows us to rebase and merge your change.
+
+Note that we are a team of volunteers; we appreciate your patience during the review process.
+-->
+
+### Description
+<!--Please explain your changes.-->
+
+%{first_multiline_commit}
+
+### Reference issue
+<!--
+You can link to a specific issue using the gitlab syntax #<issue number>.
+If the MR fixes an issue, write "Fixes #<issue number>" to have the issue automatically closed on merge.
+-->
+
+### Additional information
+<!--Any additional information you think is important.-->
--- a/.gitlab/merge_request_templates/Merge
+++ b/.gitlab/merge_request_templates/Merge
@@ -1,26 +0,0 @@
-<!-- 
-Thanks for contributing a merge request! Please name and fully describe your MR as you would for a commit message.
-If the MR fixes an issue, please include "Fixes #issue" in the commit message and the MR description.
-
-In addition, we recommend that first-time contributors read our [contribution guidelines](https://eigen.tuxfamily.org/index.php?title=Contributing_to_Eigen) and [git page](https://eigen.tuxfamily.org/index.php?title=Git), which will help you submit a more standardized MR.
-
-Before submitting the MR, you also need to complete the following checks:
- Make one PR per feature/bugfix (don't mix multiple changes into one PR). Avoid committing unrelated changes.
- Rebase before committing
- For code changes, run the test suite (at least the tests that are likely affected by the change).
-  See our [test guidelines](https://eigen.tuxfamily.org/index.php?title=Tests).
- If possible, add a test (both for bug-fixes as well as new features)
- Make sure new features are documented
-
-Note that we are a team of volunteers; we appreciate your patience during the review process.
-
-Again, thanks for contributing! -->
-
-### Reference issue
-<!-- You can link to a specific issue using the gitlab syntax #<issue number>  -->
-
-### What does this implement/fix?
-<!--Please explain your changes.-->
-
-### Additional information
-<!--Any additional information you think is important.-->
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,15 +1,6 @@
 # Changelog

-## [5.0.1] - 2025-11-11
-
-A few bug-fixes from the master branch, including
- Dirty git state [#2995]
- Failing geo_homogeneous tests [#2977]
- Alignment issues [#2982, #2984]
- Missing C++20 `<version>` header [#2986]
- BLAS/LAPACK build on windows [#2980]
-
-See the full lists of [addressed bugs](https://gitlab.com/libeigen/eigen/-/issues?state=all&label_name%5B%5D=release%3A%3A5.0.1) and [merge requests](https://gitlab.com/libeigen/eigen/-/merge_requests?state=all&label_name%5B%5D=release%3A%3A5.0.1) for more details.
+## [Unreleased]

 ## [5.0.0] - 2025-09-30

@@ -45,6 +36,14 @@ This release marks a transition to [Semantic Versioning](https://semver.org/). P
 * Euler angles are now returned in a more canonical form, potentially resulting in a change of behavior [!1301, !1314].
 * Eigen's random number generation has changed, resulting in a change of behavior.  Please do not rely on specific random numbers from Eigen - these were never guaranteed to be consistent across Eigen versions, nor are they generally consistent across platforms [!1437].

+## [3.4.1] - 2025-09-30
+
+Many bug fixes have been backported from the main branch.
+
+A list of new issues addressed can be found via the [3.4.1](https://gitlab.com/libeigen/eigen/-/issues?state=all&label_name%5B%5D=3.4.1) label on GitLab.
+
+Check the [git commit history](https://gitlab.com/libeigen/eigen/-/commits/3.4.1) for the full list of changes.  
+
 ## [3.4.0] - 2021-08-18

 **Notice:** 3.4.x will be the last major release series of Eigen that will support c++03.
@@ -344,7 +343,7 @@ Changes since 3.3.7:
  * Commit dd6de618: Fix a bug with half-precision floats on GPUs.

 ## [3.3.7] - 2018-12-11
-¸¸¸¸
+
 Changes since 3.3.6:

 * #1643: Fix compilation with GCC>=6 and compiler optimization turned off.
@@ -528,6 +527,7 @@ Changes since 3.3.2:
  * #1378: fix doc (`DiagonalIndex` vs `Diagonal`).

 ## [3.3.2] - 2017-01-18
+
 Changes since 3.3.1:

 * General:
@@ -845,6 +845,7 @@ Changes since 3.2.8:
  * #1175: fix index type conversion warnings in sparse to dense conversion.

 ## [3.2.8] - 2016-02-16
+
 Changes since 3.2.7:

 * Main fixes and improvements:
@@ -1268,6 +1269,7 @@ Main changes since 3.2-beta1:
  * Many other fixes including #230, #482, #542, #561, #564, #565, #566, #578, #581, #595, #597, #598, #599, #605, #606, #615.

 ## [3.1.3] - 2013-04-16
+
 Changes since 3.1.2:

 * #526 - Fix linear vectorized transversal in linspace.
@@ -1351,6 +1353,7 @@ Changes since 3.1.0:
 * Fixed Sparse module compilation under MSVC 2005

 ## [3.0.6] - 2012-07-09
+
 Changes since 3.0.5:
 * #447 - fix infinite recursion in `ProductBase::coeff()`
 * #478 - fix RealSchur on a zero matrix
@@ -1507,6 +1510,7 @@ Main changes since 3.0:


 ## [3.0.4] - 2011-12-06
+
 Changes since 3.0.3:

 * #363 - check for integer overflow in size computations
--- a/Eigen/Eigenvalues
+++ b/Eigen/Eigenvalues
@@ -15,6 +15,7 @@
 #include "Householder"
 #include "LU"
 #include "Geometry"
+#include "Sparse"  // Needed by ComplexQZ.

 #include "src/Core/util/DisableStupidWarnings.h"

@@ -32,8 +33,6 @@
 * \endcode
 */

-#include "src/misc/RealSvd2x2.h"
-
 // IWYU pragma: begin_exports
 #include "src/Eigenvalues/Tridiagonalization.h"
 #include "src/Eigenvalues/RealSchur.h"
@@ -44,6 +43,7 @@
 #include "src/Eigenvalues/ComplexSchur.h"
 #include "src/Eigenvalues/ComplexEigenSolver.h"
 #include "src/Eigenvalues/RealQZ.h"
+#include "src/Eigenvalues/ComplexQZ.h"
 #include "src/Eigenvalues/GeneralizedEigenSolver.h"
 #include "src/Eigenvalues/MatrixBaseEigenvalues.h"
 #ifdef EIGEN_USE_LAPACKE
--- a/Eigen/SVD
+++ b/Eigen/SVD
@@ -33,7 +33,6 @@
 */

 // IWYU pragma: begin_exports
-#include "src/misc/RealSvd2x2.h"
 #include "src/SVD/UpperBidiagonalization.h"
 #include "src/SVD/SVDBase.h"
 #include "src/SVD/JacobiSVD.h"
--- a/Eigen/Version
+++ b/Eigen/Version
@@ -7,8 +7,8 @@
 #define EIGEN_MAJOR_VERSION 5
 #define EIGEN_MINOR_VERSION 0
 #define EIGEN_PATCH_VERSION 1
-#define EIGEN_PRERELEASE_VERSION ""
-#define EIGEN_BUILD_VERSION ""
-#define EIGEN_VERSION_STRING "5.0.1"
+#define EIGEN_PRERELEASE_VERSION "dev"
+#define EIGEN_BUILD_VERSION "master"
+#define EIGEN_VERSION_STRING "5.0.1-dev+master"

 #endif  // EIGEN_VERSION_H
--- a/Eigen/src/Core/DenseCoeffsBase.h
+++ b/Eigen/src/Core/DenseCoeffsBase.h
@@ -113,6 +113,16 @@ class DenseCoeffsBase<Derived, ReadOnlyAccessors> : public EigenBase<Derived> {
    return coeff(row, col);
  }

+#ifdef EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+  /** \returns the coefficient at given the given row and column.
+   *
+   * \sa operator[](Index,Index), operator[](Index)
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr CoeffReturnType operator[](Index row, Index col) const {
+    return operator()(row, col);
+  }
+#endif
+
  /** Short version: don't use this function, use
   * \link operator[](Index) const \endlink instead.
   *
@@ -316,12 +326,21 @@ class DenseCoeffsBase<Derived, WriteAccessors> : public DenseCoeffsBase<Derived,
   *
   * \sa operator[](Index)
   */
-
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Scalar& operator()(Index row, Index col) {
    eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols());
    return coeffRef(row, col);
  }

+#ifdef EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+  /** \returns a reference to the coefficient at given the given row and column.
+   *
+   * \sa operator[](Index)
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Scalar& operator[](Index row, Index col) {
+    return operator()(row, col);
+  }
+#endif
+
  /** Short version: don't use this function, use
   * \link operator[](Index) \endlink instead.
   *
--- a/Eigen/src/Core/GenericPacketMath.h
+++ b/Eigen/src/Core/GenericPacketMath.h
@@ -63,7 +63,6 @@ struct default_packet_traits {

    HasArg = 0,
    HasAbsDiff = 0,
-    HasBlend = 0,
    // This flag is used to indicate whether packet comparison is supported.
    // pcmp_eq and pcmp_lt should be defined for it to be true.
    HasCmp = 0,
@@ -1482,21 +1481,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet, 1>& /*kernel*/) {
  // Nothing to do in the scalar case, i.e. a 1x1 matrix.
 }

-/***************************************************************************
- * Selector, i.e. vector of N boolean values used to select (i.e. blend)
- * words from 2 packets.
- ***************************************************************************/
-template <size_t N>
-struct Selector {
-  bool select[N];
-};
-
-template <typename Packet>
-EIGEN_DEVICE_FUNC inline Packet pblend(const Selector<unpacket_traits<Packet>::size>& ifPacket,
-                                       const Packet& thenPacket, const Packet& elsePacket) {
-  return ifPacket.select[0] ? thenPacket : elsePacket;
-}
-
 /** \internal \returns 1 / a (coeff-wise) */
 template <typename Packet>
 EIGEN_DEVICE_FUNC inline Packet preciprocal(const Packet& a) {
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -89,6 +89,11 @@ class TriangularBase : public EigenBase<Derived> {
    return coeffRef(row, col);
  }

+#ifdef EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+  EIGEN_DEVICE_FUNC inline Scalar operator[](Index row, Index col) const { return operator()(row, col); }
+  EIGEN_DEVICE_FUNC inline Scalar& operator[](Index row, Index col) { return operator()(row, col); }
+#endif
+
 #ifndef EIGEN_PARSED_BY_DOXYGEN
  EIGEN_DEVICE_FUNC inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
  EIGEN_DEVICE_FUNC inline Derived& derived() { return *static_cast<Derived*>(this); }
--- a/Eigen/src/Core/arch/AVX/MathFunctions.h
+++ b/Eigen/src/Core/arch/AVX/MathFunctions.h
@@ -25,8 +25,8 @@ EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_FLOAT(Packet8f)

 EIGEN_DOUBLE_PACKET_FUNCTION(atanh, Packet4d)
 EIGEN_DOUBLE_PACKET_FUNCTION(log, Packet4d)
-EIGEN_DOUBLE_PACKET_FUNCTION(log2, Packet4d)
 EIGEN_DOUBLE_PACKET_FUNCTION(exp, Packet4d)
+EIGEN_DOUBLE_PACKET_FUNCTION(log2, Packet4d)
 EIGEN_DOUBLE_PACKET_FUNCTION(tanh, Packet4d)
 EIGEN_DOUBLE_PACKET_FUNCTION(cbrt, Packet4d)
 #ifdef EIGEN_VECTORIZE_AVX2
@@ -35,6 +35,8 @@ EIGEN_DOUBLE_PACKET_FUNCTION(cos, Packet4d)
 #endif
 EIGEN_GENERIC_PACKET_FUNCTION(atan, Packet4d)
 EIGEN_GENERIC_PACKET_FUNCTION(exp2, Packet4d)
+EIGEN_GENERIC_PACKET_FUNCTION(expm1, Packet4d)
+EIGEN_GENERIC_PACKET_FUNCTION(log1p, Packet4d)

 // Notice that for newer processors, it is counterproductive to use Newton
 // iteration for square root. In particular, Skylake and Zen2 processors
--- a/Eigen/src/Core/arch/AVX/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX/PacketMath.h
@@ -115,9 +115,9 @@ struct packet_traits<float> : default_packet_traits {
    HasATan = 1,
    HasATanh = 1,
    HasLog = 1,
+    HasExp = 1,
    HasLog1p = 1,
    HasExpm1 = 1,
-    HasExp = 1,
    HasPow = 1,
    HasNdtri = 1,
    HasBessel = 1,
@@ -127,7 +127,6 @@ struct packet_traits<float> : default_packet_traits {
    HasTanh = EIGEN_FAST_MATH,
    HasErf = EIGEN_FAST_MATH,
    HasErfc = EIGEN_FAST_MATH,
-    HasBlend = 1
  };
 };
 template <>
@@ -146,17 +145,18 @@ struct packet_traits<double> : default_packet_traits {
    HasCos = EIGEN_FAST_MATH,
 #endif
    HasTanh = EIGEN_FAST_MATH,
-    HasLog = 1,
    HasErf = 1,
    HasErfc = 1,
+    HasLog = 1,
    HasExp = 1,
+    HasLog1p = 1,
+    HasExpm1 = 1,
    HasPow = 1,
    HasSqrt = 1,
    HasRsqrt = 1,
    HasCbrt = 1,
    HasATan = 1,
    HasATanh = 1,
-    HasBlend = 1
  };
 };

@@ -179,7 +179,6 @@ struct packet_traits<Eigen::half> : default_packet_traits {
    HasCos = EIGEN_FAST_MATH,
    HasNegate = 1,
    HasAbs = 1,
-    HasAbs2 = 0,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
@@ -192,7 +191,6 @@ struct packet_traits<Eigen::half> : default_packet_traits {
    HasRsqrt = 1,
    HasTanh = EIGEN_FAST_MATH,
    HasErf = EIGEN_FAST_MATH,
-    HasBlend = 0,
    HasBessel = 1,
    HasNdtri = 1
  };
@@ -218,7 +216,6 @@ struct packet_traits<bfloat16> : default_packet_traits {
    HasCos = EIGEN_FAST_MATH,
    HasNegate = 1,
    HasAbs = 1,
-    HasAbs2 = 0,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
@@ -231,7 +228,6 @@ struct packet_traits<bfloat16> : default_packet_traits {
    HasRsqrt = 1,
    HasTanh = EIGEN_FAST_MATH,
    HasErf = EIGEN_FAST_MATH,
-    HasBlend = 0,
    HasBessel = 1,
    HasNdtri = 1
  };
@@ -284,7 +280,6 @@ struct packet_traits<uint64_t> : default_packet_traits {
    // HasMin = 0,
    // HasMax = 0,
    HasDiv = 0,
-    HasBlend = 0,
    HasTranspose = 0,
    HasNegate = 0,
    HasSqrt = 0,
@@ -2070,31 +2065,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet4d, 4>& kernel) {
  kernel.packet[2] = _mm256_permute2f128_pd(T1, T3, 49);
 }

-EIGEN_STRONG_INLINE __m256i avx_blend_mask(const Selector<4>& ifPacket) {
-  return _mm256_set_epi64x(0 - ifPacket.select[3], 0 - ifPacket.select[2], 0 - ifPacket.select[1],
-                           0 - ifPacket.select[0]);
-}
-
-EIGEN_STRONG_INLINE __m256i avx_blend_mask(const Selector<8>& ifPacket) {
-  return _mm256_set_epi32(0 - ifPacket.select[7], 0 - ifPacket.select[6], 0 - ifPacket.select[5],
-                          0 - ifPacket.select[4], 0 - ifPacket.select[3], 0 - ifPacket.select[2],
-                          0 - ifPacket.select[1], 0 - ifPacket.select[0]);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet8f pblend(const Selector<8>& ifPacket, const Packet8f& thenPacket,
-                                    const Packet8f& elsePacket) {
-  const __m256 true_mask = _mm256_castsi256_ps(avx_blend_mask(ifPacket));
-  return pselect<Packet8f>(true_mask, thenPacket, elsePacket);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet4d pblend(const Selector<4>& ifPacket, const Packet4d& thenPacket,
-                                    const Packet4d& elsePacket) {
-  const __m256d true_mask = _mm256_castsi256_pd(avx_blend_mask(ifPacket));
-  return pselect<Packet4d>(true_mask, thenPacket, elsePacket);
-}
-
 // Packet math for Eigen::half
 #ifndef EIGEN_VECTORIZE_AVX512FP16
 template <>
--- a/Eigen/src/Core/arch/AVX512/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX512/PacketMath.h
@@ -84,7 +84,6 @@ struct packet_traits<half> : default_packet_traits {
    HasDiv = 1,
    HasNegate = 1,
    HasAbs = 1,
-    HasAbs2 = 0,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
@@ -160,6 +159,8 @@ struct packet_traits<double> : default_packet_traits {
    HasCos = EIGEN_FAST_MATH,
    HasLog = 1,
    HasExp = 1,
+    HasLog1p = 1,
+    HasExpm1 = 1,
    HasPow = 1,
    HasATan = 1,
    HasTanh = EIGEN_FAST_MATH,
@@ -2057,27 +2058,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet16i, 4>& kernel) {
  PACK_OUTPUT_I32_2(kernel.packet, tmp.packet, 3, 1);
 }

-template <size_t N>
-EIGEN_STRONG_INLINE int avx512_blend_mask(const Selector<N>& ifPacket) {
-  alignas(__m128i) uint8_t aux[sizeof(__m128i)];
-  for (size_t i = 0; i < N; i++) aux[i] = static_cast<uint8_t>(ifPacket.select[i]);
-  __m128i paux = _mm_sub_epi8(_mm_setzero_si128(), _mm_load_si128(reinterpret_cast<const __m128i*>(aux)));
-  return _mm_movemask_epi8(paux);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet16f pblend(const Selector<16>& ifPacket, const Packet16f& thenPacket,
-                                     const Packet16f& elsePacket) {
-  __mmask16 m = avx512_blend_mask(ifPacket);
-  return _mm512_mask_blend_ps(m, elsePacket, thenPacket);
-}
-template <>
-EIGEN_STRONG_INLINE Packet8d pblend(const Selector<8>& ifPacket, const Packet8d& thenPacket,
-                                    const Packet8d& elsePacket) {
-  __mmask8 m = avx512_blend_mask(ifPacket);
-  return _mm512_mask_blend_pd(m, elsePacket, thenPacket);
-}
-
 // Packet math for Eigen::half
 #ifndef EIGEN_VECTORIZE_AVX512FP16
 template <>
--- a/Eigen/src/Core/arch/AVX512/PacketMathFP16.h
+++ b/Eigen/src/Core/arch/AVX512/PacketMathFP16.h
@@ -42,7 +42,6 @@ struct packet_traits<half> : default_packet_traits {
    HasDiv = 1,
    HasNegate = 1,
    HasAbs = 1,
-    HasAbs2 = 0,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
--- a/Eigen/src/Core/arch/AltiVec/Complex.h
+++ b/Eigen/src/Core/arch/AltiVec/Complex.h
@@ -109,9 +109,6 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
    HasSqrt = 1,
    HasLog = 1,
    HasExp = 1,
-#ifdef EIGEN_VECTORIZE_VSX
-    HasBlend = 1,
-#endif
    HasSetLinear = 0
  };
 };
@@ -364,17 +361,6 @@ EIGEN_STRONG_INLINE Packet2cf pcmp_eq(const Packet2cf& a, const Packet2cf& b) {
  return Packet2cf(vec_and(eq, vec_perm(eq, eq, p16uc_COMPLEX32_REV)));
 }

-#ifdef EIGEN_VECTORIZE_VSX
-template <>
-EIGEN_STRONG_INLINE Packet2cf pblend(const Selector<2>& ifPacket, const Packet2cf& thenPacket,
-                                     const Packet2cf& elsePacket) {
-  Packet2cf result;
-  result.v = reinterpret_cast<Packet4f>(
-      pblend<Packet2d>(ifPacket, reinterpret_cast<Packet2d>(thenPacket.v), reinterpret_cast<Packet2d>(elsePacket.v)));
-  return result;
-}
-#endif
-
 template <>
 EIGEN_STRONG_INLINE Packet2cf psqrt<Packet2cf>(const Packet2cf& a) {
  return psqrt_complex<Packet2cf>(a);
--- a/Eigen/src/Core/arch/AltiVec/PacketMath.h
+++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h
@@ -184,6 +184,8 @@ struct packet_traits<float> : default_packet_traits {
    HasATanh = 1,
    HasLog = 1,
    HasExp = 1,
+    HasLog1p = 1,
+    HasExpm1 = 1,
 #ifdef EIGEN_VECTORIZE_VSX
    HasCmp = 1,
    HasPow = 1,
@@ -204,7 +206,6 @@ struct packet_traits<float> : default_packet_traits {
    HasErf = 0,
 #endif
    HasNegate = 1,
-    HasBlend = 1
  };
 };
 template <>
@@ -241,7 +242,6 @@ struct packet_traits<bfloat16> : default_packet_traits {
    HasTanh = 0,
    HasErf = 0,
    HasNegate = 1,
-    HasBlend = 1
  };
 };

@@ -263,7 +263,6 @@ struct packet_traits<int> : default_packet_traits {
 #else
    HasDiv = 0,
 #endif
-    HasBlend = 1,
    HasCmp = 1
  };
 };
@@ -281,7 +280,6 @@ struct packet_traits<short int> : default_packet_traits {
    HasSub = 1,
    HasMul = 1,
    HasDiv = 0,
-    HasBlend = 1,
    HasCmp = 1
  };
 };
@@ -299,7 +297,6 @@ struct packet_traits<unsigned short int> : default_packet_traits {
    HasSub = 1,
    HasMul = 1,
    HasDiv = 0,
-    HasBlend = 1,
    HasCmp = 1
  };
 };
@@ -317,7 +314,6 @@ struct packet_traits<signed char> : default_packet_traits {
    HasSub = 1,
    HasMul = 1,
    HasDiv = 0,
-    HasBlend = 1,
    HasCmp = 1
  };
 };
@@ -335,7 +331,6 @@ struct packet_traits<unsigned char> : default_packet_traits {
    HasSub = 1,
    HasMul = 1,
    HasDiv = 0,
-    HasBlend = 1,
    HasCmp = 1
  };
 };
@@ -3053,74 +3048,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet16uc, 16>& kernel) {
  kernel.packet[15] = vec_mergel(step3[7], step3[15]);
 }

-template <typename Packet>
-EIGEN_STRONG_INLINE Packet pblend4(const Selector<4>& ifPacket, const Packet& thenPacket, const Packet& elsePacket) {
-  Packet4ui select = {ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3]};
-  Packet4ui mask = reinterpret_cast<Packet4ui>(pnegate(reinterpret_cast<Packet4i>(select)));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet4i pblend(const Selector<4>& ifPacket, const Packet4i& thenPacket,
-                                    const Packet4i& elsePacket) {
-  return pblend4<Packet4i>(ifPacket, thenPacket, elsePacket);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet4f pblend(const Selector<4>& ifPacket, const Packet4f& thenPacket,
-                                    const Packet4f& elsePacket) {
-  return pblend4<Packet4f>(ifPacket, thenPacket, elsePacket);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet8s pblend(const Selector<8>& ifPacket, const Packet8s& thenPacket,
-                                    const Packet8s& elsePacket) {
-  Packet8us select = {ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3],
-                      ifPacket.select[4], ifPacket.select[5], ifPacket.select[6], ifPacket.select[7]};
-  Packet8us mask = reinterpret_cast<Packet8us>(pnegate(reinterpret_cast<Packet8s>(select)));
-  Packet8s result = vec_sel(elsePacket, thenPacket, mask);
-  return result;
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet8us pblend(const Selector<8>& ifPacket, const Packet8us& thenPacket,
-                                     const Packet8us& elsePacket) {
-  Packet8us select = {ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3],
-                      ifPacket.select[4], ifPacket.select[5], ifPacket.select[6], ifPacket.select[7]};
-  Packet8us mask = reinterpret_cast<Packet8us>(pnegate(reinterpret_cast<Packet8s>(select)));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet8bf pblend(const Selector<8>& ifPacket, const Packet8bf& thenPacket,
-                                     const Packet8bf& elsePacket) {
-  return pblend<Packet8us>(ifPacket, thenPacket, elsePacket);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet16c pblend(const Selector<16>& ifPacket, const Packet16c& thenPacket,
-                                     const Packet16c& elsePacket) {
-  Packet16uc select = {ifPacket.select[0],  ifPacket.select[1],  ifPacket.select[2],  ifPacket.select[3],
-                       ifPacket.select[4],  ifPacket.select[5],  ifPacket.select[6],  ifPacket.select[7],
-                       ifPacket.select[8],  ifPacket.select[9],  ifPacket.select[10], ifPacket.select[11],
-                       ifPacket.select[12], ifPacket.select[13], ifPacket.select[14], ifPacket.select[15]};
-
-  Packet16uc mask = reinterpret_cast<Packet16uc>(pnegate(reinterpret_cast<Packet16c>(select)));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet16uc pblend(const Selector<16>& ifPacket, const Packet16uc& thenPacket,
-                                      const Packet16uc& elsePacket) {
-  Packet16uc select = {ifPacket.select[0],  ifPacket.select[1],  ifPacket.select[2],  ifPacket.select[3],
-                       ifPacket.select[4],  ifPacket.select[5],  ifPacket.select[6],  ifPacket.select[7],
-                       ifPacket.select[8],  ifPacket.select[9],  ifPacket.select[10], ifPacket.select[11],
-                       ifPacket.select[12], ifPacket.select[13], ifPacket.select[14], ifPacket.select[15]};
-
-  Packet16uc mask = reinterpret_cast<Packet16uc>(pnegate(reinterpret_cast<Packet16c>(select)));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
 //---------- double ----------
 #ifdef EIGEN_VECTORIZE_VSX
 typedef __vector double Packet2d;
@@ -3176,9 +3103,11 @@ struct packet_traits<double> : default_packet_traits {
    HasErfc = EIGEN_FAST_MATH,
    HasATanh = 1,
    HasATan = 0,
-    HasLog = 0,
    HasCmp = 1,
+    HasLog = 1,
    HasExp = 1,
+    HasLog1p = 1,
+    HasExpm1 = 1,
    HasSqrt = 1,
    HasCbrt = 1,
 #if !EIGEN_COMP_CLANG
@@ -3187,7 +3116,6 @@ struct packet_traits<double> : default_packet_traits {
    HasRsqrt = 0,
 #endif
    HasNegate = 1,
-    HasBlend = 1
  };
 };

@@ -3713,14 +3641,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet2d, 2>& kernel) {
  kernel.packet[1] = t1;
 }

-template <>
-EIGEN_STRONG_INLINE Packet2d pblend(const Selector<2>& ifPacket, const Packet2d& thenPacket,
-                                    const Packet2d& elsePacket) {
-  Packet2l select = {ifPacket.select[0], ifPacket.select[1]};
-  Packet2ul mask = reinterpret_cast<Packet2ul>(pnegate(reinterpret_cast<Packet2l>(select)));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
 #endif  // __VSX__
 }  // end namespace internal

--- a/Eigen/src/Core/arch/Default/ConjHelper.h
+++ b/Eigen/src/Core/arch/Default/ConjHelper.h
@@ -17,6 +17,9 @@
    EIGEN_STRONG_INLINE PACKET_CPLX pmadd(const PACKET_REAL& x, const PACKET_CPLX& y, const PACKET_CPLX& c) const { \
      return padd(c, this->pmul(x, y));                                                                             \
    }                                                                                                               \
+    EIGEN_STRONG_INLINE PACKET_CPLX pmsub(const PACKET_REAL& x, const PACKET_CPLX& y, const PACKET_CPLX& c) const { \
+      return psub(this->pmul(x, y), c);                                                                             \
+    }                                                                                                               \
    EIGEN_STRONG_INLINE PACKET_CPLX pmul(const PACKET_REAL& x, const PACKET_CPLX& y) const {                        \
      return PACKET_CPLX(Eigen::internal::pmul<PACKET_REAL>(x, y.v));                                               \
    }                                                                                                               \
@@ -27,6 +30,9 @@
    EIGEN_STRONG_INLINE PACKET_CPLX pmadd(const PACKET_CPLX& x, const PACKET_REAL& y, const PACKET_CPLX& c) const { \
      return padd(c, this->pmul(x, y));                                                                             \
    }                                                                                                               \
+    EIGEN_STRONG_INLINE PACKET_CPLX pmsub(const PACKET_CPLX& x, const PACKET_REAL& y, const PACKET_CPLX& c) const { \
+      return psub(this->pmul(x, y), c);                                                                             \
+    }                                                                                                               \
    EIGEN_STRONG_INLINE PACKET_CPLX pmul(const PACKET_CPLX& x, const PACKET_REAL& y) const {                        \
      return PACKET_CPLX(Eigen::internal::pmul<PACKET_REAL>(x.v, y));                                               \
    }                                                                                                               \
@@ -76,6 +82,11 @@ struct conj_helper {
    return this->pmul(x, y) + c;
  }

+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ResultType pmsub(const LhsType& x, const RhsType& y,
+                                                         const ResultType& c) const {
+    return this->pmul(x, y) - c;
+  }
+
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ResultType pmul(const LhsType& x, const RhsType& y) const {
    return conj_if<ConjLhs>()(x) * conj_if<ConjRhs>()(y);
  }
@@ -104,6 +115,10 @@ struct conj_helper<Packet, Packet, ConjLhs, ConjRhs> {
    return Eigen::internal::pmadd(conj_if<ConjLhs>().pconj(x), conj_if<ConjRhs>().pconj(y), c);
  }

+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmsub(const Packet& x, const Packet& y, const Packet& c) const {
+    return Eigen::internal::pmsub(conj_if<ConjLhs>().pconj(x), conj_if<ConjRhs>().pconj(y), c);
+  }
+
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmul(const Packet& x, const Packet& y) const {
    return Eigen::internal::pmul(conj_if<ConjLhs>().pconj(x), conj_if<ConjRhs>().pconj(y));
  }
@@ -116,6 +131,9 @@ struct conj_helper<Packet, Packet, true, true> {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmadd(const Packet& x, const Packet& y, const Packet& c) const {
    return Eigen::internal::pmadd(pconj(x), pconj(y), c);
  }
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmsub(const Packet& x, const Packet& y, const Packet& c) const {
+    return Eigen::internal::pmsub(pconj(x), pconj(y), c);
+  }
  // We save a conjuation by using the identity conj(a)*conj(b) = conj(a*b).
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmul(const Packet& x, const Packet& y) const {
    return pconj(Eigen::internal::pmul(x, y));
--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctionsFwd.h
@@ -210,16 +210,18 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_round(const Packet& a);
  EIGEN_GENERIC_PACKET_FUNCTION(atan, PACKET)

 #define EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_DOUBLE(PACKET) \
-  EIGEN_DOUBLE_PACKET_FUNCTION(atanh, PACKET)               \
-  EIGEN_DOUBLE_PACKET_FUNCTION(log, PACKET)                 \
  EIGEN_DOUBLE_PACKET_FUNCTION(sin, PACKET)                 \
  EIGEN_DOUBLE_PACKET_FUNCTION(cos, PACKET)                 \
+  EIGEN_DOUBLE_PACKET_FUNCTION(log, PACKET)                 \
  EIGEN_DOUBLE_PACKET_FUNCTION(log2, PACKET)                \
  EIGEN_DOUBLE_PACKET_FUNCTION(exp, PACKET)                 \
  EIGEN_DOUBLE_PACKET_FUNCTION(tanh, PACKET)                \
+  EIGEN_DOUBLE_PACKET_FUNCTION(atanh, PACKET)               \
  EIGEN_DOUBLE_PACKET_FUNCTION(cbrt, PACKET)                \
-  EIGEN_GENERIC_PACKET_FUNCTION(atan, PACKET)               \
-  EIGEN_GENERIC_PACKET_FUNCTION(exp2, PACKET)
+  EIGEN_GENERIC_PACKET_FUNCTION(expm1, PACKET)              \
+  EIGEN_GENERIC_PACKET_FUNCTION(exp2, PACKET)               \
+  EIGEN_GENERIC_PACKET_FUNCTION(log1p, PACKET)              \
+  EIGEN_GENERIC_PACKET_FUNCTION(atan, PACKET)

 }  // end namespace internal
 }  // end namespace Eigen
--- a/Eigen/src/Core/arch/HVX/PacketMath.h
+++ b/Eigen/src/Core/arch/HVX/PacketMath.h
@@ -152,7 +152,6 @@ struct packet_traits<float> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 0,
    HasAbsDiff = 0,
    HasMin = 1,
    HasMax = 1,
--- a/Eigen/src/Core/arch/LSX/PacketMath.h
+++ b/Eigen/src/Core/arch/LSX/PacketMath.h
@@ -169,10 +169,8 @@ struct packet_traits<int8_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 16,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasCmp = 1,
-    HasBlend = 0
  };
 };

@@ -185,11 +183,9 @@ struct packet_traits<int16_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 8,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasCmp = 1,
    HasDiv = 1,
-    HasBlend = 0
  };
 };

@@ -202,11 +198,9 @@ struct packet_traits<int32_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 4,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasCmp = 1,
    HasDiv = 1,
-    HasBlend = 0
  };
 };

@@ -219,11 +213,9 @@ struct packet_traits<int64_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 2,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasCmp = 1,
    HasDiv = 1,
-    HasBlend = 0
  };
 };

@@ -236,11 +228,9 @@ struct packet_traits<uint8_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 16,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasNegate = 0,
    HasCmp = 1,
-    HasBlend = 0
  };
 };

@@ -253,12 +243,10 @@ struct packet_traits<uint16_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 8,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasNegate = 0,
    HasCmp = 1,
    HasDiv = 1,
-    HasBlend = 0
  };
 };

@@ -271,12 +259,10 @@ struct packet_traits<uint32_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 4,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasNegate = 0,
    HasCmp = 1,
    HasDiv = 1,
-    HasBlend = 0
  };
 };

@@ -289,12 +275,10 @@ struct packet_traits<uint64_t> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 2,

-    HasAbs2 = 0,
    HasSetLinear = 0,
    HasNegate = 0,
    HasCmp = 1,
    HasDiv = 1,
-    HasBlend = 0
  };
 };

@@ -307,9 +291,7 @@ struct packet_traits<float> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 4,

-    HasAbs2 = 0,
    HasSetLinear = 0,
-    HasBlend = 0,
    HasSign = 0,
    HasDiv = 1,
    HasExp = 1,
@@ -328,9 +310,7 @@ struct packet_traits<double> : default_packet_traits {
    AlignedOnScalar = 1,
    size = 2,

-    HasAbs2 = 0,
    HasSetLinear = 0,
-    HasBlend = 0,
    HasSign = 0,
    HasDiv = 1,
    HasSqrt = 1,
--- a/Eigen/src/Core/arch/MSA/Complex.h
+++ b/Eigen/src/Core/arch/MSA/Complex.h
@@ -105,7 +105,6 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
    HasMin = 0,
    HasMax = 0,
    HasSetLinear = 0,
-    HasBlend = 1
  };
 };

@@ -314,12 +313,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet2cf, 2>& kernel) {
  kernel.packet[1].v = tmp;
 }

-template <>
-EIGEN_STRONG_INLINE Packet2cf pblend(const Selector<2>& ifPacket, const Packet2cf& thenPacket,
-                                     const Packet2cf& elsePacket) {
-  return (Packet2cf)(Packet4f)pblend<Packet2d>(ifPacket, (Packet2d)thenPacket.v, (Packet2d)elsePacket.v);
-}
-
 //---------- double ----------

 struct Packet1cd {
--- a/Eigen/src/Core/arch/MSA/PacketMath.h
+++ b/Eigen/src/Core/arch/MSA/PacketMath.h
@@ -91,7 +91,6 @@ struct packet_traits<float> : default_packet_traits {
    HasExp = 1,
    HasSqrt = 1,
    HasRsqrt = 1,
-    HasBlend = 1
  };
 };

@@ -105,7 +104,6 @@ struct packet_traits<int32_t> : default_packet_traits {
    size = 4,
    // FIXME check the Has*
    HasDiv = 1,
-    HasBlend = 1
  };
 };

@@ -802,22 +800,6 @@ EIGEN_STRONG_INLINE Packet4f pround<Packet4f>(const Packet4f& a) {
  return v;
 }

-template <>
-EIGEN_STRONG_INLINE Packet4f pblend(const Selector<4>& ifPacket, const Packet4f& thenPacket,
-                                    const Packet4f& elsePacket) {
-  Packet4ui select = {ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3]};
-  Packet4i mask = __builtin_msa_ceqi_w((Packet4i)select, 0);
-  return (Packet4f)__builtin_msa_bsel_v((v16u8)mask, (v16u8)thenPacket, (v16u8)elsePacket);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet4i pblend(const Selector<4>& ifPacket, const Packet4i& thenPacket,
-                                    const Packet4i& elsePacket) {
-  Packet4ui select = {ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3]};
-  Packet4i mask = __builtin_msa_ceqi_w((Packet4i)select, 0);
-  return (Packet4i)__builtin_msa_bsel_v((v16u8)mask, (v16u8)thenPacket, (v16u8)elsePacket);
-}
-
 //---------- double ----------

 typedef v2f64 Packet2d;
@@ -856,7 +838,6 @@ struct packet_traits<double> : default_packet_traits {
    HasExp = 1,
    HasSqrt = 1,
    HasRsqrt = 1,
-    HasBlend = 1
  };
 };

@@ -1222,14 +1203,6 @@ EIGEN_STRONG_INLINE Packet2d pround<Packet2d>(const Packet2d& a) {
  return v;
 }

-template <>
-EIGEN_STRONG_INLINE Packet2d pblend(const Selector<2>& ifPacket, const Packet2d& thenPacket,
-                                    const Packet2d& elsePacket) {
-  Packet2ul select = {ifPacket.select[0], ifPacket.select[1]};
-  Packet2l mask = __builtin_msa_ceqi_d((Packet2l)select, 0);
-  return (Packet2d)__builtin_msa_bsel_v((v16u8)mask, (v16u8)thenPacket, (v16u8)elsePacket);
-}
-
 }  // end namespace internal

 }  // end namespace Eigen
--- a/Eigen/src/Core/arch/NEON/PacketMath.h
+++ b/Eigen/src/Core/arch/NEON/PacketMath.h
@@ -189,13 +189,11 @@ struct packet_traits<float> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0,
    HasDiv = 1,
    HasSin = EIGEN_FAST_MATH,
    HasCos = EIGEN_FAST_MATH,
@@ -205,6 +203,8 @@ struct packet_traits<float> : default_packet_traits {
    HasATanh = 1,
    HasLog = 1,
    HasExp = 1,
+    HasLog1p = 1,
+    HasExpm1 = 1,
    HasPow = 1,
    HasSqrt = 1,
    HasRsqrt = 1,
@@ -235,12 +235,10 @@ struct packet_traits<int8_t> : default_packet_traits {
    HasAbs = 1,
    HasAbsDiff = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0
  };
 };

@@ -262,12 +260,10 @@ struct packet_traits<uint8_t> : default_packet_traits {
    HasAbs = 1,
    HasAbsDiff = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0,

    HasSqrt = 1
  };
@@ -291,12 +287,10 @@ struct packet_traits<int16_t> : default_packet_traits {
    HasAbs = 1,
    HasAbsDiff = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0
  };
 };

@@ -318,12 +312,10 @@ struct packet_traits<uint16_t> : default_packet_traits {
    HasAbs = 1,
    HasAbsDiff = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0,
    HasSqrt = 1
  };
 };
@@ -345,13 +337,11 @@ struct packet_traits<int32_t> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0
  };
 };

@@ -372,13 +362,11 @@ struct packet_traits<uint32_t> : default_packet_traits {
    HasNegate = 0,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0,

    HasSqrt = 1
  };
@@ -401,13 +389,11 @@ struct packet_traits<int64_t> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0
  };
 };

@@ -428,13 +414,11 @@ struct packet_traits<uint64_t> : default_packet_traits {
    HasNegate = 0,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0
  };
 };

@@ -4631,13 +4615,11 @@ struct packet_traits<bfloat16> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0,
    HasDiv = 1,
    HasSin = EIGEN_FAST_MATH,
    HasCos = EIGEN_FAST_MATH,
@@ -5016,19 +4998,19 @@ struct packet_traits<double> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0,

    HasDiv = 1,

 #if EIGEN_ARCH_ARM64 && !EIGEN_APPLE_DOUBLE_NEON_BUG
    HasExp = 1,
    HasLog = 1,
+    HasLog1p = 1,
+    HasExpm1 = 1,
    HasPow = 1,
    HasATan = 1,
    HasATanh = 1,
@@ -5390,13 +5372,11 @@ struct packet_traits<Eigen::half> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasAbsDiff = 0,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 1,
-    HasBlend = 0,
    HasInsert = 1,
    HasReduxp = 1,
    HasDiv = 1,
--- a/Eigen/src/Core/arch/SSE/Complex.h
+++ b/Eigen/src/Core/arch/SSE/Complex.h
@@ -49,7 +49,6 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
    HasMin = 0,
    HasMax = 0,
    HasSetLinear = 0,
-    HasBlend = 1
  };
 };
 #endif
@@ -413,13 +412,6 @@ EIGEN_STRONG_INLINE Packet1cd pcmp_eq(const Packet1cd& a, const Packet1cd& b) {
  return Packet1cd(pand<Packet2d>(eq, vec2d_swizzle1(eq, 1, 0)));
 }

-template <>
-EIGEN_STRONG_INLINE Packet2cf pblend(const Selector<2>& ifPacket, const Packet2cf& thenPacket,
-                                     const Packet2cf& elsePacket) {
-  __m128d result = pblend<Packet2d>(ifPacket, _mm_castps_pd(thenPacket.v), _mm_castps_pd(elsePacket.v));
-  return Packet2cf(_mm_castpd_ps(result));
-}
-
 template <>
 EIGEN_STRONG_INLINE Packet1cd psqrt<Packet1cd>(const Packet1cd& a) {
  return psqrt_complex<Packet1cd>(a);
--- a/Eigen/src/Core/arch/SSE/PacketMath.h
+++ b/Eigen/src/Core/arch/SSE/PacketMath.h
@@ -218,10 +218,12 @@ struct packet_traits<double> : default_packet_traits {
    HasSin = EIGEN_FAST_MATH,
    HasCos = EIGEN_FAST_MATH,
    HasTanh = EIGEN_FAST_MATH,
-    HasLog = 1,
    HasErf = EIGEN_FAST_MATH,
    HasErfc = EIGEN_FAST_MATH,
+    HasLog = 1,
    HasExp = 1,
+    HasLog1p = 1,
+    HasExpm1 = 1,
    HasPow = 1,
    HasSqrt = 1,
    HasRsqrt = 1,
@@ -290,7 +292,6 @@ struct packet_traits<bool> : default_packet_traits {
    HasCmp = 1,
    HasShift = 0,
    HasAbs = 0,
-    HasAbs2 = 0,
    HasMin = 0,
    HasMax = 0,
    HasConj = 0,
@@ -1998,44 +1999,6 @@ EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet16b, 16>& kernel) {
  kernel.packet[15] = _mm_unpackhi_epi64(u7, uf);
 }

-EIGEN_STRONG_INLINE __m128i sse_blend_mask(const Selector<2>& ifPacket) {
-  return _mm_set_epi64x(0 - ifPacket.select[1], 0 - ifPacket.select[0]);
-}
-
-EIGEN_STRONG_INLINE __m128i sse_blend_mask(const Selector<4>& ifPacket) {
-  return _mm_set_epi32(0 - ifPacket.select[3], 0 - ifPacket.select[2], 0 - ifPacket.select[1], 0 - ifPacket.select[0]);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet2l pblend(const Selector<2>& ifPacket, const Packet2l& thenPacket,
-                                    const Packet2l& elsePacket) {
-  const __m128i true_mask = sse_blend_mask(ifPacket);
-  return pselect<Packet2l>(true_mask, thenPacket, elsePacket);
-}
-template <>
-EIGEN_STRONG_INLINE Packet4i pblend(const Selector<4>& ifPacket, const Packet4i& thenPacket,
-                                    const Packet4i& elsePacket) {
-  const __m128i true_mask = sse_blend_mask(ifPacket);
-  return pselect<Packet4i>(true_mask, thenPacket, elsePacket);
-}
-template <>
-EIGEN_STRONG_INLINE Packet4ui pblend(const Selector<4>& ifPacket, const Packet4ui& thenPacket,
-                                     const Packet4ui& elsePacket) {
-  return (Packet4ui)pblend(ifPacket, (Packet4i)thenPacket, (Packet4i)elsePacket);
-}
-template <>
-EIGEN_STRONG_INLINE Packet4f pblend(const Selector<4>& ifPacket, const Packet4f& thenPacket,
-                                    const Packet4f& elsePacket) {
-  const __m128i true_mask = sse_blend_mask(ifPacket);
-  return pselect<Packet4f>(_mm_castsi128_ps(true_mask), thenPacket, elsePacket);
-}
-template <>
-EIGEN_STRONG_INLINE Packet2d pblend(const Selector<2>& ifPacket, const Packet2d& thenPacket,
-                                    const Packet2d& elsePacket) {
-  const __m128i true_mask = sse_blend_mask(ifPacket);
-  return pselect<Packet2d>(_mm_castsi128_pd(true_mask), thenPacket, elsePacket);
-}
-
 // Scalar path for pmadd with FMA to ensure consistency with vectorized path.
 #if defined(EIGEN_VECTORIZE_FMA)
 template <>
@@ -2194,7 +2157,6 @@ struct packet_traits<Eigen::half> : default_packet_traits {
    HasDiv    = 1,
    HasNegate = 0,
    HasAbs    = 0,
-    HasAbs2   = 0,
    HasMin    = 0,
    HasMax    = 0,
    HasConj   = 0,
--- a/Eigen/src/Core/arch/SVE/PacketMath.h
+++ b/Eigen/src/Core/arch/SVE/PacketMath.h
@@ -49,12 +49,10 @@ struct packet_traits<numext::int32_t> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 0,
-    HasBlend = 0,
    HasReduxp = 0  // Not implemented in SVE
  };
 };
@@ -344,12 +342,10 @@ struct packet_traits<float> : default_packet_traits {
    HasNegate = 1,
    HasAbs = 1,
    HasArg = 0,
-    HasAbs2 = 1,
    HasMin = 1,
    HasMax = 1,
    HasConj = 1,
    HasSetLinear = 0,
-    HasBlend = 0,
    HasReduxp = 0,  // Not implemented in SVE

    HasDiv = 1,
--- a/Eigen/src/Core/arch/SYCL/PacketMath.h
+++ b/Eigen/src/Core/arch/SYCL/PacketMath.h
@@ -542,31 +542,6 @@ EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void ptranspose(PacketBlock<cl::sycl::cl_d
  kernel.packet[1].x() = tmp;
 }

-template <>
-EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE cl::sycl::cl_half8 pblend(
-    const Selector<unpacket_traits<cl::sycl::cl_half8>::size>& ifPacket, const cl::sycl::cl_half8& thenPacket,
-    const cl::sycl::cl_half8& elsePacket) {
-  cl::sycl::cl_short8 condition(ifPacket.select[0] ? 0 : -1, ifPacket.select[1] ? 0 : -1, ifPacket.select[2] ? 0 : -1,
-                                ifPacket.select[3] ? 0 : -1, ifPacket.select[4] ? 0 : -1, ifPacket.select[5] ? 0 : -1,
-                                ifPacket.select[6] ? 0 : -1, ifPacket.select[7] ? 0 : -1);
-  return cl::sycl::select(thenPacket, elsePacket, condition);
-}
-
-template <>
-EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE cl::sycl::cl_float4 pblend(
-    const Selector<unpacket_traits<cl::sycl::cl_float4>::size>& ifPacket, const cl::sycl::cl_float4& thenPacket,
-    const cl::sycl::cl_float4& elsePacket) {
-  cl::sycl::cl_int4 condition(ifPacket.select[0] ? 0 : -1, ifPacket.select[1] ? 0 : -1, ifPacket.select[2] ? 0 : -1,
-                              ifPacket.select[3] ? 0 : -1);
-  return cl::sycl::select(thenPacket, elsePacket, condition);
-}
-
-template <>
-inline cl::sycl::cl_double2 pblend(const Selector<unpacket_traits<cl::sycl::cl_double2>::size>& ifPacket,
-                                   const cl::sycl::cl_double2& thenPacket, const cl::sycl::cl_double2& elsePacket) {
-  cl::sycl::cl_long2 condition(ifPacket.select[0] ? 0 : -1, ifPacket.select[1] ? 0 : -1);
-  return cl::sycl::select(thenPacket, elsePacket, condition);
-}
 #endif  // SYCL_DEVICE_ONLY

 }  // end namespace internal
--- a/Eigen/src/Core/arch/ZVector/Complex.h
+++ b/Eigen/src/Core/arch/ZVector/Complex.h
@@ -72,7 +72,6 @@ struct packet_traits<std::complex<float> > : default_packet_traits {
    HasAbs2 = 0,
    HasMin = 0,
    HasMax = 0,
-    HasBlend = 1,
    HasSetLinear = 0
  };
 };
@@ -469,14 +468,6 @@ EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet2cf, 2>& kernel) {
  kernel.packet[1].cd[0] = tmp;
 }

-template <>
-EIGEN_STRONG_INLINE Packet2cf pblend(const Selector<2>& ifPacket, const Packet2cf& thenPacket,
-                                     const Packet2cf& elsePacket) {
-  Packet2cf result;
-  const Selector<4> ifPacket4 = {ifPacket.select[0], ifPacket.select[0], ifPacket.select[1], ifPacket.select[1]};
-  result.v = pblend<Packet4f>(ifPacket4, thenPacket.v, elsePacket.v);
-  return result;
-}
 #else
 template <>
 EIGEN_STRONG_INLINE Packet2cf pcmp_eq(const Packet2cf& a, const Packet2cf& b) {
@@ -553,14 +544,6 @@ EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet2cf, 2>& kernel) {
  kernel.packet[0].v = tmp;
 }

-template <>
-EIGEN_STRONG_INLINE Packet2cf pblend(const Selector<2>& ifPacket, const Packet2cf& thenPacket,
-                                     const Packet2cf& elsePacket) {
-  Packet2cf result;
-  result.v = reinterpret_cast<Packet4f>(
-      pblend<Packet2d>(ifPacket, reinterpret_cast<Packet2d>(thenPacket.v), reinterpret_cast<Packet2d>(elsePacket.v)));
-  return result;
-}
 #endif

 }  // end namespace internal
--- a/Eigen/src/Core/arch/ZVector/PacketMath.h
+++ b/Eigen/src/Core/arch/ZVector/PacketMath.h
@@ -167,7 +167,6 @@ struct packet_traits<int> : default_packet_traits {
    HasSub = 1,
    HasMul = 1,
    HasDiv = 1,
-    HasBlend = 1
  };
 };

@@ -197,7 +196,6 @@ struct packet_traits<float> : default_packet_traits {
    HasTanh = 1,
    HasErf = 1,
    HasNegate = 1,
-    HasBlend = 1
  };
 };

@@ -224,7 +222,6 @@ struct packet_traits<double> : default_packet_traits {
    HasSqrt = 1,
    HasRsqrt = 1,
    HasNegate = 1,
-    HasBlend = 1
  };
 };

@@ -594,41 +591,32 @@ EIGEN_STRONG_INLINE void prefetch<double>(const double* addr) {

 template <int N>
 EIGEN_STRONG_INLINE Packet2l parithmetic_shift_right(const Packet2l& a) {
-  return Packet2l { parithmetic_shift_right<N>(a[0]), parithmetic_shift_right<N>(a[1]) };
+  return Packet2l{parithmetic_shift_right<N>(a[0]), parithmetic_shift_right<N>(a[1])};
 }
 template <int N>
 EIGEN_STRONG_INLINE Packet4i parithmetic_shift_right(const Packet4i& a) {
-  return Packet4i {
-    parithmetic_shift_right<N>(a[0]),
-    parithmetic_shift_right<N>(a[1]),
-    parithmetic_shift_right<N>(a[2]),
-    parithmetic_shift_right<N>(a[3]) };
+  return Packet4i{parithmetic_shift_right<N>(a[0]), parithmetic_shift_right<N>(a[1]), parithmetic_shift_right<N>(a[2]),
+                  parithmetic_shift_right<N>(a[3])};
 }

 template <int N>
 EIGEN_STRONG_INLINE Packet2l plogical_shift_right(const Packet2l& a) {
-  return Packet2l { plogical_shift_right<N>(a[0]), plogical_shift_right<N>(a[1]) };
+  return Packet2l{plogical_shift_right<N>(a[0]), plogical_shift_right<N>(a[1])};
 }
 template <int N>
 EIGEN_STRONG_INLINE Packet4i plogical_shift_right(const Packet4i& a) {
-  return Packet4i {
-    plogical_shift_right<N>(a[0]),
-    plogical_shift_right<N>(a[1]),
-    plogical_shift_right<N>(a[2]),
-    plogical_shift_right<N>(a[3]) };
+  return Packet4i{plogical_shift_right<N>(a[0]), plogical_shift_right<N>(a[1]), plogical_shift_right<N>(a[2]),
+                  plogical_shift_right<N>(a[3])};
 }

 template <int N>
 EIGEN_STRONG_INLINE Packet2l plogical_shift_left(const Packet2l& a) {
-  return Packet2l { plogical_shift_left<N>(a[0]), plogical_shift_left<N>(a[1]) };
+  return Packet2l{plogical_shift_left<N>(a[0]), plogical_shift_left<N>(a[1])};
 }
 template <int N>
 EIGEN_STRONG_INLINE Packet4i plogical_shift_left(const Packet4i& a) {
-  return Packet4i {
-    plogical_shift_left<N>(a[0]),
-    plogical_shift_left<N>(a[1]),
-    plogical_shift_left<N>(a[2]),
-    plogical_shift_left<N>(a[3]) };
+  return Packet4i{plogical_shift_left<N>(a[0]), plogical_shift_left<N>(a[1]), plogical_shift_left<N>(a[2]),
+                  plogical_shift_left<N>(a[3])};
 }

 template <>
@@ -747,22 +735,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet2d, 2>& kernel) {
  kernel.packet[1] = t1;
 }

-template <>
-EIGEN_STRONG_INLINE Packet4i pblend(const Selector<4>& ifPacket, const Packet4i& thenPacket,
-                                    const Packet4i& elsePacket) {
-  Packet4ui select = {ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3]};
-  Packet4ui mask = vec_cmpeq(select, reinterpret_cast<Packet4ui>(p4i_ONE));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
-template <>
-EIGEN_STRONG_INLINE Packet2d pblend(const Selector<2>& ifPacket, const Packet2d& thenPacket,
-                                    const Packet2d& elsePacket) {
-  Packet2ul select = {ifPacket.select[0], ifPacket.select[1]};
-  Packet2ul mask = vec_cmpeq(select, reinterpret_cast<Packet2ul>(p2l_ONE));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
 /* z13 has no vector float support so we emulate that with double
   z14 has proper vector float support.
 */
@@ -1068,19 +1040,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet4f, 4>& kernel) {
  kernel.packet[3].v4f[1] = t3.packet[1];
 }

-template <>
-EIGEN_STRONG_INLINE Packet4f pblend(const Selector<4>& ifPacket, const Packet4f& thenPacket,
-                                    const Packet4f& elsePacket) {
-  Packet2ul select_hi = {ifPacket.select[0], ifPacket.select[1]};
-  Packet2ul select_lo = {ifPacket.select[2], ifPacket.select[3]};
-  Packet2ul mask_hi = vec_cmpeq(select_hi, reinterpret_cast<Packet2ul>(p2l_ONE));
-  Packet2ul mask_lo = vec_cmpeq(select_lo, reinterpret_cast<Packet2ul>(p2l_ONE));
-  Packet4f result;
-  result.v4f[0] = vec_sel(elsePacket.v4f[0], thenPacket.v4f[0], mask_hi);
-  result.v4f[1] = vec_sel(elsePacket.v4f[1], thenPacket.v4f[1], mask_lo);
-  return result;
-}
-
 template <>
 Packet4f EIGEN_STRONG_INLINE pcmp_le<Packet4f>(const Packet4f& a, const Packet4f& b) {
  Packet4f res;
@@ -1288,14 +1247,6 @@ EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet4f, 4>& kernel) {
  kernel.packet[3] = vec_mergel(t1, t3);
 }

-template <>
-EIGEN_STRONG_INLINE Packet4f pblend(const Selector<4>& ifPacket, const Packet4f& thenPacket,
-                                    const Packet4f& elsePacket) {
-  Packet4ui select = {ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3]};
-  Packet4ui mask = vec_cmpeq(select, reinterpret_cast<Packet4ui>(p4i_ONE));
-  return vec_sel(elsePacket, thenPacket, mask);
-}
-
 #endif

 template <>
@@ -1338,62 +1289,51 @@ EIGEN_STRONG_INLINE Packet4f plset<Packet4f>(const float& a) {
 }

 #if !defined(vec_float) || !defined(__ARCH__) || (defined(__ARCH__) && __ARCH__ < 13)
-#pragma GCC warning \
-    "float->int and int->float conversion is simulated. compile for z15 for improved performance"
+#pragma GCC warning "float->int and int->float conversion is simulated. compile for z15 for improved performance"
 template <>
 struct cast_impl<Packet4i, Packet4f> {
  EIGEN_DEVICE_FUNC static inline Packet4f run(const Packet4i& a) {
-    return Packet4f{float(a[0]), float(a[1]), float(a[2]), float(a[3]) };
+    return Packet4f{float(a[0]), float(a[1]), float(a[2]), float(a[3])};
  }
 };

 template <>
 struct cast_impl<Packet4f, Packet4i> {
  EIGEN_DEVICE_FUNC static inline Packet4i run(const Packet4f& a) {
-    return Packet4i{int(a[0]), int(a[1]), int(a[2]), int(a[3]) };
+    return Packet4i{int(a[0]), int(a[1]), int(a[2]), int(a[3])};
  }
 };

 template <>
 struct cast_impl<Packet2l, Packet2d> {
-  EIGEN_DEVICE_FUNC static inline Packet2d run(const Packet2l& a) {
-    return Packet2d{double(a[0]), double(a[1]) };
-  }
+  EIGEN_DEVICE_FUNC static inline Packet2d run(const Packet2l& a) { return Packet2d{double(a[0]), double(a[1])}; }
 };

 template <>
 struct cast_impl<Packet2d, Packet2l> {
  EIGEN_DEVICE_FUNC static inline Packet2l run(const Packet2d& a) {
-    return Packet2l{(long long)(a[0]), (long long)(a[1]) };
+    return Packet2l{(long long)(a[0]), (long long)(a[1])};
  }
 };
 #else
 template <>
 struct cast_impl<Packet4i, Packet4f> {
-  EIGEN_DEVICE_FUNC static inline Packet4f run(const Packet4i& a) {
-    return vec_float(a);
-  }
+  EIGEN_DEVICE_FUNC static inline Packet4f run(const Packet4i& a) { return vec_float(a); }
 };

 template <>
 struct cast_impl<Packet4f, Packet4i> {
-  EIGEN_DEVICE_FUNC static inline Packet4i run(const Packet4f& a) {
-    return vec_signed(a);
-  }
+  EIGEN_DEVICE_FUNC static inline Packet4i run(const Packet4f& a) { return vec_signed(a); }
 };

 template <>
 struct cast_impl<Packet2l, Packet2d> {
-  EIGEN_DEVICE_FUNC static inline Packet2d run(const Packet2l& a) {
-    return vec_double(a);
-  }
+  EIGEN_DEVICE_FUNC static inline Packet2d run(const Packet2l& a) { return vec_double(a); }
 };

 template <>
 struct cast_impl<Packet2d, Packet2l> {
-  EIGEN_DEVICE_FUNC static inline Packet2l run(const Packet2d& a) {
-    return vec_signed(a);
-  }
+  EIGEN_DEVICE_FUNC static inline Packet2l run(const Packet2d& a) { return vec_signed(a); }
 };
 #endif

--- a/Eigen/src/Core/functors/TernaryFunctors.h
+++ b/Eigen/src/Core/functors/TernaryFunctors.h
@@ -21,7 +21,8 @@ namespace internal {

 template <typename ThenScalar, typename ElseScalar, typename ConditionScalar>
 struct scalar_boolean_select_op {
-  static constexpr bool ThenElseAreSame = is_same<ThenScalar, ElseScalar>::value;
+  static constexpr bool ThenElseAreSame =
+      is_same<std::remove_const_t<ThenScalar>, std::remove_const_t<ElseScalar>>::value;
  EIGEN_STATIC_ASSERT(ThenElseAreSame, THEN AND ELSE MUST BE SAME TYPE)
  using Scalar = ThenScalar;
  using result_type = Scalar;
--- a/Eigen/src/Core/functors/UnaryFunctors.h
+++ b/Eigen/src/Core/functors/UnaryFunctors.h
@@ -100,7 +100,10 @@ struct scalar_abs2_op {
 };
 template <typename Scalar>
 struct functor_traits<scalar_abs2_op<Scalar>> {
-  enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasAbs2 };
+  enum {
+    Cost = NumTraits<Scalar>::MulCost,
+    PacketAccess = packet_traits<Scalar>::HasMul && !NumTraits<Scalar>::IsComplex
+  };
 };

 template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -830,6 +830,11 @@
 #endif
 #endif

+// Multidimensional subscript operator feature test
+#if defined(__cpp_multidimensional_subscript) && __cpp_multidimensional_subscript >= 202110L
+#define EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+#endif
+
 //------------------------------------------------------------------------------------------
 // Preprocessor programming helpers
 //------------------------------------------------------------------------------------------
--- a/Eigen/src/Eigenvalues/ComplexQZ.h
+++ b/Eigen/src/Eigenvalues/ComplexQZ.h
@@ -0,0 +1,651 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Alexey Korepanov
+// Copyright (C) 2025 Ludwig Striet <ludwig.striet@mathematik.uni-freiburg.de>
+//
+// 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
+// https://mozilla.org/MPL/2.0/.
+//
+// Derived from: Eigen/src/Eigenvalues/RealQZ.h
+
+#ifndef EIGEN_COMPLEX_QZ_H_
+#define EIGEN_COMPLEX_QZ_H_
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+ *
+ *
+ * \class ComplexQZ
+ *
+ * \brief Performs a QZ decomposition of a pair of matrices A, B
+ *
+ * \tparam MatrixType_ the type input type of the matrix.
+ *
+ * Given to complex square matrices A and B, this class computes the QZ decomposition
+ * \f$ A = Q S Z \f$, \f$ B = Q T Z\f$ where Q and Z are unitary matrices and
+ * S and T a re upper-triangular matrices. More precisely, Q and Z fulfill
+ * \f$ Q Q* = Id\f$ and \f$ Z Z* = Id\f$. The generalized Eigenvalues are then
+ * obtained as ratios of corresponding diagonal entries, lambda(i) = S(i,i) / T(i, i).
+ *
+ * The QZ algorithm was introduced in the seminal work "An Algorithm for
+ * Generalized Matrix Eigenvalue Problems" by Moler & Stewart in 1973. The matrix
+ * pair S = A, T = B is first transformed to Hessenberg-Triangular form where S is an
+ * upper Hessenberg matrix and T is an upper Triangular matrix.
+ *
+ * This pair is subsequently reduced to the desired form using implicit QZ shifts as
+ * described in the original paper. The algorithms to find small entries on the
+ * diagonals and subdiagonals are based on the variants in the implementation
+ * for Real matrices in the RealQZ class.
+ *
+ * \sa class RealQZ
+ */
+
+namespace Eigen {
+
+template <typename MatrixType_>
+class ComplexQZ {
+ public:
+  using MatrixType = MatrixType_;
+  using Scalar = typename MatrixType_::Scalar;
+  using RealScalar = typename MatrixType_::RealScalar;
+
+  enum {
+    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+    ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+    Options = internal::traits<MatrixType>::Options,
+    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+  };
+
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using Vec2 = Matrix<Scalar, 2, 1>;
+  using Vec3 = Matrix<Scalar, 3, 1>;
+  using Row2 = Matrix<Scalar, 1, 2>;
+  using Mat2 = Matrix<Scalar, 2, 2>;
+
+  /** \brief Returns matrix Q in the QZ decomposition.
+   *
+   * \returns A const reference to the matrix Q.
+   */
+  const MatrixType& matrixQ() const {
+    eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
+    eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
+    return m_Q;
+  }
+
+  /** \brief Returns matrix Z in the QZ decomposition.
+   *
+   * \returns A const reference to the matrix Z.
+   */
+  const MatrixType& matrixZ() const {
+    eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
+    eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
+    return m_Z;
+  }
+
+  /** \brief Returns matrix S in the QZ decomposition.
+   *
+   * \returns A const reference to the matrix S.
+   */
+  const MatrixType& matrixS() const {
+    eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
+    return m_S;
+  }
+
+  /** \brief Returns matrix S in the QZ decomposition.
+   *
+   * \returns A const reference to the matrix S.
+   */
+  const MatrixType& matrixT() const {
+    eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
+    return m_T;
+  }
+
+  /** \brief Constructor
+   *
+   * \param[in] n size of the matrices whose QZ decomposition we compute
+   *
+   * This constructor is used when we use the compute(...) method later,
+   * especially when we aim to compute the decomposition of two sparse
+   * matrices.
+   */
+  ComplexQZ(Index n, bool computeQZ = true, unsigned int maxIters = 400)
+      : m_n(n),
+        m_S(n, n),
+        m_T(n, n),
+        m_Q(computeQZ ? n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
+            computeQZ ? n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
+        m_Z(computeQZ ? n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
+            computeQZ ? n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
+        m_ws(2 * n),
+        m_computeQZ(computeQZ),
+        m_maxIters(maxIters) {}
+
+  /** \brief Constructor. computes the QZ decomposition of given matrices
+   * upon creation
+   *
+   * \param[in] A         input matrix A
+   * \param[in] B         input matrix B
+   * \param[in] computeQZ If false, the matrices Q and Z are not computed
+   *
+   * This constructor calls the compute() method to compute the QZ decomposition.
+   * If input matrices are sparse, call the constructor that uses only the
+   * size as input the computeSparse(...) method.
+   */
+  ComplexQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true, unsigned int maxIters = 400)
+      : m_n(A.rows()),
+        m_maxIters(maxIters),
+        m_computeQZ(computeQZ),
+        m_S(A.rows(), A.cols()),
+        m_T(A.rows(), A.cols()),
+        m_Q(computeQZ ? m_n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
+            computeQZ ? m_n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
+        m_Z(computeQZ ? m_n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
+            computeQZ ? m_n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
+        m_ws(2 * m_n) {
+    compute(A, B, computeQZ);
+  }
+
+  /** \brief Compute the QZ decomposition of complex input matrices
+   *
+   * \param[in] A         Matrix A.
+   * \param[in] B         Matrix B.
+   * \param[in] computeQZ If false, the matrices Q and Z are not computed.
+   */
+  void compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true);
+
+  /** \brief Compute the decomposition of sparse complex input matrices.
+   * Main difference to the compute(...) method is that it computes a
+   * SparseQR decomposition of B
+   *
+   * \param[in] A         Matrix A.
+   * \param[in] B         Matrix B.
+   * \param[in] computeQZ If false, the matrices Q and Z are not computed.
+   */
+  template <typename SparseMatrixType_>
+  void computeSparse(const SparseMatrixType_& A, const SparseMatrixType_& B, bool computeQZ = true);
+
+  /** \brief Reports whether the last computation was successfull.
+   *
+   * \returns \c Success if computation was successfull, \c NoConvergence otherwise.
+   */
+  ComputationInfo info() const { return m_info; }
+
+  /** \brief number of performed QZ steps
+   */
+  unsigned int iterations() const {
+    eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
+    return m_global_iter;
+  }
+
+ private:
+  Index m_n;
+  const unsigned int m_maxIters;
+  unsigned int m_global_iter;
+  bool m_isInitialized;
+  bool m_computeQZ;
+  ComputationInfo m_info;
+  MatrixType m_S, m_T, m_Q, m_Z;
+  RealScalar m_normOfT, m_normOfS;
+  Vec m_ws;
+
+  // Test if a Scalar is 0 up to a certain tolerance
+  static bool is_negligible(const Scalar x, const RealScalar tol = NumTraits<RealScalar>::epsilon()) {
+    return numext::abs(x) <= tol;
+  }
+
+  void do_QZ_step(Index p, Index q);
+
+  inline Mat2 computeZk2(const Row2& b);
+
+  // This is basically taken from from Eigen3::RealQZ
+  void hessenbergTriangular(const MatrixType& A, const MatrixType& B);
+
+  // This function can be called when m_Q and m_Z are initialized and m_S, m_T
+  // are in hessenberg-triangular form
+  void reduceHessenbergTriangular();
+
+  // Sparse variant of the above method.
+  template <typename SparseMatrixType_>
+  void hessenbergTriangularSparse(const SparseMatrixType_& A, const SparseMatrixType_& B);
+
+  void computeNorms();
+
+  Index findSmallSubdiagEntry(Index l);
+  Index findSmallDiagEntry(Index f, Index l);
+
+  void push_down_zero_ST(Index k, Index l);
+
+  void reduceDiagonal2x2block(Index i);
+};
+
+template <typename MatrixType_>
+void ComplexQZ<MatrixType_>::compute(const MatrixType& A, const MatrixType& B, bool computeQZ) {
+  m_computeQZ = computeQZ;
+  m_n = A.rows();
+
+  eigen_assert(m_n == A.cols() && "A is not a square matrix");
+  eigen_assert(m_n == B.rows() && m_n == B.cols() && "B is not a square matrix or B is not of the same size as A");
+
+  m_isInitialized = true;
+  m_global_iter = 0;
+
+  // This will initialize m_Q and m_Z and bring m_S, m_T to hessenberg-triangular form
+  hessenbergTriangular(A, B);
+
+  // We assume that we already have that S is upper-Hessenberg and T is
+  // upper-triangular. This is what the hessenbergTriangular(...) method does
+  reduceHessenbergTriangular();
+}
+
+// This is basically taken from from Eigen3::RealQZ
+template <typename MatrixType_>
+void ComplexQZ<MatrixType_>::hessenbergTriangular(const MatrixType& A, const MatrixType& B) {
+  // Copy A and B, these will be the matrices on which we operate later
+  m_S = A;
+  m_T = B;
+
+  // Perform QR decomposition of the matrix Q
+  HouseholderQR<MatrixType> qr(m_T);
+  m_T = qr.matrixQR();
+  m_T.template triangularView<StrictlyLower>().setZero();
+
+  if (m_computeQZ) m_Q = qr.householderQ();
+
+  // overwrite S with Q* x S
+  m_S.applyOnTheLeft(qr.householderQ().adjoint());
+
+  if (m_computeQZ) m_Z = MatrixType::Identity(m_n, m_n);
+
+  // reduce S to upper Hessenberg with Givens rotations
+  for (Index j = 0; j <= m_n - 3; j++) {
+    for (Index i = m_n - 1; i >= j + 2; i--) {
+      JacobiRotation<Scalar> G;
+      // delete S(i,j)
+      if (!numext::is_exactly_zero(m_S.coeff(i, j))) {
+        G.makeGivens(m_S.coeff(i - 1, j), m_S.coeff(i, j), &m_S.coeffRef(i - 1, j));
+        m_S.coeffRef(i, j) = Scalar(0);
+        m_T.rightCols(m_n - i + 1).applyOnTheLeft(i - 1, i, G.adjoint());
+        m_S.rightCols(m_n - j - 1).applyOnTheLeft(i - 1, i, G.adjoint());
+        // This is what we want to achieve
+        if (!is_negligible(m_S(i, j)))
+          m_info = ComputationInfo::NumericalIssue;
+        else
+          m_S(i, j) = Scalar(0);
+        // update Q
+        if (m_computeQZ) m_Q.applyOnTheRight(i - 1, i, G);
+      }
+
+      if (!numext::is_exactly_zero(m_T.coeff(i, i - 1))) {
+        // Compute rotation and update matrix T
+        G.makeGivens(m_T.coeff(i, i), m_T.coeff(i, i - 1), &m_T.coeffRef(i, i));
+        m_T.topRows(i).applyOnTheRight(i - 1, i, G.adjoint());
+        m_T.coeffRef(i, i - 1) = Scalar(0);
+        // Update matrix S
+        m_S.applyOnTheRight(i - 1, i, G.adjoint());
+        // update Z
+        if (m_computeQZ) m_Z.applyOnTheLeft(i - 1, i, G);
+      }
+    }
+  }
+}
+
+template <typename MatrixType>
+template <typename SparseMatrixType_>
+void ComplexQZ<MatrixType>::hessenbergTriangularSparse(const SparseMatrixType_& A, const SparseMatrixType_& B) {
+  m_S = A.toDense();
+
+  SparseQR<SparseMatrix<Scalar, ColMajor>, NaturalOrdering<Index>> sparseQR;
+
+  eigen_assert(B.isCompressed() &&
+               "SparseQR requires a sparse matrix in compressed mode."
+               "Call .makeCompressed() before passing it to SparseQR");
+
+  // Computing QR decomposition of T...
+  sparseQR.setPivotThreshold(RealScalar(0));  // This prevends algorithm from doing pivoting
+  sparseQR.compute(B);
+  // perform QR decomposition of T, overwrite T with R, save Q
+  // HouseholderQR<Mat> qrT(m_T);
+  m_T = sparseQR.matrixR();
+  m_T.template triangularView<StrictlyLower>().setZero();
+
+  if (m_computeQZ) m_Q = sparseQR.matrixQ();
+
+  // overwrite S with Q* S
+  m_S = sparseQR.matrixQ().adjoint() * m_S;
+
+  if (m_computeQZ) m_Z = MatrixType::Identity(m_n, m_n);
+
+  // reduce S to upper Hessenberg with Givens rotations
+  for (Index j = 0; j <= m_n - 3; j++) {
+    for (Index i = m_n - 1; i >= j + 2; i--) {
+      JacobiRotation<Scalar> G;
+      // kill S(i,j)
+      // if(!numext::is_exactly_zero(_S.coeff(i, j)))
+      if (m_S.coeff(i, j) != Scalar(0)) {
+        // This is the adapted code
+        G.makeGivens(m_S.coeff(i - 1, j), m_S.coeff(i, j), &m_S.coeffRef(i - 1, j));
+        m_S.coeffRef(i, j) = Scalar(0);
+        m_T.rightCols(m_n - i + 1).applyOnTheLeft(i - 1, i, G.adjoint());
+        m_S.rightCols(m_n - j - 1).applyOnTheLeft(i - 1, i, G.adjoint());
+        // This is what we want to achieve
+        if (!is_negligible(m_S(i, j))) {
+          m_info = ComputationInfo::NumericalIssue;
+        }
+        m_S(i, j) = Scalar(0);
+        // update Q
+        if (m_computeQZ) m_Q.applyOnTheRight(i - 1, i, G);
+      }
+
+      if (!numext::is_exactly_zero(m_T.coeff(i, i - 1))) {
+        // Compute rotation and update matrix T
+        G.makeGivens(m_T.coeff(i, i), m_T.coeff(i, i - 1), &m_T.coeffRef(i, i));
+        m_T.topRows(i).applyOnTheRight(i - 1, i, G.adjoint());
+        m_T.coeffRef(i, i - 1) = Scalar(0);
+        // Update matrix S
+        m_S.applyOnTheRight(i - 1, i, G.adjoint());
+        // update Z
+        if (m_computeQZ) m_Z.applyOnTheLeft(i - 1, i, G);
+      }
+    }
+  }
+}
+
+template <typename MatrixType>
+template <typename SparseMatrixType_>
+void ComplexQZ<MatrixType>::computeSparse(const SparseMatrixType_& A, const SparseMatrixType_& B, bool computeQZ) {
+  m_computeQZ = computeQZ;
+  m_n = A.rows();
+  eigen_assert(m_n == A.cols() && "A is not a square matrix");
+  eigen_assert(m_n == B.rows() && m_n == B.cols() && "B is not a square matrix or B is not of the same size as A");
+  m_isInitialized = true;
+  m_global_iter = 0;
+  hessenbergTriangularSparse(A, B);
+
+  // We assume that we already have that A is upper-Hessenberg and B is
+  // upper-triangular. This is what the hessenbergTriangular(...) method does
+  reduceHessenbergTriangular();
+}
+
+template <typename MatrixType_>
+void ComplexQZ<MatrixType_>::reduceHessenbergTriangular() {
+  Index l = m_n - 1, f;
+  unsigned int local_iter = 0;
+  computeNorms();
+
+  while (l > 0 && local_iter < m_maxIters) {
+    f = findSmallSubdiagEntry(l);
+
+    // Subdiag entry is small -> can be safely set to 0
+    if (f > 0) {
+      m_S.coeffRef(f, f - 1) = Scalar(0);
+    }
+    if (f == l) {  // One root found
+      l--;
+      local_iter = 0;
+    } else if (f == l - 1) {  // Two roots found
+      // We found an undesired non-zero at (f+1,f) in S and eliminate it immediately
+      reduceDiagonal2x2block(f);
+      l -= 2;
+      local_iter = 0;
+    } else {
+      Index z = findSmallDiagEntry(f, l);
+      if (z >= f) {
+        push_down_zero_ST(z, l);
+      } else {
+        do_QZ_step(f, m_n - l - 1);
+        local_iter++;
+        m_global_iter++;
+      }
+    }
+  }
+
+  m_info = (local_iter < m_maxIters) ? Success : NoConvergence;
+}
+
+template <typename MatrixType_>
+inline typename ComplexQZ<MatrixType_>::Mat2 ComplexQZ<MatrixType_>::computeZk2(const Row2& b) {
+  Mat2 S;
+  S << Scalar(0), Scalar(1), Scalar(1), Scalar(0);
+  Vec2 bprime = S * b.adjoint();
+  JacobiRotation<Scalar> J;
+  J.makeGivens(bprime(0), bprime(1));
+  Mat2 Z = S;
+  Z.applyOnTheLeft(0, 1, J);
+  Z = S * Z;
+  return Z;
+}
+
+template <typename MatrixType_>
+void ComplexQZ<MatrixType_>::do_QZ_step(Index p, Index q) {
+  // This is certainly not the most efficient way of doing this,
+  // but a readable one.
+  const auto a = [p, this](Index i, Index j) { return m_S(p + i - 1, p + j - 1); };
+  const auto b = [p, this](Index i, Index j) { return m_T(p + i - 1, p + j - 1); };
+  const Index m = m_n - p - q;  // Size of the inner block
+  Scalar x, y, z;
+  // We could introduce doing exceptional shifts from time to time.
+  Scalar W1 = a(m - 1, m - 1) / b(m - 1, m - 1) - a(1, 1) / b(1, 1), W2 = a(m, m) / b(m, m) - a(1, 1) / b(1, 1),
+         W3 = a(m, m - 1) / b(m - 1, m - 1);
+
+  x = (W1 * W2 - a(m - 1, m) / b(m, m) * W3 + W3 * b(m - 1, m) / b(m, m) * a(1, 1) / b(1, 1)) * b(1, 1) / a(2, 1) +
+      a(1, 2) / b(2, 2) - a(1, 1) / b(1, 1) * b(1, 2) / b(2, 2);
+  y = (a(2, 2) / b(2, 2) - a(1, 1) / b(1, 1)) - a(2, 1) / b(1, 1) * b(1, 2) / b(2, 2) - W1 - W2 +
+      W3 * (b(m - 1, m) / b(m, m));
+  z = a(3, 2) / b(2, 2);
+  Vec3 X;
+  const PermutationMatrix<3, 3, int> S3(Vector3i(2, 0, 1));
+  for (Index k = p; k < p + m - 2; k++) {
+    X << x, y, z;
+    Vec2 ess;
+    Scalar tau;
+    RealScalar beta;
+    X.makeHouseholder(ess, tau, beta);
+    // The permutations are needed because the makeHouseHolder-method computes
+    // the householder transformation in a way that the vector is reflected to
+    // (1 0 ... 0) instead of (0 ... 0 1)
+    m_S.template middleRows<3>(k)
+        .rightCols((std::min)(m_n, m_n - k + 1))
+        .applyHouseholderOnTheLeft(ess, tau, m_ws.data());
+    m_T.template middleRows<3>(k).rightCols(m_n - k).applyHouseholderOnTheLeft(ess, tau, m_ws.data());
+    if (m_computeQZ) m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(ess, std::conj(tau), m_ws.data());
+
+    // Compute Matrix Zk1 s.t. (b(k+2,k) ... b(k+2, k+2)) Zk1 = (0,0,*)
+    Vec3 bprime = (m_T.template block<1, 3>(k + 2, k) * S3).adjoint();
+    bprime.makeHouseholder(ess, tau, beta);
+    m_S.template middleCols<3>(k).topRows((std::min)(k + 4, m_n)).applyOnTheRight(S3);
+    m_S.template middleCols<3>(k)
+        .topRows((std::min)(k + 4, m_n))
+        .applyHouseholderOnTheRight(ess, std::conj(tau), m_ws.data());
+    m_S.template middleCols<3>(k).topRows((std::min)(k + 4, m_n)).applyOnTheRight(S3.transpose());
+    m_T.template middleCols<3>(k).topRows((std::min)(k + 3, m_n)).applyOnTheRight(S3);
+    m_T.template middleCols<3>(k)
+        .topRows((std::min)(k + 3, m_n))
+        .applyHouseholderOnTheRight(ess, std::conj(tau), m_ws.data());
+    m_T.template middleCols<3>(k).topRows((std::min)(k + 3, m_n)).applyOnTheRight(S3.transpose());
+    if (m_computeQZ) {
+      m_Z.template middleRows<3>(k).applyOnTheLeft(S3.transpose());
+      m_Z.template middleRows<3>(k).applyHouseholderOnTheLeft(ess, tau, m_ws.data());
+      m_Z.template middleRows<3>(k).applyOnTheLeft(S3);
+    }
+    Mat2 Zk2 = computeZk2(m_T.template block<1, 2>(k + 1, k));
+    m_S.template middleCols<2>(k).topRows((std::min)(k + 4, m_n)).applyOnTheRight(Zk2);
+    m_T.template middleCols<2>(k).topRows((std::min)(k + 3, m_n)).applyOnTheRight(Zk2);
+
+    if (m_computeQZ) m_Z.template middleRows<2>(k).applyOnTheLeft(Zk2.adjoint());
+
+    x = m_S(k + 1, k);
+    y = m_S(k + 2, k);
+    if (k < p + m - 3) {
+      z = m_S(k + 3, k);
+    }
+  };
+
+  // Find a Householdermartirx Qn1 s.t. Qn1 (x y)^T = (* 0)
+  JacobiRotation<Scalar> J;
+  J.makeGivens(x, y);
+  m_S.template middleRows<2>(p + m - 2).applyOnTheLeft(0, 1, J.adjoint());
+  m_T.template middleRows<2>(p + m - 2).applyOnTheLeft(0, 1, J.adjoint());
+
+  if (m_computeQZ) m_Q.template middleCols<2>(p + m - 2).applyOnTheRight(0, 1, J);
+
+  // Find a Householdermatrix Zn1 s.t. (b(n,n-1) b(n,n)) * Zn1 = (0 *)
+  Mat2 Zn1 = computeZk2(m_T.template block<1, 2>(p + m - 1, p + m - 2));
+  m_S.template middleCols<2>(p + m - 2).applyOnTheRight(Zn1);
+  m_T.template middleCols<2>(p + m - 2).applyOnTheRight(Zn1);
+
+  if (m_computeQZ) m_Z.template middleRows<2>(p + m - 2).applyOnTheLeft(Zn1.adjoint());
+}
+
+/** \internal we found an undesired non-zero at (i+1,i) on the subdiagonal of S and reduce the block */
+template <typename MatrixType_>
+void ComplexQZ<MatrixType_>::reduceDiagonal2x2block(Index i) {
+  // We have found a non-zero on the subdiagonal and want to eliminate it
+  Mat2 Si = m_S.template block<2, 2>(i, i), Ti = m_T.template block<2, 2>(i, i);
+  if (is_negligible(Ti(0, 0)) && !is_negligible(Ti(1, 1))) {
+    Eigen::JacobiRotation<Scalar> G;
+    G.makeGivens(m_S(i, i), m_S(i + 1, i));
+    m_S.applyOnTheLeft(i, i + 1, G.adjoint());
+    m_T.applyOnTheLeft(i, i + 1, G.adjoint());
+
+    if (m_computeQZ) m_Q.applyOnTheRight(i, i + 1, G);
+
+  } else if (!is_negligible(Ti(0, 0)) && is_negligible(Ti(1, 1))) {
+    Eigen::JacobiRotation<Scalar> G;
+    G.makeGivens(m_S(i + 1, i + 1), m_S(i + 1, i));
+    m_S.applyOnTheRight(i, i + 1, G.adjoint());
+    m_T.applyOnTheRight(i, i + 1, G.adjoint());
+    if (m_computeQZ) m_Z.applyOnTheLeft(i, i + 1, G);
+  } else if (!is_negligible(Ti(0, 0)) && !is_negligible((Ti(1, 1)))) {
+    Scalar mu = Si(0, 0) / Ti(0, 0);
+    Scalar a12_bar = Si(0, 1) - mu * Ti(0, 1);
+    Scalar a22_bar = Si(1, 1) - mu * Ti(1, 1);
+    Scalar p = Scalar(0.5) * (a22_bar / Ti(1, 1) - Ti(0, 1) * Si(1, 0) / (Ti(0, 0) * Ti(1, 1)));
+    RealScalar sgn_p = p.real() >= RealScalar(0) ? RealScalar(1) : RealScalar(-1);
+    Scalar q = Si(1, 0) * a12_bar / (Ti(0, 0) * Ti(1, 1));
+    Scalar r = p * p + q;
+    Scalar lambda = mu + p + sgn_p * numext::sqrt(r);
+    Mat2 E = Si - lambda * Ti;
+    Index l;
+    E.rowwise().norm().maxCoeff(&l);
+    JacobiRotation<Scalar> G;
+    G.makeGivens(E(l, 1), E(l, 0));
+    m_S.applyOnTheRight(i, i + 1, G.adjoint());
+    m_T.applyOnTheRight(i, i + 1, G.adjoint());
+
+    if (m_computeQZ) m_Z.applyOnTheLeft(i, i + 1, G);
+
+    Mat2 tildeSi = m_S.template block<2, 2>(i, i), tildeTi = m_T.template block<2, 2>(i, i);
+    Mat2 C = tildeSi.norm() < (lambda * tildeTi).norm() ? tildeSi : lambda * tildeTi;
+    G.makeGivens(C(0, 0), C(1, 0));
+    m_S.applyOnTheLeft(i, i + 1, G.adjoint());
+    m_T.applyOnTheLeft(i, i + 1, G.adjoint());
+
+    if (m_computeQZ) m_Q.applyOnTheRight(i, i + 1, G);
+  }
+
+  if (!is_negligible(m_S(i + 1, i), m_normOfS * NumTraits<RealScalar>::epsilon())) {
+    m_info = ComputationInfo::NumericalIssue;
+  } else {
+    m_S(i + 1, i) = Scalar(0);
+  }
+}
+
+/** \internal We found a zero at T(k,k) and want to "push it down" to T(l,l) */
+template <typename MatrixType_>
+void ComplexQZ<MatrixType_>::push_down_zero_ST(Index k, Index l) {
+  // Test Preconditions
+
+  JacobiRotation<Scalar> J;
+  for (Index j = k + 1; j <= l; j++) {
+    // Create a 0 at _T(j, j)
+    J.makeGivens(m_T(j - 1, j), m_T(j, j), &m_T.coeffRef(j - 1, j));
+    if (m_n - j - 1 > 0) {
+      m_T.rightCols(m_n - j - 1).applyOnTheLeft(j - 1, j, J.adjoint());
+    }
+    m_T.coeffRef(j, j) = Scalar(0);
+
+    m_S.applyOnTheLeft(j - 1, j, J.adjoint());
+
+    if (m_computeQZ) m_Q.applyOnTheRight(j - 1, j, J);
+
+    // Delete the non-desired non-zero at _S(j, j-2)
+    if (j > 1) {
+      J.makeGivens(std::conj(m_S(j, j - 1)), std::conj(m_S(j, j - 2)));
+      m_S.applyOnTheRight(j - 1, j - 2, J);
+      m_S(j, j - 2) = Scalar(0);
+      m_T.applyOnTheRight(j - 1, j - 2, J);
+      if (m_computeQZ) m_Z.applyOnTheLeft(j - 1, j - 2, J.adjoint());
+    }
+  }
+
+  // Assume we have the desired structure now, up to the non-zero entry at
+  // _S(l, l-1) which we will delete through a last right-jacobi-rotation
+  J.makeGivens(std::conj(m_S(l, l)), std::conj(m_S(l, l - 1)));
+  m_S.topRows(l + 1).applyOnTheRight(l, l - 1, J);
+
+  if (!is_negligible(m_S(l, l - 1), m_normOfS * NumTraits<Scalar>::epsilon())) {
+    m_info = ComputationInfo::NumericalIssue;
+  } else {
+    m_S(l, l - 1) = Scalar(0);
+  }
+  m_T.topRows(l + 1).applyOnTheRight(l, l - 1, J);
+
+  if (m_computeQZ) m_Z.applyOnTheLeft(l, l - 1, J.adjoint());
+
+  // Ensure postconditions
+  if (!is_negligible(m_T(l, l)) || !is_negligible(m_S(l, l - 1))) {
+    m_info = ComputationInfo::NumericalIssue;
+  } else {
+    m_T(l, l) = Scalar(0);
+    m_S(l, l - 1) = Scalar(0);
+  }
+}
+
+/** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */
+template <typename MatrixType_>
+void ComplexQZ<MatrixType_>::computeNorms() {
+  const Index size = m_S.cols();
+  m_normOfS = RealScalar(0);
+  m_normOfT = RealScalar(0);
+  for (Index j = 0; j < size; ++j) {
+    m_normOfS += m_S.col(j).segment(0, (std::min)(size, j + 2)).cwiseAbs().sum();
+    m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum();
+  }
+}
+
+/** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0). Copied from Eigen3 RealQZ
+ * implementation */
+template <typename MatrixType_>
+inline Index ComplexQZ<MatrixType_>::findSmallSubdiagEntry(Index iu) {
+  Index res = iu;
+  while (res > 0) {
+    RealScalar s = numext::abs(m_S.coeff(res - 1, res - 1)) + numext::abs(m_S.coeff(res, res));
+    if (s == Scalar(0)) s = m_normOfS;
+    if (numext::abs(m_S.coeff(res, res - 1)) < NumTraits<RealScalar>::epsilon() * s) break;
+    res--;
+  }
+  return res;
+}
+
+//
+/** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1).
+ * Copied from Eigen3 RealQZ implementation. */
+template <typename MatrixType_>
+inline Index ComplexQZ<MatrixType_>::findSmallDiagEntry(Index f, Index l) {
+  Index res = l;
+  while (res >= f) {
+    if (numext::abs(m_T.coeff(res, res)) <= NumTraits<RealScalar>::epsilon() * m_normOfT) break;
+    res--;
+  }
+  return res;
+}
+
+}  // namespace Eigen
+
+#endif  // _COMPLEX_QZ_H_
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -325,6 +325,22 @@ class SelfAdjointEigenSolver {
    return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
  }

+  /** \brief Computes the matrix exponential the matrix.
+   *
+   * \returns the matrix exponential the matrix.
+   *
+   * \pre The eigenvalues and eigenvectors of a positive-definite matrix
+   * have been computed before.
+   *
+   * \sa operatorInverseSqrt(), operatorSqrt(),
+   * <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a>
+   */
+  EIGEN_DEVICE_FUNC MatrixType operatorExp() const {
+    eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
+    eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
+    return m_eivec * m_eivalues.array().exp().matrix().asDiagonal() * m_eivec.adjoint();
+  }
+
  /** \brief Computes the inverse square root of the matrix.
   *
   * \returns the inverse positive-definite square root of the matrix
--- a/Eigen/src/Geometry/Transform.h
+++ b/Eigen/src/Geometry/Transform.h
@@ -365,6 +365,15 @@ class Transform {
   * \sa MatrixBase::operator(Index,Index) */
  EIGEN_DEVICE_FUNC inline Scalar& operator()(Index row, Index col) { return m_matrix(row, col); }

+#ifdef EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+  /** shortcut for m_matrix(row,col);
+   * \sa MatrixBase::operator(Index,Index) const */
+  EIGEN_DEVICE_FUNC inline Scalar operator[](Index row, Index col) const { return m_matrix[row, col]; }
+  /** shortcut for m_matrix(row,col);
+   * \sa MatrixBase::operator(Index,Index) */
+  EIGEN_DEVICE_FUNC inline Scalar& operator[](Index row, Index col) { return m_matrix[row, col]; }
+#endif
+
  /** \returns a read-only expression of the transformation matrix */
  EIGEN_DEVICE_FUNC inline const MatrixType& matrix() const { return m_matrix; }
  /** \returns a writable expression of the transformation matrix */
--- a/Eigen/src/Jacobi/Jacobi.h
+++ b/Eigen/src/Jacobi/Jacobi.h
@@ -335,8 +335,8 @@ struct apply_rotation_in_the_plane_selector<Scalar, OtherScalar, SizeAtCompileTi
        for (Index i = alignedStart; i < alignedEnd; i += PacketSize) {
          Packet xi = pload<Packet>(px);
          Packet yi = pload<Packet>(py);
-          pstore(px, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
-          pstore(py, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
+          pstore(px, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
+          pstore(py, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
          px += PacketSize;
          py += PacketSize;
        }
@@ -347,18 +347,18 @@ struct apply_rotation_in_the_plane_selector<Scalar, OtherScalar, SizeAtCompileTi
          Packet xi1 = ploadu<Packet>(px + PacketSize);
          Packet yi = pload<Packet>(py);
          Packet yi1 = pload<Packet>(py + PacketSize);
-          pstoreu(px, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
-          pstoreu(px + PacketSize, padd(pm.pmul(pc, xi1), pcj.pmul(ps, yi1)));
-          pstore(py, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
-          pstore(py + PacketSize, psub(pcj.pmul(pc, yi1), pm.pmul(ps, xi1)));
+          pstoreu(px, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
+          pstoreu(px + PacketSize, pm.pmadd(pc, xi1, pcj.pmul(ps, yi1)));
+          pstore(py, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
+          pstore(py + PacketSize, pcj.pmsub(pc, yi1, pm.pmul(ps, xi1)));
          px += Peeling * PacketSize;
          py += Peeling * PacketSize;
        }
        if (alignedEnd != peelingEnd) {
          Packet xi = ploadu<Packet>(x + peelingEnd);
          Packet yi = pload<Packet>(y + peelingEnd);
-          pstoreu(x + peelingEnd, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
-          pstore(y + peelingEnd, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
+          pstoreu(x + peelingEnd, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
+          pstore(y + peelingEnd, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
        }
      }

@@ -381,8 +381,8 @@ struct apply_rotation_in_the_plane_selector<Scalar, OtherScalar, SizeAtCompileTi
      for (Index i = 0; i < size; i += PacketSize) {
        Packet xi = pload<Packet>(px);
        Packet yi = pload<Packet>(py);
-        pstore(px, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
-        pstore(py, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
+        pstore(px, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
+        pstore(py, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
        px += PacketSize;
        py += PacketSize;
      }
@@ -420,6 +420,34 @@ EIGEN_DEVICE_FUNC void inline apply_rotation_in_the_plane(DenseBase<VectorX>& xp
      x, incrx, y, incry, size, c, s);
 }

+template <typename MatrixType, typename RealScalar, typename Index>
+void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, JacobiRotation<RealScalar>* j_left,
+                         JacobiRotation<RealScalar>* j_right) {
+  using std::abs;
+  using std::sqrt;
+  Matrix<RealScalar, 2, 2> m;
+  m << numext::real(matrix.coeff(p, p)), numext::real(matrix.coeff(p, q)), numext::real(matrix.coeff(q, p)),
+      numext::real(matrix.coeff(q, q));
+  JacobiRotation<RealScalar> rot1;
+  RealScalar t = m.coeff(0, 0) + m.coeff(1, 1);
+  RealScalar d = m.coeff(1, 0) - m.coeff(0, 1);
+
+  if (abs(d) < (std::numeric_limits<RealScalar>::min)()) {
+    rot1.s() = RealScalar(0);
+    rot1.c() = RealScalar(1);
+  } else {
+    // If d!=0, then t/d cannot overflow because the magnitude of the
+    // entries forming d are not too small compared to the ones forming t.
+    RealScalar u = t / d;
+    RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u));
+    rot1.s() = RealScalar(1) / tmp;
+    rot1.c() = u / tmp;
+  }
+  m.applyOnTheLeft(0, 1, rot1);
+  j_right->makeJacobi(m, 0, 1);
+  *j_left = rot1 * j_right->transpose();
+}
+
 }  // end namespace internal

 }  // end namespace Eigen
--- a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
+++ b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
@@ -89,7 +89,7 @@ struct simpl_chol_helper {
        m_set[u] = v;
        u = next;
      }
-    };
+    }
  };

  // Computes the higher adjacency pattern by transposing the input lower adjacency matrix.
--- a/Eigen/src/misc/RealSvd2x2.h
+++ b/Eigen/src/misc/RealSvd2x2.h
@@ -1,53 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
-// Copyright (C) 2013-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_REALSVD2X2_H
-#define EIGEN_REALSVD2X2_H
-
-// IWYU pragma: private
-#include "./InternalHeaderCheck.h"
-
-namespace Eigen {
-
-namespace internal {
-
-template <typename MatrixType, typename RealScalar, typename Index>
-void real_2x2_jacobi_svd(const MatrixType &matrix, Index p, Index q, JacobiRotation<RealScalar> *j_left,
-                         JacobiRotation<RealScalar> *j_right) {
-  using std::abs;
-  using std::sqrt;
-  Matrix<RealScalar, 2, 2> m;
-  m << numext::real(matrix.coeff(p, p)), numext::real(matrix.coeff(p, q)), numext::real(matrix.coeff(q, p)),
-      numext::real(matrix.coeff(q, q));
-  JacobiRotation<RealScalar> rot1;
-  RealScalar t = m.coeff(0, 0) + m.coeff(1, 1);
-  RealScalar d = m.coeff(1, 0) - m.coeff(0, 1);
-
-  if (abs(d) < (std::numeric_limits<RealScalar>::min)()) {
-    rot1.s() = RealScalar(0);
-    rot1.c() = RealScalar(1);
-  } else {
-    // If d!=0, then t/d cannot overflow because the magnitude of the
-    // entries forming d are not too small compared to the ones forming t.
-    RealScalar u = t / d;
-    RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u));
-    rot1.s() = RealScalar(1) / tmp;
-    rot1.c() = u / tmp;
-  }
-  m.applyOnTheLeft(0, 1, rot1);
-  j_right->makeJacobi(m, 0, 1);
-  *j_left = rot1 * j_right->transpose();
-}
-
-}  // end namespace internal
-
-}  // end namespace Eigen
-
-#endif  // EIGEN_REALSVD2X2_H
--- a/ci/build.linux.gitlab-ci.yml
+++ b/ci/build.linux.gitlab-ci.yml
@@ -12,7 +12,6 @@
  rules:
    - if: $CI_PIPELINE_SOURCE == "schedule" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "web" && $CI_PROJECT_NAMESPACE == "libeigen"
-    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "merge_request_event" && $CI_PROJECT_NAMESPACE == "libeigen" && $CI_MERGE_REQUEST_LABELS =~ "/all-tests/"
  cache:
    key: "$CI_JOB_NAME_SLUG-$CI_COMMIT_REF_SLUG-BUILD"
@@ -110,7 +109,7 @@ build:linux:docs:
  rules:
    - if: $CI_PIPELINE_SOURCE == "schedule" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "web" && $CI_PROJECT_NAMESPACE == "libeigen"
-    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen"
+    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen" && $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
    - if: $CI_PIPELINE_SOURCE == "merge_request_event" && $CI_PROJECT_NAMESPACE == "libeigen" && $CI_MERGE_REQUEST_LABELS =~ "/all-tests/"


--- a/ci/build.windows.gitlab-ci.yml
+++ b/ci/build.windows.gitlab-ci.yml
@@ -17,7 +17,6 @@
  rules:
    - if: $CI_PIPELINE_SOURCE == "schedule" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "web" && $CI_PROJECT_NAMESPACE == "libeigen"
-    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "merge_request_event" && $CI_PROJECT_NAMESPACE == "libeigen" && $CI_MERGE_REQUEST_LABELS =~ "/all-tests/"

  cache:
--- a/ci/deploy.gitlab-ci.yml
+++ b/ci/deploy.gitlab-ci.yml
@@ -36,6 +36,6 @@ deploy:docs:
  rules:
    - if: $CI_PIPELINE_SOURCE == "schedule" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "web" && $CI_PROJECT_NAMESPACE == "libeigen"
-    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen"
+    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen" && $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
      variables:
        PAGES_PREFIX: docs-$CI_COMMIT_REF_NAME
--- a/ci/test.linux.gitlab-ci.yml
+++ b/ci/test.linux.gitlab-ci.yml
@@ -8,7 +8,6 @@
  rules:
    - if: $CI_PIPELINE_SOURCE == "schedule" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "web" && $CI_PROJECT_NAMESPACE == "libeigen"
-    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "merge_request_event" && $CI_PROJECT_NAMESPACE == "libeigen" && $CI_MERGE_REQUEST_LABELS =~ "/all-tests/"
  tags:
    - saas-linux-2xlarge-amd64
--- a/ci/test.windows.gitlab-ci.yml
+++ b/ci/test.windows.gitlab-ci.yml
@@ -8,7 +8,6 @@
  rules:
    - if: $CI_PIPELINE_SOURCE == "schedule" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "web" && $CI_PROJECT_NAMESPACE == "libeigen"
-    - if: $CI_PIPELINE_SOURCE == "push" && $CI_PROJECT_NAMESPACE == "libeigen"
    - if: $CI_PIPELINE_SOURCE == "merge_request_event" && $CI_PROJECT_NAMESPACE == "libeigen" && $CI_MERGE_REQUEST_LABELS =~ "/all-tests/"
  tags:
    - eigen-runner
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -257,6 +257,7 @@ ei_add_test(eigensolver_selfadjoint)
 ei_add_test(eigensolver_generic)
 ei_add_test(eigensolver_complex)
 ei_add_test(real_qz)
+ei_add_test(complex_qz)
 ei_add_test(eigensolver_generalized_real)
 ei_add_test(jacobi)
 ei_add_test(jacobisvd)
--- a/test/complex_qz.cpp
+++ b/test/complex_qz.cpp
@@ -0,0 +1,84 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 The Eigen Authors
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#define EIGEN_RUNTIME_NO_MALLOC
+#include "main.h"
+
+#include <Eigen/Eigenvalues>
+
+/* this test covers the following files:
+   ComplexQZ.h
+*/
+
+template <typename MatrixType>
+void generate_random_matrix_pair(const Index dim, MatrixType& A, MatrixType& B) {
+  A.setRandom(dim, dim);
+  B.setRandom(dim, dim);
+  // Zero out each row of B to with a probability of 10%.
+  for (int i = 0; i < dim; i++) {
+    if (internal::random<int>(0, 10) == 0) B.row(i).setZero();
+  }
+}
+
+template <typename MatrixType>
+void complex_qz(const MatrixType& A, const MatrixType& B) {
+  using std::abs;
+  const Index dim = A.rows();
+  ComplexQZ<MatrixType> qz(A, B);
+  VERIFY_IS_EQUAL(qz.info(), Success);
+  auto T = qz.matrixT(), S = qz.matrixS();
+  bool is_all_zero_T = true, is_all_zero_S = true;
+  using RealScalar = typename MatrixType::RealScalar;
+  RealScalar tol = dim * 10 * NumTraits<RealScalar>::epsilon();
+  for (Index j = 0; j < dim; j++) {
+    for (Index i = j + 1; i < dim; i++) {
+      if (std::abs(T(i, j)) > tol) {
+        std::cerr << std::abs(T(i, j)) << std::endl;
+        is_all_zero_T = false;
+      }
+      if (std::abs(S(i, j)) > tol) {
+        std::cerr << std::abs(S(i, j)) << std::endl;
+        is_all_zero_S = false;
+      }
+    }
+  }
+  VERIFY_IS_EQUAL(is_all_zero_T, true);
+  VERIFY_IS_EQUAL(is_all_zero_S, true);
+  VERIFY_IS_APPROX(qz.matrixQ() * qz.matrixS() * qz.matrixZ(), A);
+  VERIFY_IS_APPROX(qz.matrixQ() * qz.matrixT() * qz.matrixZ(), B);
+  VERIFY_IS_APPROX(qz.matrixQ() * qz.matrixQ().adjoint(), MatrixType::Identity(dim, dim));
+  VERIFY_IS_APPROX(qz.matrixZ() * qz.matrixZ().adjoint(), MatrixType::Identity(dim, dim));
+}
+
+EIGEN_DECLARE_TEST(complex_qz) {
+  for (int i = 0; i < g_repeat; i++) {
+    // Check for very small, fixed-sized double- and float complex matrices
+    Eigen::Matrix2cd A_2x2, B_2x2;
+    A_2x2.setRandom();
+    B_2x2.setRandom();
+    B_2x2.row(1).setZero();
+    Eigen::Matrix3cf A_3x3, B_3x3;
+    A_3x3.setRandom();
+    B_3x3.setRandom();
+    B_3x3.col(i % 3).setRandom();
+    CALL_SUBTEST_1(complex_qz(A_2x2, B_2x2));
+    CALL_SUBTEST_2(complex_qz(A_3x3, B_3x3));
+
+    // Test for float complex matrices
+    const Index dim = internal::random<Index>(15, 80);
+    Eigen::MatrixXcf A_float, B_float;
+    generate_random_matrix_pair(dim, A_float, B_float);
+    CALL_SUBTEST_3(complex_qz(A_float, B_float));
+
+    // Test for double complex matrices
+    Eigen::MatrixXcd A_double, B_double;
+    generate_random_matrix_pair(dim, A_double, B_double);
+    CALL_SUBTEST_4(complex_qz(A_double, B_double));
+  }
+}
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -13,6 +13,7 @@
 #include <limits>
 #include <Eigen/Eigenvalues>
 #include <Eigen/SparseCore>
+#include <unsupported/Eigen/MatrixFunctions>

 template <typename MatrixType>
 void selfadjointeigensolver_essential_check(const MatrixType& m) {
@@ -135,11 +136,13 @@ void selfadjointeigensolver(const MatrixType& m) {
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
+  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorExp());

  eiSymmUninitialized.compute(symmA, false);
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
+  VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorExp());

  // test Tridiagonalization's methods
  Tridiagonalization<MatrixType> tridiag(symmC);
@@ -167,6 +170,14 @@ void selfadjointeigensolver(const MatrixType& m) {
                                            eiSymmTridiag.eigenvectors().real().transpose());
  }

+  // Test matrix expponential from eigendecomposition.
+  // First scale to avoid overflow.
+  symmB = symmB / symmB.norm();
+  eiSymm.compute(symmB);
+  MatrixType expSymmB = eiSymm.operatorExp();
+  symmB = symmB.template selfadjointView<Lower>();
+  VERIFY_IS_APPROX(expSymmB, symmB.exp());
+
  if (rows > 1 && rows < 20) {
    // Test matrix with NaN
    symmC(0, 0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -747,22 +747,6 @@ void packetmath() {
    }
  }

-  if (PacketTraits::HasBlend) {
-    Packet thenPacket = internal::pload<Packet>(data1);
-    Packet elsePacket = internal::pload<Packet>(data2);
-    EIGEN_ALIGN_MAX internal::Selector<PacketSize> selector;
-    for (int i = 0; i < PacketSize; ++i) {
-      selector.select[i] = i;
-    }
-
-    Packet blend = internal::pblend(selector, thenPacket, elsePacket);
-    EIGEN_ALIGN_MAX Scalar result[size];
-    internal::pstore(result, blend);
-    for (int i = 0; i < PacketSize; ++i) {
-      VERIFY(test::isApproxAbs(result[i], (selector.select[i] ? data1[i] : data2[i]), refvalue));
-    }
-  }
-
  {
    for (int i = 0; i < PacketSize; ++i) {
      // "if" mask
--- a/unsupported/Eigen/CXX11/src/Tensor/TensorContractionMapper.h
+++ b/unsupported/Eigen/CXX11/src/Tensor/TensorContractionMapper.h
@@ -117,6 +117,10 @@ class SimpleTensorContractionMapper {
    return m_tensor.coeff(computeIndex(row, col));
  }

+#ifdef EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar operator[](Index row, Index col) const { return operator()(row, col); }
+#endif
+
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index computeIndex(Index row, Index col) const {
    const bool left = (side == Lhs);
    EIGEN_UNUSED_VARIABLE(left);  // annoying bug in g++8.1: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=85963
--- a/unsupported/Eigen/CXX11/src/Tensor/TensorEvaluator.h
+++ b/unsupported/Eigen/CXX11/src/Tensor/TensorEvaluator.h
@@ -691,9 +691,9 @@ struct TensorEvaluator<const TensorSelectOp<IfArgType, ThenArgType, ElseArgType>
  static constexpr int Layout = TensorEvaluator<IfArgType, Device>::Layout;
  enum {
    IsAligned = TensorEvaluator<ThenArgType, Device>::IsAligned & TensorEvaluator<ElseArgType, Device>::IsAligned,
-    PacketAccess = (TensorEvaluator<ThenArgType, Device>::PacketAccess &&
-                    TensorEvaluator<ElseArgType, Device>::PacketAccess && PacketType<Scalar, Device>::HasBlend) ||
-                   TernaryPacketAccess,
+    PacketAccess =
+        (TensorEvaluator<ThenArgType, Device>::PacketAccess && TensorEvaluator<ElseArgType, Device>::PacketAccess) ||
+        TernaryPacketAccess,
    BlockAccess = TensorEvaluator<IfArgType, Device>::BlockAccess &&
                  TensorEvaluator<ThenArgType, Device>::BlockAccess &&
                  TensorEvaluator<ElseArgType, Device>::BlockAccess,
@@ -789,13 +789,14 @@ struct TensorEvaluator<const TensorSelectOp<IfArgType, ThenArgType, ElseArgType>

  template <int LoadMode, bool UseTernary = TernaryPacketAccess, std::enable_if_t<!UseTernary, bool> = true>
  EIGEN_DEVICE_FUNC PacketReturnType packet(Index index) const {
-    internal::Selector<PacketSize> select;
+    EIGEN_ALIGN_TO_BOUNDARY(sizeof(PacketReturnType)) std::remove_const_t<Scalar> arr[PacketSize];
    EIGEN_UNROLL_LOOP
    for (Index i = 0; i < PacketSize; ++i) {
-      select.select[i] = m_condImpl.coeff(index + i);
+      arr[i] = m_condImpl.coeff(index + i) ? Scalar(-1) : Scalar(0);
    }
-    return internal::pblend(select, m_thenImpl.template packet<LoadMode>(index),
-                            m_elseImpl.template packet<LoadMode>(index));
+    return TernarySelectOp().template packetOp<PacketReturnType>(m_thenImpl.template packet<LoadMode>(index),
+                                                                 m_elseImpl.template packet<LoadMode>(index),
+                                                                 internal::pload<PacketReturnType>(arr));
  }

  template <int LoadMode, bool UseTernary = TernaryPacketAccess, std::enable_if_t<UseTernary, bool> = true>
--- a/unsupported/Eigen/src/Polynomials/Companion.h
+++ b/unsupported/Eigen/src/Polynomials/Companion.h
@@ -63,6 +63,10 @@ class companion {
    }
  }

+#ifdef EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+  EIGEN_STRONG_INLINE const Scalar_ operator[](Index row, Index col) const { return operator()(row, col); }
+#endif
+
 public:
  template <typename VectorType>
  void setPolynomial(const VectorType& poly) {
--- a/unsupported/Eigen/src/SparseExtra/RandomSetter.h
+++ b/unsupported/Eigen/src/SparseExtra/RandomSetter.h
@@ -280,6 +280,11 @@ class RandomSetter {
    return m_hashmaps[outerMajor][key].value;
  }

+#ifdef EIGEN_MULTIDIMENSIONAL_SUBSCRIPT
+  /** \returns a reference to the coefficient at given coordinates \a row, \a col */
+  Scalar& operator[](Index row, Index col) { return operator()(row, col); }
+#endif
+
  /** \returns the number of non zero coefficients
   *
   * \note According to the underlying map/hash_map implementation,