diff --git a/Eigen/src/AccelerateSupport/AccelerateSupport.h b/Eigen/src/AccelerateSupport/AccelerateSupport.h
index 13a26dfbb..c944aeabd 100644
--- a/Eigen/src/AccelerateSupport/AccelerateSupport.h
+++ b/Eigen/src/AccelerateSupport/AccelerateSupport.h
@@ -110,7 +110,7 @@ using AccelerateCholeskyAtA = AccelerateImpl<MatrixType, 0, SparseFactorizationC
 namespace internal {
 template <typename T>
 struct AccelFactorizationDeleter {
-  void operator()(T* sym) {
+  void operator()(T* sym) const {
     if (sym) {
       SparseCleanup(*sym);
       delete sym;
diff --git a/Eigen/src/Core/arch/AVX512/GemmKernel.h b/Eigen/src/Core/arch/AVX512/GemmKernel.h
index fbc2d5822..556894c14 100644
--- a/Eigen/src/Core/arch/AVX512/GemmKernel.h
+++ b/Eigen/src/Core/arch/AVX512/GemmKernel.h
@@ -977,12 +977,12 @@ struct gemm_pack_rhs<Scalar, Index, DataMapper, 8, ColMajor, Conjugate, PanelMod
   typedef typename DataMapper::LinearMapper LinearMapper;
   enum { PacketSize = packet_traits<Scalar>::size };
   EIGEN_DONT_INLINE void operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                                    Index offset = 0);
+                                    Index offset = 0) const;
 };
 
 template <typename Scalar, typename Index, typename DataMapper, bool Conjugate, bool PanelMode>
 EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, DataMapper, 8, ColMajor, Conjugate, PanelMode>::operator()(
-    Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   constexpr int nr = 8;
   EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS COLMAJOR");
   EIGEN_UNUSED_VARIABLE(stride);
@@ -1111,7 +1111,7 @@ struct gemm_pack_rhs<Scalar, Index, DataMapper, 8, RowMajor, Conjugate, PanelMod
     QuarterPacketSize = unpacket_traits<QuarterPacket>::size
   };
   EIGEN_DONT_INLINE void operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                                    Index offset = 0) {
+                                    Index offset = 0) const {
     constexpr int nr = 8;
     EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS ROWMAJOR");
     EIGEN_UNUSED_VARIABLE(stride);
@@ -1211,13 +1211,13 @@ template <typename Scalar, typename Index, typename DataMapper, int mr, bool Con
 struct gebp_kernel<Scalar, Scalar, Index, DataMapper, mr, 8, ConjugateLhs, ConjugateRhs> {
   EIGEN_ALWAYS_INLINE void operator()(const DataMapper& res, const Scalar* blockA, const Scalar* blockB, Index rows,
                                       Index depth, Index cols, Scalar alpha, Index strideA = -1, Index strideB = -1,
-                                      Index offsetA = 0, Index offsetB = 0);
+                                      Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Scalar, typename Index, typename DataMapper, int mr, bool ConjugateLhs, bool ConjugateRhs>
 EIGEN_ALWAYS_INLINE void gebp_kernel<Scalar, Scalar, Index, DataMapper, mr, 8, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const Scalar* blockA, const Scalar* blockB, Index rows, Index depth, Index cols,
-    Scalar alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    Scalar alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   if (res.incr() == 1) {
     if (alpha == 1) {
       gemm_kern_avx512<Scalar, mr, 8, true, false, true>(rows, cols, depth, &alpha, blockA, blockB, (Scalar*)res.data(),
diff --git a/Eigen/src/Core/arch/AltiVec/MatrixProduct.h b/Eigen/src/Core/arch/AltiVec/MatrixProduct.h
index 472552610..e31213ee7 100644
--- a/Eigen/src/Core/arch/AltiVec/MatrixProduct.h
+++ b/Eigen/src/Core/arch/AltiVec/MatrixProduct.h
@@ -290,7 +290,7 @@ EIGEN_STRONG_INLINE void symm_pack_lhs_helper(Scalar* blockA, const Scalar* _lhs
 template <typename Index, int nr, int StorageOrder>
 struct symm_pack_rhs<std::complex<float>, Index, nr, StorageOrder> {
   void operator()(std::complex<float>* blockB, const std::complex<float>* _rhs, Index rhsStride, Index rows, Index cols,
-                  Index k2) {
+                  Index k2) const {
     symm_pack_complex_rhs_helper<float, StorageOrder, 1>(blockB, _rhs, rhsStride, rows, cols, k2);
   }
 };
@@ -298,7 +298,7 @@ struct symm_pack_rhs<std::complex<float>, Index, nr, StorageOrder> {
 template <typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
 struct symm_pack_lhs<std::complex<float>, Index, Pack1, Pack2_dummy, StorageOrder> {
   void operator()(std::complex<float>* blockA, const std::complex<float>* _lhs, Index lhsStride, Index cols,
-                  Index rows) {
+                  Index rows) const {
     symm_pack_complex_lhs_helper<float, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
   }
 };
@@ -308,7 +308,7 @@ struct symm_pack_lhs<std::complex<float>, Index, Pack1, Pack2_dummy, StorageOrde
 template <typename Index, int nr, int StorageOrder>
 struct symm_pack_rhs<std::complex<double>, Index, nr, StorageOrder> {
   void operator()(std::complex<double>* blockB, const std::complex<double>* _rhs, Index rhsStride, Index rows,
-                  Index cols, Index k2) {
+                  Index cols, Index k2) const {
     symm_pack_complex_rhs_helper<double, StorageOrder, 2>(blockB, _rhs, rhsStride, rows, cols, k2);
   }
 };
@@ -316,7 +316,7 @@ struct symm_pack_rhs<std::complex<double>, Index, nr, StorageOrder> {
 template <typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
 struct symm_pack_lhs<std::complex<double>, Index, Pack1, Pack2_dummy, StorageOrder> {
   void operator()(std::complex<double>* blockA, const std::complex<double>* _lhs, Index lhsStride, Index cols,
-                  Index rows) {
+                  Index rows) const {
     symm_pack_complex_lhs_helper<double, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
   }
 };
@@ -324,14 +324,14 @@ struct symm_pack_lhs<std::complex<double>, Index, Pack1, Pack2_dummy, StorageOrd
 // *********** symm_pack float32 ***********
 template <typename Index, int nr, int StorageOrder>
 struct symm_pack_rhs<float, Index, nr, StorageOrder> {
-  void operator()(float* blockB, const float* _rhs, Index rhsStride, Index rows, Index cols, Index k2) {
+  void operator()(float* blockB, const float* _rhs, Index rhsStride, Index rows, Index cols, Index k2) const {
     symm_pack_rhs_helper<float, StorageOrder, 1>(blockB, _rhs, rhsStride, rows, cols, k2);
   }
 };
 
 template <typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
 struct symm_pack_lhs<float, Index, Pack1, Pack2_dummy, StorageOrder> {
-  void operator()(float* blockA, const float* _lhs, Index lhsStride, Index cols, Index rows) {
+  void operator()(float* blockA, const float* _lhs, Index lhsStride, Index cols, Index rows) const {
     symm_pack_lhs_helper<float, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
   }
 };
@@ -339,14 +339,14 @@ struct symm_pack_lhs<float, Index, Pack1, Pack2_dummy, StorageOrder> {
 // *********** symm_pack float64 ***********
 template <typename Index, int nr, int StorageOrder>
 struct symm_pack_rhs<double, Index, nr, StorageOrder> {
-  void operator()(double* blockB, const double* _rhs, Index rhsStride, Index rows, Index cols, Index k2) {
+  void operator()(double* blockB, const double* _rhs, Index rhsStride, Index rows, Index cols, Index k2) const {
     symm_pack_rhs_helper<double, StorageOrder, 2>(blockB, _rhs, rhsStride, rows, cols, k2);
   }
 };
 
 template <typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
 struct symm_pack_lhs<double, Index, Pack1, Pack2_dummy, StorageOrder> {
-  void operator()(double* blockA, const double* _lhs, Index lhsStride, Index cols, Index rows) {
+  void operator()(double* blockA, const double* _lhs, Index lhsStride, Index cols, Index rows) const {
     symm_pack_lhs_helper<double, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
   }
 };
@@ -381,7 +381,7 @@ template <typename Scalar, typename DataMapper, typename Packet, typename Packet
 struct dhs_cpack {
   template <bool transpose>
   EIGEN_ALWAYS_INLINE void dhs_cblock(PacketBlock<PacketC, 8>& cblock, PacketBlock<Packet, 4>& block,
-                                      Packet16uc permute) {
+                                      Packet16uc permute) const {
     if (transpose) {
       block.packet[0] = vec_perm(cblock.packet[0].v, cblock.packet[1].v, permute);
       block.packet[1] = vec_perm(cblock.packet[2].v, cblock.packet[3].v, permute);
@@ -453,7 +453,7 @@ struct dhs_cpack {
   }
 
   EIGEN_STRONG_INLINE void operator()(std::complex<Scalar>* blockA, const DataMapper& lhs, Index depth, Index rows,
-                                      Index stride, Index offset) {
+                                      Index stride, Index offset) const {
     const Index vectorSize = quad_traits<Scalar>::vectorsize;
     const Index vectorDelta = vectorSize * ((PanelMode) ? stride : depth);
     Index rir = ((PanelMode) ? (vectorSize * offset) : 0), rii;
@@ -557,7 +557,7 @@ template <typename Scalar, typename DataMapper, typename Packet, int StorageOrde
 struct dhs_pack {
   template <Index n>
   EIGEN_ALWAYS_INLINE void dhs_copy(Scalar* blockA, const DataMapper& lhs2, Index& i, Index& ri, Index depth,
-                                    const Index vectorSize) {
+                                    const Index vectorSize) const {
     PacketBlock<Packet, 4> block[n];
 
     for (; i + n * vectorSize <= depth; i += n * vectorSize) {
@@ -584,7 +584,7 @@ struct dhs_pack {
   }
 
   EIGEN_STRONG_INLINE void operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride,
-                                      Index offset) {
+                                      Index offset) const {
     const Index vectorSize = quad_traits<Scalar>::vectorsize;
     Index ri = 0, j = 0;
 
@@ -660,7 +660,7 @@ template <typename DataMapper, int StorageOrder, bool PanelMode>
 struct dhs_pack<double, DataMapper, Packet2d, StorageOrder, PanelMode, true> {
   template <Index n>
   EIGEN_ALWAYS_INLINE void dhs_copy(double* blockA, const DataMapper& lhs2, Index& i, Index& ri, Index depth,
-                                    const Index vectorSize) {
+                                    const Index vectorSize) const {
     PacketBlock<Packet2d, 2> block[n];
 
     for (; i + n * vectorSize <= depth; i += n * vectorSize) {
@@ -689,7 +689,7 @@ struct dhs_pack<double, DataMapper, Packet2d, StorageOrder, PanelMode, true> {
   }
 
   EIGEN_STRONG_INLINE void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride,
-                                      Index offset) {
+                                      Index offset) const {
     const Index vectorSize = quad_traits<double>::vectorsize;
     Index ri = 0, j = 0;
 
@@ -737,7 +737,7 @@ template <typename DataMapper, int StorageOrder, bool PanelMode>
 struct dhs_pack<double, DataMapper, Packet2d, StorageOrder, PanelMode, false> {
   template <Index n>
   EIGEN_ALWAYS_INLINE void dhs_copy(double* blockB, const DataMapper& rhs2, Index& i, Index& ri, Index depth,
-                                    const Index vectorSize) {
+                                    const Index vectorSize) const {
     PacketBlock<Packet2d, 2> block1[n], block2[n];
     PacketBlock<Packet2d, 4> block3[n];
 
@@ -779,7 +779,7 @@ struct dhs_pack<double, DataMapper, Packet2d, StorageOrder, PanelMode, false> {
   }
 
   EIGEN_STRONG_INLINE void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride,
-                                      Index offset) {
+                                      Index offset) const {
     const Index vectorSize = quad_traits<double>::vectorsize;
     Index ri = 0, j = 0;
 
@@ -835,7 +835,7 @@ struct dhs_pack<double, DataMapper, Packet2d, StorageOrder, PanelMode, false> {
 template <typename DataMapper, int StorageOrder, bool PanelMode>
 struct dhs_pack<bfloat16, DataMapper, Packet8bf, StorageOrder, PanelMode, true> {
   EIGEN_STRONG_INLINE void operator()(bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride,
-                                      Index offset) {
+                                      Index offset) const {
     const Index vectorSize = quad_traits<bfloat16>::vectorsize;
     Index ri = 0, j = 0;
 
@@ -1178,7 +1178,7 @@ struct dhs_pack<bfloat16, DataMapper, Packet8bf, StorageOrder, PanelMode, true>
 template <typename DataMapper, int StorageOrder, bool PanelMode>
 struct dhs_pack<bfloat16, DataMapper, Packet8bf, StorageOrder, PanelMode, false> {
   EIGEN_STRONG_INLINE void operator()(bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride,
-                                      Index offset) {
+                                      Index offset) const {
     const Index vectorSize = quad_traits<bfloat16>::vectorsize;
     Index ri = 0, j = 0;
 
@@ -1350,7 +1350,7 @@ struct dhs_cpack<double, DataMapper, Packet, PacketC, StorageOrder, Conjugate, P
   }
 
   EIGEN_STRONG_INLINE void operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows,
-                                      Index stride, Index offset) {
+                                      Index stride, Index offset) const {
     const Index vectorSize = quad_traits<double>::vectorsize;
     const Index vectorDelta = vectorSize * ((PanelMode) ? stride : depth);
     Index rir = ((PanelMode) ? (vectorSize * offset) : 0), rii;
@@ -1442,7 +1442,7 @@ struct dhs_cpack<double, DataMapper, Packet, PacketC, StorageOrder, Conjugate, P
   }
 
   EIGEN_STRONG_INLINE void operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols,
-                                      Index stride, Index offset) {
+                                      Index stride, Index offset) const {
     const Index vectorSize = quad_traits<double>::vectorsize;
     const Index vectorDelta = 2 * vectorSize * ((PanelMode) ? stride : depth);
     Index rir = ((PanelMode) ? (2 * vectorSize * offset) : 0), rii;
@@ -3162,24 +3162,26 @@ void gemmbfloat16(const DataMapper& res, const bfloat16* indexA, const bfloat16*
  * **********************************/
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode> {
-  void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0, Index offset = 0);
+  void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>::operator()(
-    double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) {
+    double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) const {
   dhs_pack<double, DataMapper, Packet2d, ColMajor, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode> {
-  void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0, Index offset = 0);
+  void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>::operator()(
-    double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) {
+    double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) const {
   dhs_pack<double, DataMapper, Packet2d, RowMajor, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
@@ -3187,48 +3189,52 @@ void gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Co
 #if EIGEN_ALTIVEC_USE_CUSTOM_PACK
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<double, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode> {
-  void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0, Index offset = 0);
+  void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<double, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>::operator()(
-    double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_pack<double, DataMapper, Packet2d, ColMajor, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<double, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode> {
-  void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0, Index offset = 0);
+  void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<double, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>::operator()(
-    double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_pack<double, DataMapper, Packet2d, RowMajor, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<bfloat16, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode> {
-  void operator()(bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0, Index offset = 0);
+  void operator()(bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<bfloat16, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>::operator()(
-    bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_pack<bfloat16, DataMapper, Packet8bf, ColMajor, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<bfloat16, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode> {
-  void operator()(bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0, Index offset = 0);
+  void operator()(bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<bfloat16, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>::operator()(
-    bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    bfloat16* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_pack<bfloat16, DataMapper, Packet8bf, RowMajor, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
@@ -3236,48 +3242,52 @@ void gemm_pack_rhs<bfloat16, Index, DataMapper, nr, RowMajor, Conjugate, PanelMo
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<bfloat16, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode> {
-  void operator()(bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0, Index offset = 0);
+  void operator()(bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<bfloat16, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>::operator()(
-    bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) {
+    bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) const {
   dhs_pack<bfloat16, DataMapper, Packet8bf, ColMajor, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<bfloat16, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode> {
-  void operator()(bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0, Index offset = 0);
+  void operator()(bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<bfloat16, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>::operator()(
-    bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) {
+    bfloat16* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) const {
   dhs_pack<bfloat16, DataMapper, Packet8bf, RowMajor, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode> {
-  void operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0, Index offset = 0);
+  void operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>::operator()(
-    float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) {
+    float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) const {
   dhs_pack<float, DataMapper, Packet4f, RowMajor, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode> {
-  void operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0, Index offset = 0);
+  void operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>::operator()(
-    float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) {
+    float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset) const {
   dhs_pack<float, DataMapper, Packet4f, ColMajor, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
@@ -3285,13 +3295,13 @@ void gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Con
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode> {
   void operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate,
                    PanelMode>::operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows,
-                                          Index stride, Index offset) {
+                                          Index stride, Index offset) const {
   dhs_cpack<float, DataMapper, Packet4f, Packet2cf, RowMajor, Conjugate, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
@@ -3299,13 +3309,13 @@ void gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet,
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode> {
   void operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate,
                    PanelMode>::operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows,
-                                          Index stride, Index offset) {
+                                          Index stride, Index offset) const {
   dhs_cpack<float, DataMapper, Packet4f, Packet2cf, ColMajor, Conjugate, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
@@ -3313,24 +3323,26 @@ void gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet,
 #if EIGEN_ALTIVEC_USE_CUSTOM_PACK
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<float, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode> {
-  void operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0, Index offset = 0);
+  void operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<float, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>::operator()(
-    float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_pack<float, DataMapper, Packet4f, ColMajor, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<float, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode> {
-  void operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0, Index offset = 0);
+  void operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<float, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>::operator()(
-    float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_pack<float, DataMapper, Packet4f, RowMajor, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
@@ -3339,12 +3351,12 @@ void gemm_pack_rhs<float, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode> {
   void operator()(std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>::operator()(
-    std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_cpack<float, DataMapper, Packet4f, Packet2cf, ColMajor, Conjugate, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
@@ -3352,12 +3364,12 @@ void gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, ColMajor, Conjuga
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode> {
   void operator()(std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>::operator()(
-    std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_cpack<float, DataMapper, Packet4f, Packet2cf, RowMajor, Conjugate, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
@@ -3365,13 +3377,13 @@ void gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, RowMajor, Conjuga
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode> {
   void operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate,
                    PanelMode>::operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows,
-                                          Index stride, Index offset) {
+                                          Index stride, Index offset) const {
   dhs_cpack<double, DataMapper, Packet2d, Packet1cd, RowMajor, Conjugate, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
@@ -3379,13 +3391,13 @@ void gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 struct gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode> {
   void operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
 void gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate,
                    PanelMode>::operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows,
-                                          Index stride, Index offset) {
+                                          Index stride, Index offset) const {
   dhs_cpack<double, DataMapper, Packet2d, Packet1cd, ColMajor, Conjugate, PanelMode, true> pack;
   pack(blockA, lhs, depth, rows, stride, offset);
 }
@@ -3393,12 +3405,12 @@ void gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode> {
   void operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>::operator()(
-    std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_cpack<double, DataMapper, Packet2d, Packet1cd, ColMajor, Conjugate, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
@@ -3406,12 +3418,12 @@ void gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, ColMajor, Conjug
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 struct gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode> {
   void operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                  Index offset = 0);
+                  Index offset = 0) const;
 };
 
 template <typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 void gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>::operator()(
-    std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   dhs_cpack<double, DataMapper, Packet2d, Packet1cd, RowMajor, Conjugate, PanelMode, false> pack;
   pack(blockB, rhs, depth, cols, stride, offset);
 }
@@ -3423,13 +3435,13 @@ struct gebp_kernel<float, float, Index, DataMapper, mr, nr, ConjugateLhs, Conjug
   typedef typename quad_traits<float>::rhstype RhsPacket;
 
   void operator()(const DataMapper& res, const float* blockA, const float* blockB, Index rows, Index depth, Index cols,
-                  float alpha, Index strideA = -1, Index strideB = -1, Index offsetA = 0, Index offsetB = 0);
+                  float alpha, Index strideA = -1, Index strideB = -1, Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
 void gebp_kernel<float, float, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const float* blockA, const float* blockB, Index rows, Index depth, Index cols, float alpha,
-    Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   const Index accRows = quad_traits<float>::rows;
   const Index accCols = quad_traits<float>::size;
   static void (*gemm_function)(const DataMapper&, const float*, const float*, Index, Index, Index, float, Index, Index,
@@ -3449,7 +3461,7 @@ struct gebp_kernel<std::complex<float>, std::complex<float>, Index, DataMapper,
 
   void operator()(const DataMapper& res, const std::complex<float>* blockA, const std::complex<float>* blockB,
                   Index rows, Index depth, Index cols, std::complex<float> alpha, Index strideA = -1,
-                  Index strideB = -1, Index offsetA = 0, Index offsetB = 0);
+                  Index strideB = -1, Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
@@ -3457,7 +3469,7 @@ void gebp_kernel<std::complex<float>, std::complex<float>, Index, DataMapper, mr
                  ConjugateRhs>::operator()(const DataMapper& res, const std::complex<float>* blockA,
                                            const std::complex<float>* blockB, Index rows, Index depth, Index cols,
                                            std::complex<float> alpha, Index strideA, Index strideB, Index offsetA,
-                                           Index offsetB) {
+                                           Index offsetB) const {
   const Index accRows = quad_traits<float>::rows;
   const Index accCols = quad_traits<float>::size;
   static void (*gemm_function)(const DataMapper&, const std::complex<float>*, const std::complex<float>*, Index, Index,
@@ -3482,13 +3494,13 @@ struct gebp_kernel<float, std::complex<float>, Index, DataMapper, mr, nr, Conjug
 
   void operator()(const DataMapper& res, const float* blockA, const std::complex<float>* blockB, Index rows,
                   Index depth, Index cols, std::complex<float> alpha, Index strideA = -1, Index strideB = -1,
-                  Index offsetA = 0, Index offsetB = 0);
+                  Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
 void gebp_kernel<float, std::complex<float>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const float* blockA, const std::complex<float>* blockB, Index rows, Index depth, Index cols,
-    std::complex<float> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    std::complex<float> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   const Index accRows = quad_traits<float>::rows;
   const Index accCols = quad_traits<float>::size;
   static void (*gemm_function)(const DataMapper&, const float*, const std::complex<float>*, Index, Index, Index,
@@ -3513,13 +3525,13 @@ struct gebp_kernel<std::complex<float>, float, Index, DataMapper, mr, nr, Conjug
 
   void operator()(const DataMapper& res, const std::complex<float>* blockA, const float* blockB, Index rows,
                   Index depth, Index cols, std::complex<float> alpha, Index strideA = -1, Index strideB = -1,
-                  Index offsetA = 0, Index offsetB = 0);
+                  Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
 void gebp_kernel<std::complex<float>, float, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const std::complex<float>* blockA, const float* blockB, Index rows, Index depth, Index cols,
-    std::complex<float> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    std::complex<float> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   const Index accRows = quad_traits<float>::rows;
   const Index accCols = quad_traits<float>::size;
   static void (*gemm_function)(const DataMapper&, const std::complex<float>*, const float*, Index, Index, Index,
@@ -3543,13 +3555,13 @@ struct gebp_kernel<double, double, Index, DataMapper, mr, nr, ConjugateLhs, Conj
 
   void operator()(const DataMapper& res, const double* blockA, const double* blockB, Index rows, Index depth,
                   Index cols, double alpha, Index strideA = -1, Index strideB = -1, Index offsetA = 0,
-                  Index offsetB = 0);
+                  Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
 void gebp_kernel<double, double, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const double* blockA, const double* blockB, Index rows, Index depth, Index cols,
-    double alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    double alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   const Index accRows = quad_traits<double>::rows;
   const Index accCols = quad_traits<double>::size;
   static void (*gemm_function)(const DataMapper&, const double*, const double*, Index, Index, Index, double, Index,
@@ -3569,7 +3581,7 @@ struct gebp_kernel<std::complex<double>, std::complex<double>, Index, DataMapper
 
   void operator()(const DataMapper& res, const std::complex<double>* blockA, const std::complex<double>* blockB,
                   Index rows, Index depth, Index cols, std::complex<double> alpha, Index strideA = -1,
-                  Index strideB = -1, Index offsetA = 0, Index offsetB = 0);
+                  Index strideB = -1, Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
@@ -3577,7 +3589,7 @@ void gebp_kernel<std::complex<double>, std::complex<double>, Index, DataMapper,
                  ConjugateRhs>::operator()(const DataMapper& res, const std::complex<double>* blockA,
                                            const std::complex<double>* blockB, Index rows, Index depth, Index cols,
                                            std::complex<double> alpha, Index strideA, Index strideB, Index offsetA,
-                                           Index offsetB) {
+                                           Index offsetB) const {
   const Index accRows = quad_traits<double>::rows;
   const Index accCols = quad_traits<double>::size;
   static void (*gemm_function)(const DataMapper&, const std::complex<double>*, const std::complex<double>*, Index,
@@ -3603,13 +3615,13 @@ struct gebp_kernel<std::complex<double>, double, Index, DataMapper, mr, nr, Conj
 
   void operator()(const DataMapper& res, const std::complex<double>* blockA, const double* blockB, Index rows,
                   Index depth, Index cols, std::complex<double> alpha, Index strideA = -1, Index strideB = -1,
-                  Index offsetA = 0, Index offsetB = 0);
+                  Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
 void gebp_kernel<std::complex<double>, double, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const std::complex<double>* blockA, const double* blockB, Index rows, Index depth,
-    Index cols, std::complex<double> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    Index cols, std::complex<double> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   const Index accRows = quad_traits<double>::rows;
   const Index accCols = quad_traits<double>::size;
   static void (*gemm_function)(const DataMapper&, const std::complex<double>*, const double*, Index, Index, Index,
@@ -3634,13 +3646,13 @@ struct gebp_kernel<double, std::complex<double>, Index, DataMapper, mr, nr, Conj
 
   void operator()(const DataMapper& res, const double* blockA, const std::complex<double>* blockB, Index rows,
                   Index depth, Index cols, std::complex<double> alpha, Index strideA = -1, Index strideB = -1,
-                  Index offsetA = 0, Index offsetB = 0);
+                  Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
 void gebp_kernel<double, std::complex<double>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const double* blockA, const std::complex<double>* blockB, Index rows, Index depth,
-    Index cols, std::complex<double> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    Index cols, std::complex<double> alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   const Index accRows = quad_traits<double>::rows;
   const Index accCols = quad_traits<double>::size;
   static void (*gemm_function)(const DataMapper&, const double*, const std::complex<double>*, Index, Index, Index,
@@ -3664,13 +3676,13 @@ struct gebp_kernel<bfloat16, bfloat16, Index, DataMapper, mr, nr, ConjugateLhs,
 
   void operator()(const DataMapper& res, const bfloat16* blockA, const bfloat16* blockB, Index rows, Index depth,
                   Index cols, bfloat16 alpha, Index strideA = -1, Index strideB = -1, Index offsetA = 0,
-                  Index offsetB = 0);
+                  Index offsetB = 0) const;
 };
 
 template <typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
 void gebp_kernel<bfloat16, bfloat16, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(
     const DataMapper& res, const bfloat16* blockA, const bfloat16* blockB, Index rows, Index depth, Index cols,
-    bfloat16 alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) {
+    bfloat16 alpha, Index strideA, Index strideB, Index offsetA, Index offsetB) const {
   static void (*gemm_function)(const DataMapper&, const bfloat16*, const bfloat16*, Index, Index, Index, bfloat16,
                                Index, Index, Index, Index) =
 #ifdef EIGEN_MATRIX_PRODUCT_MMA_ALTIVEC_H
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathPow.h b/Eigen/src/Core/arch/Default/GenericPacketMathPow.h
index 08686f86b..c24b8f381 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathPow.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathPow.h
@@ -164,7 +164,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcbrt_double(const Pa
 template <typename Scalar>
 struct accurate_log2 {
   template <typename Packet>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) const {
     log2_x_hi = plog2(x);
     log2_x_lo = pzero(x);
   }
@@ -185,7 +185,7 @@ struct accurate_log2 {
 template <>
 struct accurate_log2<float> {
   template <typename Packet>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& z, Packet& log2_x_hi, Packet& log2_x_lo) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& z, Packet& log2_x_hi, Packet& log2_x_lo) const {
     // Split the two lowest order constant coefficient into double-word representation.
     constexpr double kC0 = 1.442695041742110273474963832995854318141937255859375e+00;
     constexpr float kC0_hi = static_cast<float>(kC0);
@@ -231,7 +231,7 @@ struct accurate_log2<float> {
 template <>
 struct accurate_log2<double> {
   template <typename Packet>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Packet& x, Packet& log2_x_hi, Packet& log2_x_lo) const {
     // We use a transformation of variables:
     //    r = c * (x-1) / (x+1),
     // such that
diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h
index 15a3bbf0f..8adfdd6f9 100644
--- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h
+++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h
@@ -459,7 +459,7 @@ class gebp_traits {
   typedef QuadPacket<RhsPacket> RhsPacketx4;
   typedef ResPacket AccPacket;
 
-  EIGEN_STRONG_INLINE void initAcc(AccPacket& p) { p = pset1<ResPacket>(ResScalar(0)); }
+  EIGEN_STRONG_INLINE void initAcc(AccPacket& p) const { p = pset1<ResPacket>(ResScalar(0)); }
 
   template <typename RhsPacketType>
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketType& dest) const {
@@ -564,7 +564,7 @@ class gebp_traits<std::complex<RealScalar>, RealScalar, ConjLhs_, false, Arch, P
 
   typedef ResPacket AccPacket;
 
-  EIGEN_STRONG_INLINE void initAcc(AccPacket& p) { p = pset1<ResPacket>(ResScalar(0)); }
+  EIGEN_STRONG_INLINE void initAcc(AccPacket& p) const { p = pset1<ResPacket>(ResScalar(0)); }
 
   template <typename RhsPacketType>
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketType& dest) const {
@@ -752,9 +752,9 @@ class gebp_traits<std::complex<RealScalar>, std::complex<RealScalar>, ConjLhs_,
   // this actually holds 8 packets!
   typedef QuadPacket<RhsPacket> RhsPacketx4;
 
-  EIGEN_STRONG_INLINE void initAcc(Scalar& p) { p = Scalar(0); }
+  EIGEN_STRONG_INLINE void initAcc(Scalar& p) const { p = Scalar(0); }
 
-  EIGEN_STRONG_INLINE void initAcc(DoublePacketType& p) {
+  EIGEN_STRONG_INLINE void initAcc(DoublePacketType& p) const {
     p.first = pset1<RealPacket>(RealScalar(0));
     p.second = pset1<RealPacket>(RealScalar(0));
   }
@@ -896,7 +896,7 @@ class gebp_traits<RealScalar, std::complex<RealScalar>, false, ConjRhs_, Arch, P
   typedef QuadPacket<RhsPacket> RhsPacketx4;
   typedef ResPacket AccPacket;
 
-  EIGEN_STRONG_INLINE void initAcc(AccPacket& p) { p = pset1<ResPacket>(ResScalar(0)); }
+  EIGEN_STRONG_INLINE void initAcc(AccPacket& p) const { p = pset1<ResPacket>(ResScalar(0)); }
 
   template <typename RhsPacketType>
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketType& dest) const {
@@ -1021,7 +1021,7 @@ struct gebp_kernel {
 
   EIGEN_DONT_INLINE void operator()(const DataMapper& res, const LhsScalar* blockA, const RhsScalar* blockB, Index rows,
                                     Index depth, Index cols, ResScalar alpha, Index strideA = -1, Index strideB = -1,
-                                    Index offsetA = 0, Index offsetB = 0);
+                                    Index offsetA = 0, Index offsetB = 0) const;
 };
 
 template <typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr,
@@ -1040,7 +1040,7 @@ struct last_row_process_16_packets {
 
   EIGEN_STRONG_INLINE void operator()(const DataMapper& res, SwappedTraits& straits, const LhsScalar* blA,
                                       const RhsScalar* blB, Index depth, const Index endk, Index i, Index j2,
-                                      ResScalar alpha, SAccPacket& C0) {
+                                      ResScalar alpha, SAccPacket& C0) const {
     EIGEN_UNUSED_VARIABLE(res);
     EIGEN_UNUSED_VARIABLE(straits);
     EIGEN_UNUSED_VARIABLE(blA);
@@ -1068,7 +1068,7 @@ struct last_row_process_16_packets<LhsScalar, RhsScalar, Index, DataMapper, mr,
 
   EIGEN_STRONG_INLINE void operator()(const DataMapper& res, SwappedTraits& straits, const LhsScalar* blA,
                                       const RhsScalar* blB, Index depth, const Index endk, Index i, Index j2,
-                                      ResScalar alpha, SAccPacket& C0) {
+                                      ResScalar alpha, SAccPacket& C0) const {
     typedef typename unpacket_traits<typename unpacket_traits<SResPacket>::half>::half SResPacketQuarter;
     typedef typename unpacket_traits<typename unpacket_traits<SLhsPacket>::half>::half SLhsPacketQuarter;
     typedef typename unpacket_traits<typename unpacket_traits<SRhsPacket>::half>::half SRhsPacketQuarter;
@@ -1373,7 +1373,7 @@ EIGEN_DONT_INLINE void gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr,
                                    ConjugateRhs>::operator()(const DataMapper& res, const LhsScalar* blockA,
                                                              const RhsScalar* blockB, Index rows, Index depth,
                                                              Index cols, ResScalar alpha, Index strideA, Index strideB,
-                                                             Index offsetA, Index offsetB) {
+                                                             Index offsetA, Index offsetB) const {
   Traits traits;
   SwappedTraits straits;
 
@@ -1781,14 +1781,14 @@ template <typename Scalar, typename Index, typename DataMapper, int Pack1, int P
 struct gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode> {
   typedef typename DataMapper::LinearMapper LinearMapper;
   EIGEN_DONT_INLINE void operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
-                                    Index offset = 0);
+                                    Index offset = 0) const;
 };
 
 template <typename Scalar, typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate,
           bool PanelMode>
 EIGEN_DONT_INLINE void gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate,
                                      PanelMode>::operator()(Scalar* blockA, const DataMapper& lhs, Index depth,
-                                                            Index rows, Index stride, Index offset) {
+                                                            Index rows, Index stride, Index offset) const {
   typedef typename unpacket_traits<Packet>::half HalfPacket;
   typedef typename unpacket_traits<typename unpacket_traits<Packet>::half>::half QuarterPacket;
   enum {
@@ -1931,14 +1931,14 @@ template <typename Scalar, typename Index, typename DataMapper, int Pack1, int P
 struct gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode> {
   typedef typename DataMapper::LinearMapper LinearMapper;
   EIGEN_DONT_INLINE void operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride = 0,
-                                    Index offset = 0);
+                                    Index offset = 0) const;
 };
 
 template <typename Scalar, typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate,
           bool PanelMode>
 EIGEN_DONT_INLINE void gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate,
                                      PanelMode>::operator()(Scalar* blockA, const DataMapper& lhs, Index depth,
-                                                            Index rows, Index stride, Index offset) {
+                                                            Index rows, Index stride, Index offset) const {
   typedef typename unpacket_traits<Packet>::half HalfPacket;
   typedef typename unpacket_traits<typename unpacket_traits<Packet>::half>::half QuarterPacket;
   enum {
@@ -2057,12 +2057,12 @@ struct gemm_pack_rhs<Scalar, Index, DataMapper, nr, ColMajor, Conjugate, PanelMo
   typedef typename DataMapper::LinearMapper LinearMapper;
   enum { PacketSize = packet_traits<Scalar>::size };
   EIGEN_DONT_INLINE void operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                                    Index offset = 0);
+                                    Index offset = 0) const;
 };
 
 template <typename Scalar, typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
 EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>::operator()(
-    Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) {
+    Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset) const {
   EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS COLMAJOR");
   EIGEN_UNUSED_VARIABLE(stride);
   EIGEN_UNUSED_VARIABLE(offset);
@@ -2244,7 +2244,7 @@ struct gemm_pack_rhs<Scalar, Index, DataMapper, nr, RowMajor, Conjugate, PanelMo
     QuarterPacketSize = unpacket_traits<QuarterPacket>::size
   };
   EIGEN_DONT_INLINE void operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride = 0,
-                                    Index offset = 0) {
+                                    Index offset = 0) const {
     EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS ROWMAJOR");
     EIGEN_UNUSED_VARIABLE(stride);
     EIGEN_UNUSED_VARIABLE(offset);
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index c0aae46ff..c84895554 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -222,7 +222,7 @@ class ComplexEigenSolver {
   }
 
   /** \brief Returns the maximum number of iterations. */
-  Index getMaxIterations() { return m_schur.getMaxIterations(); }
+  Index getMaxIterations() const { return m_schur.getMaxIterations(); }
 
  protected:
   EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index b9193a2f6..86d6349c1 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -233,7 +233,7 @@ class ComplexSchur {
   }
 
   /** \brief Returns the maximum number of iterations. */
-  Index getMaxIterations() { return m_maxIters; }
+  Index getMaxIterations() const { return m_maxIters; }
 
   /** \brief Maximum number of iterations per row.
    *
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 9dba7bd18..f43991207 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -292,7 +292,7 @@ class EigenSolver {
   }
 
   /** \brief Returns the maximum number of iterations. */
-  Index getMaxIterations() { return m_realSchur.getMaxIterations(); }
+  Index getMaxIterations() const { return m_realSchur.getMaxIterations(); }
 
  private:
   void doComputeEigenvectors();
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index fda3b2eb7..6e50a432c 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -207,7 +207,7 @@ class RealSchur {
   }
 
   /** \brief Returns the maximum number of iterations. */
-  Index getMaxIterations() { return m_maxIters; }
+  Index getMaxIterations() const { return m_maxIters; }
 
   /** \brief Maximum number of iterations per row.
    *
diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
index 2be9d436a..d4a187ea3 100644
--- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
+++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
@@ -93,7 +93,7 @@ class DiagonalPreconditioner {
     return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived());
   }
 
-  ComputationInfo info() { return Success; }
+  ComputationInfo info() const { return Success; }
 
  protected:
   Vector m_invdiag;
@@ -166,7 +166,7 @@ class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<Scalar_>
     return factorize(mat);
   }
 
-  ComputationInfo info() { return Success; }
+  ComputationInfo info() const { return Success; }
 
  protected:
 };
@@ -205,7 +205,7 @@ class IdentityPreconditioner {
     return b;
   }
 
-  ComputationInfo info() { return Success; }
+  ComputationInfo info() const { return Success; }
 };
 
 }  // end namespace Eigen
diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h
index 282d09f7d..8ede4de13 100644
--- a/Eigen/src/OrderingMethods/Ordering.h
+++ b/Eigen/src/OrderingMethods/Ordering.h
@@ -53,7 +53,7 @@ class AMDOrdering {
    * This routine is much faster if the input matrix is column-major
    */
   template <typename MatrixType>
-  void operator()(const MatrixType& mat, PermutationType& perm) {
+  void operator()(const MatrixType& mat, PermutationType& perm) const {
     // Compute the symmetric pattern
     SparseMatrix<typename MatrixType::Scalar, ColMajor, StorageIndex> symm;
     internal::ordering_helper_at_plus_a(mat, symm);
@@ -65,7 +65,7 @@ class AMDOrdering {
 
   /** Compute the permutation with a selfadjoint matrix */
   template <typename SrcType, unsigned int SrcUpLo>
-  void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm) {
+  void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm) const {
     SparseMatrix<typename SrcType::Scalar, ColMajor, StorageIndex> C;
     C = mat;
 
@@ -90,7 +90,7 @@ class NaturalOrdering {
 
   /** Compute the permutation vector from a column-major sparse matrix */
   template <typename MatrixType>
-  void operator()(const MatrixType& /*mat*/, PermutationType& perm) {
+  void operator()(const MatrixType& /*mat*/, PermutationType& perm) const {
     perm.resize(0);
   }
 };
@@ -113,7 +113,7 @@ class COLAMDOrdering {
    * \warning The input sparse matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()).
    */
   template <typename MatrixType>
-  void operator()(const MatrixType& mat, PermutationType& perm) {
+  void operator()(const MatrixType& mat, PermutationType& perm) const {
     eigen_assert(mat.isCompressed() &&
                  "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it "
                  "to COLAMDOrdering");
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index f1662c51d..9b12f4902 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -368,7 +368,7 @@ class SparseLU : public SparseSolverBase<SparseLU<MatrixType_, OrderingType_>>,
    *
    * \sa logAbsDeterminant(), signDeterminant()
    */
-  Scalar absDeterminant() {
+  Scalar absDeterminant() const {
     using std::abs;
     eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
     // Initialize with the determinant of the row matrix
@@ -420,7 +420,7 @@ class SparseLU : public SparseSolverBase<SparseLU<MatrixType_, OrderingType_>>,
    *
    * \sa absDeterminant(), logAbsDeterminant()
    */
-  Scalar signDeterminant() {
+  Scalar signDeterminant() const {
     eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
     // Initialize with the determinant of the row matrix
     Index det = 1;
@@ -446,7 +446,7 @@ class SparseLU : public SparseSolverBase<SparseLU<MatrixType_, OrderingType_>>,
    *
    * \sa absDeterminant(), logAbsDeterminant()
    */
-  Scalar determinant() {
+  Scalar determinant() const {
     eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
     // Initialize with the determinant of the row matrix
     Scalar det = Scalar(1.);
diff --git a/unsupported/Eigen/src/BVH/KdBVH.h b/unsupported/Eigen/src/BVH/KdBVH.h
index d421e6f6b..631aef749 100644
--- a/unsupported/Eigen/src/BVH/KdBVH.h
+++ b/unsupported/Eigen/src/BVH/KdBVH.h
@@ -33,7 +33,7 @@ struct vector_int_pair {
 // iterator range or using bounding_box in a unified way
 template <typename ObjectList, typename VolumeList, typename BoxIter>
 struct get_boxes_helper {
-  void operator()(const ObjectList &objects, BoxIter boxBegin, BoxIter boxEnd, VolumeList &outBoxes) {
+  void operator()(const ObjectList &objects, BoxIter boxBegin, BoxIter boxEnd, VolumeList &outBoxes) const {
     outBoxes.insert(outBoxes.end(), boxBegin, boxEnd);
     eigen_assert(outBoxes.size() == objects.size());
     EIGEN_ONLY_USED_FOR_DEBUG(objects);
@@ -42,7 +42,7 @@ struct get_boxes_helper {
 
 template <typename ObjectList, typename VolumeList>
 struct get_boxes_helper<ObjectList, VolumeList, int> {
-  void operator()(const ObjectList &objects, int, int, VolumeList &outBoxes) {
+  void operator()(const ObjectList &objects, int, int, VolumeList &outBoxes) const {
     outBoxes.reserve(objects.size());
     for (int i = 0; i < (int)objects.size(); ++i) outBoxes.push_back(bounding_box(objects[i]));
   }
diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
index 74edd740b..2c1050944 100644
--- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
@@ -160,7 +160,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<MatrixType_, Preconditioner_> >
   /**
    * Get the restart value
    */
-  Index restart() { return m_restart; }
+  Index restart() const { return m_restart; }
 
   /**
    * Set the restart value (default is 30)
@@ -178,7 +178,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<MatrixType_, Preconditioner_> >
   /**
    * Get the size of the deflation subspace size
    */
-  Index deflSize() { return m_r; }
+  Index deflSize() const { return m_r; }
 
   /**
    * Set the maximum size of the deflation subspace
diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
index 238d0ec6e..f2c85ffdd 100644
--- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
@@ -291,7 +291,7 @@ class GMRES : public IterativeSolverBase<GMRES<MatrixType_, Preconditioner_> > {
 
   /** Get the number of iterations after that a restart is performed.
    */
-  Index get_restart() { return m_restart; }
+  Index get_restart() const { return m_restart; }
 
   /** Set the number of iterations after that a restart is performed.
    *  \param restart   number of iterations for a restarti, default is 30.
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h
index c14ad63be..eea424caa 100644
--- a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h
+++ b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h
@@ -185,22 +185,22 @@ class LevenbergMarquardt : internal::no_assignment_operator {
   FVectorType &diag() { return m_diag; }
 
   /** \returns the number of iterations performed */
-  Index iterations() { return m_iter; }
+  Index iterations() const { return m_iter; }
 
   /** \returns the number of functions evaluation */
-  Index nfev() { return m_nfev; }
+  Index nfev() const { return m_nfev; }
 
   /** \returns the number of jacobian evaluation */
-  Index njev() { return m_njev; }
+  Index njev() const { return m_njev; }
 
   /** \returns the norm of current vector function */
-  RealScalar fnorm() { return m_fnorm; }
+  RealScalar fnorm() const { return m_fnorm; }
 
   /** \returns the norm of the gradient of the error */
-  RealScalar gnorm() { return m_gnorm; }
+  RealScalar gnorm() const { return m_gnorm; }
 
   /** \returns the LevenbergMarquardt parameter */
-  RealScalar lm_param(void) { return m_par; }
+  RealScalar lm_param(void) const { return m_par; }
 
   /** \returns a reference to the  current vector function
    */
@@ -217,7 +217,7 @@ class LevenbergMarquardt : internal::no_assignment_operator {
 
   /** the permutation used in the QR factorization
    */
-  PermutationType permutation() { return m_permutation; }
+  PermutationType permutation() const { return m_permutation; }
 
   /**
    * \brief Reports whether the minimization was successful
diff --git a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
index 5f156ac85..6503e6226 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
@@ -105,7 +105,7 @@ class LevenbergMarquardt {
   Scalar fnorm, gnorm;
   bool useExternalScaling;
 
-  Scalar lm_param(void) { return par; }
+  Scalar lm_param(void) const { return par; }
 
  private:
   FunctorType &functor;
diff --git a/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
index f718d78ab..d991e7562 100644
--- a/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
+++ b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
@@ -146,14 +146,14 @@ class MatrixMarketIterator {
 
   inline std::string& matname() { return m_matname; }
 
-  inline int sym() { return m_sym; }
+  inline int sym() const { return m_sym; }
 
-  bool hasRhs() { return m_hasRhs; }
-  bool hasrefX() { return m_hasrefX; }
-  bool isFolderValid() { return bool(m_folder_id); }
+  bool hasRhs() const { return m_hasRhs; }
+  bool hasrefX() const { return m_hasrefX; }
+  bool isFolderValid() const { return bool(m_folder_id); }
 
  protected:
-  inline bool Fileexists(std::string file) {
+  inline bool Fileexists(std::string file) const {
     std::ifstream file_id(file.c_str());
     if (!file_id.good()) {
       return false;
diff --git a/unsupported/Eigen/src/SparseExtra/SparseInverse.h b/unsupported/Eigen/src/SparseExtra/SparseInverse.h
index 142cc8f12..b9e1c55d7 100644
--- a/unsupported/Eigen/src/SparseExtra/SparseInverse.h
+++ b/unsupported/Eigen/src/SparseExtra/SparseInverse.h
@@ -38,7 +38,7 @@ class KahanSum {
   Scalar _correction{};
 
  public:
-  Scalar value() { return _sum; }
+  Scalar value() const { return _sum; }
 
   void operator+=(Scalar increment) {
     const Scalar correctedIncrement = increment + _correction;
@@ -59,7 +59,7 @@ class FABSum {
   Index _blockUsed{};
 
  public:
-  Scalar value() { return _block.topRows(_blockUsed).sum() + _totalSum.value(); }
+  Scalar value() const { return _block.topRows(_blockUsed).sum() + _totalSum.value(); }
 
   void operator+=(Scalar increment) {
     _block(_blockUsed++, 0) = increment;
diff --git a/unsupported/Eigen/src/Tensor/TensorContractionThreadPool.h b/unsupported/Eigen/src/Tensor/TensorContractionThreadPool.h
index 4ca98fd06..6a6cc35a5 100644
--- a/unsupported/Eigen/src/Tensor/TensorContractionThreadPool.h
+++ b/unsupported/Eigen/src/Tensor/TensorContractionThreadPool.h
@@ -281,7 +281,7 @@ struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgT
   // Dummy struct to represent an empty DoneCallback.
 
   struct NoCallback {
-    void operator()() { eigen_assert(false && "NoCallback should never be called"); }
+    void operator()() const { eigen_assert(false && "NoCallback should never be called"); }
   };
 
   // ------------------------------------------------------------------------ //
diff --git a/unsupported/Eigen/src/Tensor/TensorForcedEval.h b/unsupported/Eigen/src/Tensor/TensorForcedEval.h
index 40f5296da..e945dad51 100644
--- a/unsupported/Eigen/src/Tensor/TensorForcedEval.h
+++ b/unsupported/Eigen/src/Tensor/TensorForcedEval.h
@@ -72,7 +72,7 @@ namespace internal {
 template <typename Device, typename CoeffReturnType>
 struct non_integral_type_placement_new {
   template <typename StorageType>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(Index numValues, StorageType m_buffer) {
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(Index numValues, StorageType m_buffer) const {
     // Initialize non-trivially constructible types.
     if (!internal::is_arithmetic<CoeffReturnType>::value) {
       for (Index i = 0; i < numValues; ++i) new (m_buffer + i) CoeffReturnType();
@@ -86,7 +86,7 @@ struct non_integral_type_placement_new {
 template <typename CoeffReturnType>
 struct non_integral_type_placement_new<Eigen::SyclDevice, CoeffReturnType> {
   template <typename StorageType>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(Index, StorageType) {}
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(Index, StorageType) const {}
 };
 }  // end namespace internal
 
diff --git a/unsupported/Eigen/src/Tensor/TensorScan.h b/unsupported/Eigen/src/Tensor/TensorScan.h
index 39d3558b8..f3e5bd4ce 100644
--- a/unsupported/Eigen/src/Tensor/TensorScan.h
+++ b/unsupported/Eigen/src/Tensor/TensorScan.h
@@ -143,7 +143,7 @@ EIGEN_STRONG_INLINE void ReducePacket(Self& self, Index offset, typename Self::C
 
 template <typename Self, bool Vectorize, bool Parallel>
 struct ReduceBlock {
-  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) {
+  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) const {
     for (Index idx2 = 0; idx2 < self.stride(); idx2++) {
       // Calculate the starting offset for the scan
       Index offset = idx1 + idx2;
@@ -155,7 +155,7 @@ struct ReduceBlock {
 // Specialization for vectorized reduction.
 template <typename Self>
 struct ReduceBlock<Self, /*Vectorize=*/true, /*Parallel=*/false> {
-  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) {
+  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) const {
     using Packet = typename Self::PacketReturnType;
     const int PacketSize = internal::unpacket_traits<Packet>::size;
     Index idx2 = 0;
@@ -204,7 +204,7 @@ EIGEN_STRONG_INLINE Index AdjustBlockSize(Index item_size, Index block_size) {
 
 template <typename Self>
 struct ReduceBlock<Self, /*Vectorize=*/true, /*Parallel=*/true> {
-  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) {
+  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) const {
     using Scalar = typename Self::CoeffReturnType;
     using Packet = typename Self::PacketReturnType;
     const int PacketSize = internal::unpacket_traits<Packet>::size;
@@ -243,7 +243,7 @@ struct ReduceBlock<Self, /*Vectorize=*/true, /*Parallel=*/true> {
 
 template <typename Self>
 struct ReduceBlock<Self, /*Vectorize=*/false, /*Parallel=*/true> {
-  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) {
+  EIGEN_STRONG_INLINE void operator()(Self& self, Index idx1, typename Self::CoeffReturnType* data) const {
     using Scalar = typename Self::CoeffReturnType;
     self.device().parallelFor(
         self.stride(), TensorOpCost(self.size(), self.size(), 16 * self.size()),
@@ -261,7 +261,7 @@ struct ReduceBlock<Self, /*Vectorize=*/false, /*Parallel=*/true> {
 // Specialization for multi-threaded execution.
 template <typename Self, typename Reducer, bool Vectorize>
 struct ScanLauncher<Self, Reducer, ThreadPoolDevice, Vectorize> {
-  void operator()(Self& self, typename Self::CoeffReturnType* data) {
+  void operator()(Self& self, typename Self::CoeffReturnType* data) const {
     using Scalar = typename Self::CoeffReturnType;
     using Packet = typename Self::PacketReturnType;
     const int PacketSize = internal::unpacket_traits<Packet>::size;
@@ -333,7 +333,7 @@ __global__ EIGEN_HIP_LAUNCH_BOUNDS_1024 void ScanKernel(Self self, Index total_s
 
 template <typename Self, typename Reducer, bool Vectorize>
 struct ScanLauncher<Self, Reducer, GpuDevice, Vectorize> {
-  void operator()(const Self& self, typename Self::CoeffReturnType* data) {
+  void operator()(const Self& self, typename Self::CoeffReturnType* data) const {
     Index total_size = internal::array_prod(self.dimensions());
     Index num_blocks = (total_size / self.size() + 63) / 64;
     Index block_size = 64;