diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 8767e1f50..605a05e84 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -328,10 +328,10 @@ template<typename Derived> class MatrixBase
 
 /////////// LU module ///////////
 
-    inline const FullPivLU<PlainObject> fullPivLu() const;
-    inline const PartialPivLU<PlainObject> partialPivLu() const;
+    template<typename PermutationIndex = DefaultPermutationIndex> inline const FullPivLU<PlainObject, PermutationIndex> fullPivLu() const;
+    template<typename PermutationIndex = DefaultPermutationIndex> inline const PartialPivLU<PlainObject, PermutationIndex> partialPivLu() const;
 
-    inline const PartialPivLU<PlainObject> lu() const;
+    template<typename PermutationIndex = DefaultPermutationIndex> inline const PartialPivLU<PlainObject, PermutationIndex> lu() const;
 
     EIGEN_DEVICE_FUNC
     inline const Inverse<Derived> inverse() const;
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index a00b14914..8963019a9 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -270,8 +270,8 @@ typedef lapack_int DefaultPermutationIndex;
 typedef int DefaultPermutationIndex;
 #endif
 
-template<typename MatrixType> class FullPivLU;
-template<typename MatrixType> class PartialPivLU;
+template<typename MatrixType, typename PermutationIndex = DefaultPermutationIndex> class FullPivLU;
+template<typename MatrixType, typename PermutationIndex = DefaultPermutationIndex> class PartialPivLU;
 namespace internal {
 template<typename MatrixType> struct inverse_impl;
 }
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 259b549b9..d02dcc10f 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -15,12 +15,12 @@
 namespace Eigen {
 
 namespace internal {
-template<typename MatrixType_> struct traits<FullPivLU<MatrixType_> >
+template<typename MatrixType_, typename PermutationIndex_> struct traits<FullPivLU<MatrixType_, PermutationIndex_> >
  : traits<MatrixType_>
 {
   typedef MatrixXpr XprKind;
   typedef SolverStorage StorageKind;
-  typedef int StorageIndex;
+  typedef PermutationIndex_ StorageIndex;
   enum { Flags = 0 };
 };
 
@@ -59,8 +59,8 @@ template<typename MatrixType_> struct traits<FullPivLU<MatrixType_> >
   *
   * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
   */
-template<typename MatrixType_> class FullPivLU
-  : public SolverBase<FullPivLU<MatrixType_> >
+template<typename MatrixType_, typename PermutationIndex_> class FullPivLU
+  : public SolverBase<FullPivLU<MatrixType_, PermutationIndex_> >
 {
   public:
     typedef MatrixType_ MatrixType;
@@ -72,10 +72,11 @@ template<typename MatrixType_> class FullPivLU
       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
-    typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
-    typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
-    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
-    typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
+    using PermutationIndex = PermutationIndex_;
+    typedef typename internal::plain_row_type<MatrixType, PermutationIndex>::type IntRowVectorType;
+    typedef typename internal::plain_col_type<MatrixType, PermutationIndex>::type IntColVectorType;
+    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, PermutationIndex> PermutationQType;
+    typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> PermutationPType;
     typedef typename MatrixType::PlainObject PlainObject;
 
     /**
@@ -437,14 +438,14 @@ template<typename MatrixType_> class FullPivLU
     bool m_isInitialized, m_usePrescribedThreshold;
 };
 
-template<typename MatrixType>
-FullPivLU<MatrixType>::FullPivLU()
+template<typename MatrixType, typename PermutationIndex>
+FullPivLU<MatrixType, PermutationIndex>::FullPivLU()
   : m_isInitialized(false), m_usePrescribedThreshold(false)
 {
 }
 
-template<typename MatrixType>
-FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
+template<typename MatrixType, typename PermutationIndex>
+FullPivLU<MatrixType, PermutationIndex>::FullPivLU(Index rows, Index cols)
   : m_lu(rows, cols),
     m_p(rows),
     m_q(cols),
@@ -455,9 +456,9 @@ FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
 {
 }
 
-template<typename MatrixType>
+template<typename MatrixType, typename PermutationIndex>
 template<typename InputType>
-FullPivLU<MatrixType>::FullPivLU(const EigenBase<InputType>& matrix)
+FullPivLU<MatrixType, PermutationIndex>::FullPivLU(const EigenBase<InputType>& matrix)
   : m_lu(matrix.rows(), matrix.cols()),
     m_p(matrix.rows()),
     m_q(matrix.cols()),
@@ -469,9 +470,9 @@ FullPivLU<MatrixType>::FullPivLU(const EigenBase<InputType>& matrix)
   compute(matrix.derived());
 }
 
-template<typename MatrixType>
+template<typename MatrixType, typename PermutationIndex>
 template<typename InputType>
-FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix)
+FullPivLU<MatrixType, PermutationIndex>::FullPivLU(EigenBase<InputType>& matrix)
   : m_lu(matrix.derived()),
     m_p(matrix.rows()),
     m_q(matrix.cols()),
@@ -483,11 +484,10 @@ FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix)
   computeInPlace();
 }
 
-template<typename MatrixType>
-void FullPivLU<MatrixType>::computeInPlace()
+template<typename MatrixType, typename PermutationIndex>
+void FullPivLU<MatrixType, PermutationIndex>::computeInPlace()
 {
-  // the permutations are stored as int indices, so just to be sure:
-  eigen_assert(m_lu.rows()<=NumTraits<int>::highest() && m_lu.cols()<=NumTraits<int>::highest());
+  eigen_assert(m_lu.rows()<=NumTraits<PermutationIndex>::highest() && m_lu.cols()<=NumTraits<PermutationIndex>::highest());
 
   m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
 
@@ -574,8 +574,8 @@ void FullPivLU<MatrixType>::computeInPlace()
   m_isInitialized = true;
 }
 
-template<typename MatrixType>
-typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
+template<typename MatrixType, typename PermutationIndex>
+typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType, PermutationIndex>::determinant() const
 {
   eigen_assert(m_isInitialized && "LU is not initialized.");
   eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
@@ -585,8 +585,8 @@ typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant
 /** \returns the matrix represented by the decomposition,
  * i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$.
  * This function is provided for debug purposes. */
-template<typename MatrixType>
-MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
+template<typename MatrixType, typename PermutationIndex>
+MatrixType FullPivLU<MatrixType, PermutationIndex>::reconstructedMatrix() const
 {
   eigen_assert(m_isInitialized && "LU is not initialized.");
   const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
@@ -610,11 +610,12 @@ MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
 /********* Implementation of kernel() **************************************************/
 
 namespace internal {
-template<typename MatrixType_>
-struct kernel_retval<FullPivLU<MatrixType_> >
-  : kernel_retval_base<FullPivLU<MatrixType_> >
+template<typename MatrixType_, typename PermutationIndex_>
+struct kernel_retval<FullPivLU<MatrixType_, PermutationIndex_> >
+  : kernel_retval_base<FullPivLU<MatrixType_, PermutationIndex_> >
 {
-  EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<MatrixType_>)
+  using DecompositionType = FullPivLU<MatrixType_, PermutationIndex_>;
+  EIGEN_MAKE_KERNEL_HELPERS(DecompositionType)
 
   enum { MaxSmallDimAtCompileTime = min_size_prefer_fixed(
               MatrixType::MaxColsAtCompileTime,
@@ -696,11 +697,12 @@ struct kernel_retval<FullPivLU<MatrixType_> >
 
 /***** Implementation of image() *****************************************************/
 
-template<typename MatrixType_>
-struct image_retval<FullPivLU<MatrixType_> >
-  : image_retval_base<FullPivLU<MatrixType_> >
+template<typename MatrixType_, typename PermutationIndex_>
+struct image_retval<FullPivLU<MatrixType_, PermutationIndex_> >
+  : image_retval_base<FullPivLU<MatrixType_, PermutationIndex_> >
 {
-  EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<MatrixType_>)
+  using DecompositionType = FullPivLU<MatrixType_, PermutationIndex_>;
+  EIGEN_MAKE_IMAGE_HELPERS(DecompositionType)
 
   enum { MaxSmallDimAtCompileTime = min_size_prefer_fixed(
               MatrixType::MaxColsAtCompileTime,
@@ -737,9 +739,9 @@ struct image_retval<FullPivLU<MatrixType_> >
 } // end namespace internal
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-template<typename MatrixType_>
+template<typename MatrixType_, typename PermutationIndex_>
 template<typename RhsType, typename DstType>
-void FullPivLU<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
+void FullPivLU<MatrixType_, PermutationIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
 {
   /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
   * So we proceed as follows:
@@ -784,9 +786,9 @@ void FullPivLU<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
     dst.row(permutationQ().indices().coeff(i)).setZero();
 }
 
-template<typename MatrixType_>
+template<typename MatrixType_, typename PermutationIndex_>
 template<bool Conjugate, typename RhsType, typename DstType>
-void FullPivLU<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
+void FullPivLU<MatrixType_, PermutationIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
 {
   /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1},
    * and since permutations are real and unitary, we can write this
@@ -842,10 +844,10 @@ namespace internal {
 
 
 /***** Implementation of inverse() *****************************************************/
-template<typename DstXprType, typename MatrixType>
-struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType>::Scalar>, Dense2Dense>
+template<typename DstXprType, typename MatrixType, typename PermutationIndex>
+struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType, PermutationIndex> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType, PermutationIndex>::Scalar>, Dense2Dense>
 {
-  typedef FullPivLU<MatrixType> LuType;
+  typedef FullPivLU<MatrixType, PermutationIndex> LuType;
   typedef Inverse<LuType> SrcXprType;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename MatrixType::Scalar> &)
   {
@@ -863,10 +865,11 @@ struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_
   * \sa class FullPivLU
   */
 template<typename Derived>
-inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
+template<typename PermutationIndex>
+inline const FullPivLU<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
 MatrixBase<Derived>::fullPivLu() const
 {
-  return FullPivLU<PlainObject>(eval());
+  return FullPivLU<PlainObject, PermutationIndex>(eval());
 }
 
 } // end namespace Eigen
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 137739892..bee743794 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -16,12 +16,12 @@
 namespace Eigen {
 
 namespace internal {
-template<typename MatrixType_> struct traits<PartialPivLU<MatrixType_> >
+template<typename MatrixType_, typename PermutationIndex_> struct traits<PartialPivLU<MatrixType_, PermutationIndex_> >
  : traits<MatrixType_>
 {
   typedef MatrixXpr XprKind;
   typedef SolverStorage StorageKind;
-  typedef int StorageIndex;
+  typedef PermutationIndex_ StorageIndex;
   typedef traits<MatrixType_> BaseTraits;
   enum {
     Flags = BaseTraits::Flags & RowMajorBit,
@@ -75,8 +75,8 @@ struct enable_if_ref<Ref<T>,Derived> {
   *
   * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
   */
-template<typename MatrixType_> class PartialPivLU
-  : public SolverBase<PartialPivLU<MatrixType_> >
+template<typename MatrixType_, typename PermutationIndex_> class PartialPivLU
+  : public SolverBase<PartialPivLU<MatrixType_, PermutationIndex_> >
 {
   public:
 
@@ -89,8 +89,9 @@ template<typename MatrixType_> class PartialPivLU
       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
-    typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
-    typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
+    using PermutationIndex = PermutationIndex_;
+    typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> PermutationType;
+    typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime, PermutationIndex> TranspositionType;
     typedef typename MatrixType::PlainObject PlainObject;
 
     /**
@@ -279,8 +280,8 @@ template<typename MatrixType_> class PartialPivLU
     bool m_isInitialized;
 };
 
-template<typename MatrixType>
-PartialPivLU<MatrixType>::PartialPivLU()
+template<typename MatrixType, typename PermutationIndex>
+PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU()
   : m_lu(),
     m_p(),
     m_rowsTranspositions(),
@@ -290,8 +291,8 @@ PartialPivLU<MatrixType>::PartialPivLU()
 {
 }
 
-template<typename MatrixType>
-PartialPivLU<MatrixType>::PartialPivLU(Index size)
+template<typename MatrixType, typename PermutationIndex>
+PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU(Index size)
   : m_lu(size, size),
     m_p(size),
     m_rowsTranspositions(size),
@@ -301,9 +302,9 @@ PartialPivLU<MatrixType>::PartialPivLU(Index size)
 {
 }
 
-template<typename MatrixType>
+template<typename MatrixType, typename PermutationIndex>
 template<typename InputType>
-PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix)
+PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU(const EigenBase<InputType>& matrix)
   : m_lu(matrix.rows(),matrix.cols()),
     m_p(matrix.rows()),
     m_rowsTranspositions(matrix.rows()),
@@ -314,9 +315,9 @@ PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix)
   compute(matrix.derived());
 }
 
-template<typename MatrixType>
+template<typename MatrixType, typename PermutationIndex>
 template<typename InputType>
-PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix)
+PartialPivLU<MatrixType, PermutationIndex>::PartialPivLU(EigenBase<InputType>& matrix)
   : m_lu(matrix.derived()),
     m_p(matrix.rows()),
     m_rowsTranspositions(matrix.rows()),
@@ -520,11 +521,10 @@ void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, t
 
 } // end namespace internal
 
-template<typename MatrixType>
-void PartialPivLU<MatrixType>::compute()
+template<typename MatrixType, typename PermutationIndex>
+void PartialPivLU<MatrixType, PermutationIndex>::compute()
 {
-  // the row permutation is stored as int indices, so just to be sure:
-  eigen_assert(m_lu.rows()<NumTraits<int>::highest());
+  eigen_assert(m_lu.rows()<NumTraits<PermutationIndex>::highest());
 
   if(m_lu.cols()>0)
     m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
@@ -545,8 +545,8 @@ void PartialPivLU<MatrixType>::compute()
   m_isInitialized = true;
 }
 
-template<typename MatrixType>
-typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
+template<typename MatrixType, typename PermutationIndex>
+typename PartialPivLU<MatrixType, PermutationIndex>::Scalar PartialPivLU<MatrixType, PermutationIndex>::determinant() const
 {
   eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
   return Scalar(m_det_p) * m_lu.diagonal().prod();
@@ -555,8 +555,8 @@ typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant(
 /** \returns the matrix represented by the decomposition,
  * i.e., it returns the product: P^{-1} L U.
  * This function is provided for debug purpose. */
-template<typename MatrixType>
-MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
+template<typename MatrixType, typename PermutationIndex>
+MatrixType PartialPivLU<MatrixType, PermutationIndex>::reconstructedMatrix() const
 {
   eigen_assert(m_isInitialized && "LU is not initialized.");
   // LU
@@ -574,10 +574,10 @@ MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
 namespace internal {
 
 /***** Implementation of inverse() *****************************************************/
-template<typename DstXprType, typename MatrixType>
-struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
+template<typename DstXprType, typename MatrixType, typename PermutationIndex>
+struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType, PermutationIndex> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType, PermutationIndex>::Scalar>, Dense2Dense>
 {
-  typedef PartialPivLU<MatrixType> LuType;
+  typedef PartialPivLU<MatrixType, PermutationIndex> LuType;
   typedef Inverse<LuType> SrcXprType;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
   {
@@ -595,10 +595,11 @@ struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assi
   * \sa class PartialPivLU
   */
 template<typename Derived>
-inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
+template<typename PermutationIndex>
+inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
 MatrixBase<Derived>::partialPivLu() const
 {
-  return PartialPivLU<PlainObject>(eval());
+  return PartialPivLU<PlainObject, PermutationIndex>(eval());
 }
 
 /** \lu_module
@@ -610,10 +611,11 @@ MatrixBase<Derived>::partialPivLu() const
   * \sa class PartialPivLU
   */
 template<typename Derived>
-inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
+template<typename PermutationIndex>
+inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject, PermutationIndex>
 MatrixBase<Derived>::lu() const
 {
-  return PartialPivLU<PlainObject>(eval());
+  return PartialPivLU<PlainObject, PermutationIndex>(eval());
 }
 
 } // end namespace Eigen
diff --git a/test/lu.cpp b/test/lu.cpp
index 1bbadcbf0..3369f82a2 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -19,7 +19,6 @@ typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
 
 template<typename MatrixType> void lu_non_invertible()
 {
-  STATIC_CHECK(( internal::is_same<typename FullPivLU<MatrixType>::StorageIndex,int>::value ));
 
   typedef typename MatrixType::RealScalar RealScalar;
   /* this test covers the following files:
@@ -163,8 +162,6 @@ template<typename MatrixType> void lu_partial_piv(Index size = MatrixType::ColsA
   m1.setRandom();
   PartialPivLU<MatrixType> plu(m1);
 
-  STATIC_CHECK(( internal::is_same<typename PartialPivLU<MatrixType>::StorageIndex,int>::value ));
-
   VERIFY_IS_APPROX(m1, plu.reconstructedMatrix());
 
   check_solverbase<MatrixType, MatrixType>(m1, plu, size, size, size);