diff --git a/Eigen/src/Cholesky/LLT_LAPACKE.h b/Eigen/src/Cholesky/LLT_LAPACKE.h
index 62bc679..a32ff23 100644
--- a/Eigen/src/Cholesky/LLT_LAPACKE.h
+++ b/Eigen/src/Cholesky/LLT_LAPACKE.h
@@ -70,6 +70,7 @@
 
   template<typename Scalar, UpLoType Mode>
   struct lapacke_llt {
+    EIGEN_STATIC_ASSERT(((Mode == Lower) || (Mode == Upper)),MODE_MUST_BE_UPPER_OR_LOWER)
     template<typename MatrixType>
     static Index blocked(MatrixType& m)
     {
@@ -80,10 +81,11 @@
       /* Set up parameters for ?potrf */
       lapack_int size = to_lapack(m.rows());
       lapack_int matrix_order = lapack_storage_of(m);
+      constexpr char uplo = Mode == Upper ? 'U' : 'L';
       Scalar* a = &(m.coeffRef(0,0));
       lapack_int lda = to_lapack(m.outerStride());
 
-      lapack_int info = potrf(matrix_order, translate_mode<Mode>, size, to_lapack(a), lda );
+      lapack_int info = potrf(matrix_order, uplo, size, to_lapack(a), lda );
       info = (info==0) ? -1 : info>0 ? info-1 : size;
       return info;
     }
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index ea2178f..c6e81a7 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -362,9 +362,9 @@
 /////////// QR module ///////////
 
     inline const HouseholderQR<PlainObject> householderQr() const;
-    inline const ColPivHouseholderQR<PlainObject> colPivHouseholderQr() const;
-    inline const FullPivHouseholderQR<PlainObject> fullPivHouseholderQr() const;
-    inline const CompleteOrthogonalDecomposition<PlainObject> completeOrthogonalDecomposition() const;
+    template<typename StorageIndex = DefaultStorageIndex> inline const ColPivHouseholderQR<PlainObject, StorageIndex> colPivHouseholderQr() const;
+    template<typename StorageIndex = int> inline const FullPivHouseholderQR<PlainObject, StorageIndex> fullPivHouseholderQr() const;
+    template<typename StorageIndex = DefaultStorageIndex> inline const CompleteOrthogonalDecomposition<PlainObject, StorageIndex> completeOrthogonalDecomposition() const;
 
 /////////// Eigenvalues module ///////////
 
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 8f87c4a..1ebe0cd 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -252,15 +252,22 @@
 template<typename MatrixType,int RowFactor,int ColFactor> class Replicate;
 template<typename MatrixType, int Direction = BothDirections> class Reverse;
 
+#if defined(EIGEN_USE_LAPACKE)
+// Lapacke interface requires StorageIndex to be lapack_int
+typedef lapack_int DefaultStorageIndex;
+#else
+typedef int DefaultStorageIndex;
+#endif
+
 template<typename MatrixType> class FullPivLU;
 template<typename MatrixType> class PartialPivLU;
 namespace internal {
 template<typename MatrixType> struct inverse_impl;
 }
 template<typename MatrixType> class HouseholderQR;
-template<typename MatrixType> class ColPivHouseholderQR;
-template<typename MatrixType> class FullPivHouseholderQR;
-template<typename MatrixType> class CompleteOrthogonalDecomposition;
+template<typename MatrixType, typename StorageIndex = DefaultStorageIndex> class ColPivHouseholderQR;
+template<typename MatrixType, typename StorageIndex = int> class FullPivHouseholderQR;
+template<typename MatrixType, typename StorageIndex = DefaultStorageIndex> class CompleteOrthogonalDecomposition;
 template<typename MatrixType> class SVDBase;
 template<typename MatrixType, int Options = 0> class JacobiSVD;
 template<typename MatrixType, int Options = 0> class BDCSVD;
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index c906997..32494df 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -16,12 +16,12 @@
 namespace Eigen {
 
 namespace internal {
-template<typename MatrixType_> struct traits<ColPivHouseholderQR<MatrixType_> >
+template<typename MatrixType_, typename StorageIndex_> struct traits<ColPivHouseholderQR<MatrixType_, StorageIndex_> >
  : traits<MatrixType_>
 {
   typedef MatrixXpr XprKind;
   typedef SolverStorage StorageKind;
-  typedef int StorageIndex;
+  typedef StorageIndex_ StorageIndex;
   enum { Flags = 0 };
 };
 
@@ -50,8 +50,8 @@
   * 
   * \sa MatrixBase::colPivHouseholderQr()
   */
-template<typename MatrixType_> class ColPivHouseholderQR
-        : public SolverBase<ColPivHouseholderQR<MatrixType_> >
+template<typename MatrixType_, typename StorageIndex_> class ColPivHouseholderQR
+        : public SolverBase<ColPivHouseholderQR<MatrixType_, StorageIndex_> >
 {
   public:
 
@@ -65,17 +65,13 @@
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
-    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
+    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, StorageIndex> PermutationType;
     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
     typedef HouseholderSequence<MatrixType,internal::remove_all_t<typename HCoeffsType::ConjugateReturnType>> HouseholderSequenceType;
     typedef typename MatrixType::PlainObject PlainObject;
 
-  private:
-
-    typedef typename PermutationType::StorageIndex PermIndexType;
-
   public:
 
     /**
@@ -104,7 +100,7 @@
     ColPivHouseholderQR(Index rows, Index cols)
       : m_qr(rows, cols),
         m_hCoeffs((std::min)(rows,cols)),
-        m_colsPermutation(PermIndexType(cols)),
+        m_colsPermutation(cols),
         m_colsTranspositions(cols),
         m_temp(cols),
         m_colNormsUpdated(cols),
@@ -128,7 +124,7 @@
     explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
-        m_colsPermutation(PermIndexType(matrix.cols())),
+        m_colsPermutation(matrix.cols()),
         m_colsTranspositions(matrix.cols()),
         m_temp(matrix.cols()),
         m_colNormsUpdated(matrix.cols()),
@@ -149,7 +145,7 @@
     explicit ColPivHouseholderQR(EigenBase<InputType>& matrix)
       : m_qr(matrix.derived()),
         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
-        m_colsPermutation(PermIndexType(matrix.cols())),
+        m_colsPermutation(matrix.cols()),
         m_colsTranspositions(matrix.cols()),
         m_temp(matrix.cols()),
         m_colNormsUpdated(matrix.cols()),
@@ -441,7 +437,7 @@
 
   protected:
 
-    friend class CompleteOrthogonalDecomposition<MatrixType>;
+    friend class CompleteOrthogonalDecomposition<MatrixType, StorageIndex>;
 
     EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
 
@@ -460,8 +456,8 @@
     Index m_det_p;
 };
 
-template<typename MatrixType>
-typename MatrixType::Scalar ColPivHouseholderQR<MatrixType>::determinant() const
+template<typename MatrixType, typename StorageIndex>
+typename MatrixType::Scalar ColPivHouseholderQR<MatrixType, StorageIndex>::determinant() const
 {
   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@@ -470,8 +466,8 @@
   return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
 }
 
-template<typename MatrixType>
-typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
+template<typename MatrixType, typename StorageIndex>
+typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, StorageIndex>::absDeterminant() const
 {
   using std::abs;
   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
@@ -479,8 +475,8 @@
   return abs(m_qr.diagonal().prod());
 }
 
-template<typename MatrixType>
-typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
+template<typename MatrixType, typename StorageIndex>
+typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType, StorageIndex>::logAbsDeterminant() const
 {
   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@@ -493,20 +489,20 @@
   *
   * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
   */
-template<typename MatrixType>
+template<typename MatrixType, typename StorageIndex>
 template<typename InputType>
-ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
+ColPivHouseholderQR<MatrixType, StorageIndex>& ColPivHouseholderQR<MatrixType, StorageIndex>::compute(const EigenBase<InputType>& matrix)
 {
   m_qr = matrix.derived();
   computeInPlace();
   return *this;
 }
 
-template<typename MatrixType>
-void ColPivHouseholderQR<MatrixType>::computeInPlace()
+template<typename MatrixType, typename StorageIndex>
+void ColPivHouseholderQR<MatrixType, StorageIndex>::computeInPlace()
 {
-  // the column permutation is stored as int indices, so just to be sure:
-  eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
+
+  eigen_assert(m_qr.cols()<=NumTraits<StorageIndex>::highest());
 
   using std::abs;
 
@@ -595,18 +591,18 @@
     }
   }
 
-  m_colsPermutation.setIdentity(PermIndexType(cols));
-  for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
-    m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
+  m_colsPermutation.setIdentity(cols);
+  for(Index k = 0; k < size/*m_nonzero_pivots*/; ++k)
+    m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
 
   m_det_p = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;
 }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-template<typename MatrixType_>
+template<typename MatrixType_, typename StorageIndex_>
 template<typename RhsType, typename DstType>
-void ColPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
+void ColPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
 {
   const Index nonzero_pivots = nonzeroPivots();
 
@@ -628,9 +624,9 @@
   for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
 }
 
-template<typename MatrixType_>
+template<typename MatrixType_, typename StorageIndex_>
 template<bool Conjugate, typename RhsType, typename DstType>
-void ColPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
+void ColPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
 {
   const Index nonzero_pivots = nonzeroPivots();
 
@@ -656,10 +652,10 @@
 
 namespace internal {
 
-template<typename DstXprType, typename MatrixType>
-struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
+template<typename DstXprType, typename MatrixType, typename StorageIndex>
+struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType, StorageIndex> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType, StorageIndex>::Scalar>, Dense2Dense>
 {
-  typedef ColPivHouseholderQR<MatrixType> QrType;
+  typedef ColPivHouseholderQR<MatrixType, StorageIndex> QrType;
   typedef Inverse<QrType> SrcXprType;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
   {
@@ -672,8 +668,8 @@
 /** \returns the matrix Q as a sequence of householder transformations.
   * You can extract the meaningful part only by using:
   * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/
-template<typename MatrixType>
-typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
+template<typename MatrixType, typename StorageIndex>
+typename ColPivHouseholderQR<MatrixType, StorageIndex>::HouseholderSequenceType ColPivHouseholderQR<MatrixType, StorageIndex>
   ::householderQ() const
 {
   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
@@ -685,10 +681,11 @@
   * \sa class ColPivHouseholderQR
   */
 template<typename Derived>
-const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
+template<typename StorageIndex>
+const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject, StorageIndex>
 MatrixBase<Derived>::colPivHouseholderQr() const
 {
-  return ColPivHouseholderQR<PlainObject>(eval());
+  return ColPivHouseholderQR<PlainObject, StorageIndex>(eval());
 }
 
 } // end namespace Eigen
diff --git a/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h b/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
index 7652d31..15f0978 100644
--- a/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
+++ b/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
@@ -36,64 +36,110 @@
 
 #include "./InternalHeaderCheck.h"
 
-namespace Eigen { 
+namespace Eigen {
 
-/** \internal Specialization for the data types supported by LAPACKe */
+#if defined(EIGEN_USE_LAPACKE)
 
-#define EIGEN_LAPACKE_QR_COLPIV(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW) \
-template<> template<typename InputType> inline \
-ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >& \
-ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >::compute( \
-              const EigenBase<InputType>& matrix) \
-\
-{ \
-  using std::abs; \
-  typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> MatrixType; \
-  typedef MatrixType::RealScalar RealScalar; \
-  Index rows = matrix.rows();\
-  Index cols = matrix.cols();\
-\
-  m_qr = matrix;\
-  Index size = m_qr.diagonalSize();\
-  m_hCoeffs.resize(size);\
-\
-  m_colsTranspositions.resize(cols);\
-  /*Index number_of_transpositions = 0;*/ \
-\
-  m_nonzero_pivots = 0; \
-  m_maxpivot = RealScalar(0);\
-  m_colsPermutation.resize(cols); \
-  m_colsPermutation.indices().setZero(); \
-\
-  lapack_int lda = internal::convert_index<lapack_int,Index>(m_qr.outerStride()); \
-  lapack_int matrix_order = LAPACKE_COLROW; \
-  LAPACKE_##LAPACKE_PREFIX##geqp3( matrix_order, internal::convert_index<lapack_int,Index>(rows), internal::convert_index<lapack_int,Index>(cols), \
-                              (LAPACKE_TYPE*)m_qr.data(), lda, (lapack_int*)m_colsPermutation.indices().data(), (LAPACKE_TYPE*)m_hCoeffs.data()); \
-  m_isInitialized = true; \
-  m_maxpivot=m_qr.diagonal().cwiseAbs().maxCoeff(); \
-  m_hCoeffs.adjointInPlace(); \
-  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); \
-  lapack_int *perm = m_colsPermutation.indices().data(); \
-  for(Index i=0;i<size;i++) { \
-    m_nonzero_pivots += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);\
-  } \
-  for(Index i=0;i<cols;i++) perm[i]--;\
-\
-  /*m_det_pq = (number_of_transpositions%2) ? -1 : 1;  // TODO: It's not needed now; fix upon availability in Eigen */ \
-\
-  return *this; \
-}
+    template<typename Scalar>
+    inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, Scalar* a, lapack_int lda, lapack_int* jpvt, Scalar* tau);
+    template<>
+    inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, float* a, lapack_int lda, lapack_int* jpvt, float* tau)
+    { return LAPACKE_sgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
+    template<>
+    inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, double* a, lapack_int lda, lapack_int* jpvt, double* tau)
+    { return LAPACKE_dgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
+    template<>
+    inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, lapack_complex_float* a, lapack_int lda, lapack_int* jpvt, lapack_complex_float* tau)
+    { return LAPACKE_cgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
+    template<>
+    inline lapack_int call_geqp3(int matrix_layout, lapack_int m, lapack_int n, lapack_complex_double* a, lapack_int lda, lapack_int* jpvt, lapack_complex_double* tau)
+    { return LAPACKE_zgeqp3(matrix_layout, m, n, a, lda, jpvt, tau); }
 
-EIGEN_LAPACKE_QR_COLPIV(double,   double,        d, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_LAPACKE_QR_COLPIV(float,    float,         s, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_LAPACKE_QR_COLPIV(dcomplex, lapack_complex_double, z, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_LAPACKE_QR_COLPIV(scomplex, lapack_complex_float,  c, ColMajor, LAPACK_COL_MAJOR)
+    template <typename MatrixType>
+    struct ColPivHouseholderQR_LAPACKE_impl {
+      typedef typename MatrixType::Scalar Scalar;
+      typedef typename MatrixType::RealScalar RealScalar;
+      typedef typename internal::lapacke_helpers::translate_type_imp<Scalar>::type LapackeType;
+      static constexpr int LapackeStorage = MatrixType::IsRowMajor ? (LAPACK_ROW_MAJOR) : (LAPACK_COL_MAJOR);
 
-EIGEN_LAPACKE_QR_COLPIV(double,   double,        d, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_LAPACKE_QR_COLPIV(float,    float,         s, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_LAPACKE_QR_COLPIV(dcomplex, lapack_complex_double, z, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_LAPACKE_QR_COLPIV(scomplex, lapack_complex_float,  c, RowMajor, LAPACK_ROW_MAJOR)
+      typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
+      typedef PermutationMatrix<Dynamic, Dynamic, lapack_int> PermutationType;
+      typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
 
-} // end namespace Eigen
+      static void run(MatrixType& qr, HCoeffsType& hCoeffs, PermutationType& colsPermutation,
+                      IntRowVectorType& colsTranspositions, Index& nonzero_pivots, RealScalar& maxpivot,
+                      bool usePrescribedThreshold, RealScalar prescribedThreshold, Index& det_p, bool& isInitialized) {
+        using std::abs;
+        
+        isInitialized = false;
+        hCoeffs.resize(qr.diagonalSize());
+        colsTranspositions.resize(qr.cols());
+        nonzero_pivots = 0;
+        maxpivot = RealScalar(0);
+        colsPermutation.resize(qr.cols());
+        colsPermutation.indices().setZero();
 
-#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H
+        lapack_int rows = internal::lapacke_helpers::to_lapack(qr.rows());
+        lapack_int cols = internal::lapacke_helpers::to_lapack(qr.cols());
+        LapackeType* qr_data = (LapackeType*)(qr.data());
+        lapack_int lda = internal::lapacke_helpers::to_lapack(qr.outerStride());
+        lapack_int* perm_data = colsPermutation.indices().data();
+        LapackeType* hCoeffs_data = (LapackeType*)(hCoeffs.data());
+
+        lapack_int info = call_geqp3(LapackeStorage, rows, cols, qr_data, lda, perm_data, hCoeffs_data);
+        if (info != 0) return;
+
+        maxpivot = qr.diagonal().cwiseAbs().maxCoeff();
+        hCoeffs.adjointInPlace();
+        RealScalar defaultThreshold = NumTraits<RealScalar>::epsilon() * RealScalar(qr.diagonalSize());
+        RealScalar threshold = usePrescribedThreshold ? prescribedThreshold : defaultThreshold;
+        RealScalar premultiplied_threshold = abs(maxpivot) * threshold;
+        nonzero_pivots = (qr.diagonal().cwiseAbs().array() > premultiplied_threshold).count();
+        colsPermutation.indices().array() -= 1;
+        det_p = colsPermutation.determinant();
+        isInitialized = true;
+      };
+    };
+
+    typedef Matrix<float, Dynamic, Dynamic, RowMajor> MatrixXfR;
+    typedef Matrix<double, Dynamic, Dynamic, RowMajor> MatrixXdR;
+    typedef Matrix<scomplex, Dynamic, Dynamic, RowMajor> MatrixXcfR;
+    typedef Matrix<dcomplex, Dynamic, Dynamic, RowMajor> MatrixXcdR;
+
+    template <> inline void ColPivHouseholderQR<MatrixXf, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+    template <> inline void ColPivHouseholderQR<MatrixXd, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+    template <> inline void ColPivHouseholderQR<MatrixXcf, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+    template <> inline void ColPivHouseholderQR<MatrixXcd, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+    template <> inline void ColPivHouseholderQR<MatrixXfR, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+    template <> inline void ColPivHouseholderQR<MatrixXdR, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+    template <> inline void ColPivHouseholderQR<MatrixXcfR, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+    template <> inline void ColPivHouseholderQR<MatrixXcdR, lapack_int>::computeInPlace() {
+      ColPivHouseholderQR_LAPACKE_impl<MatrixType>::run(m_qr, m_hCoeffs, m_colsPermutation, m_colsTranspositions,
+                                                        m_nonzero_pivots, m_maxpivot, m_usePrescribedThreshold,
+                                                        m_prescribedThreshold, m_det_p, m_isInitialized); }
+
+#endif
+}  // end namespace Eigen
+
+#endif  // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H
diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
index 02583a2..b3b07c8 100644
--- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h
+++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
@@ -15,12 +15,12 @@
 namespace Eigen {
 
 namespace internal {
-template <typename MatrixType_>
-struct traits<CompleteOrthogonalDecomposition<MatrixType_> >
+template <typename MatrixType_, typename StorageIndex_>
+struct traits<CompleteOrthogonalDecomposition<MatrixType_, StorageIndex_> >
     : traits<MatrixType_> {
   typedef MatrixXpr XprKind;
   typedef SolverStorage StorageKind;
-  typedef int StorageIndex;
+  typedef StorageIndex_ StorageIndex;
   enum { Flags = 0 };
 };
 
@@ -49,8 +49,8 @@
   * 
   * \sa MatrixBase::completeOrthogonalDecomposition()
   */
-template <typename MatrixType_> class CompleteOrthogonalDecomposition
-          : public SolverBase<CompleteOrthogonalDecomposition<MatrixType_> >
+template <typename MatrixType_, typename StorageIndex_> class CompleteOrthogonalDecomposition
+          : public SolverBase<CompleteOrthogonalDecomposition<MatrixType_, StorageIndex_> >
 {
  public:
   typedef MatrixType_ MatrixType;
@@ -65,7 +65,7 @@
     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
   };
   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
-  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime>
+  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, StorageIndex>
       PermutationType;
   typedef typename internal::plain_row_type<MatrixType, Index>::type
       IntRowVectorType;
@@ -78,9 +78,6 @@
       HouseholderSequenceType;
   typedef typename MatrixType::PlainObject PlainObject;
 
- private:
-  typedef typename PermutationType::Index PermIndexType;
-
  public:
   /**
    * \brief Default Constructor.
@@ -417,26 +414,26 @@
   template <typename Rhs>
   void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
 
-  ColPivHouseholderQR<MatrixType> m_cpqr;
+  ColPivHouseholderQR<MatrixType, StorageIndex> m_cpqr;
   HCoeffsType m_zCoeffs;
   RowVectorType m_temp;
 };
 
-template <typename MatrixType>
+template <typename MatrixType, typename StorageIndex>
 typename MatrixType::Scalar
-CompleteOrthogonalDecomposition<MatrixType>::determinant() const {
+CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::determinant() const {
   return m_cpqr.determinant();
 }
 
-template <typename MatrixType>
+template <typename MatrixType, typename StorageIndex>
 typename MatrixType::RealScalar
-CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const {
+CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::absDeterminant() const {
   return m_cpqr.absDeterminant();
 }
 
-template <typename MatrixType>
+template <typename MatrixType, typename StorageIndex>
 typename MatrixType::RealScalar
-CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const {
+CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::logAbsDeterminant() const {
   return m_cpqr.logAbsDeterminant();
 }
 
@@ -447,11 +444,10 @@
  * \sa class CompleteOrthogonalDecomposition,
  * CompleteOrthogonalDecomposition(const MatrixType&)
  */
-template <typename MatrixType>
-void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
+template <typename MatrixType, typename StorageIndex>
+void CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::computeInPlace()
 {
-  // the column permutation is stored as int indices, so just to be sure:
-  eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
+  eigen_assert(m_cpqr.cols() <= NumTraits<StorageIndex>::highest());
 
   const Index rank = m_cpqr.rank();
   const Index cols = m_cpqr.cols();
@@ -503,9 +499,9 @@
   }
 }
 
-template <typename MatrixType>
+template <typename MatrixType, typename StorageIndex>
 template <bool Conjugate, typename Rhs>
-void CompleteOrthogonalDecomposition<MatrixType>::applyZOnTheLeftInPlace(
+void CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::applyZOnTheLeftInPlace(
     Rhs& rhs) const {
   const Index cols = this->cols();
   const Index nrhs = rhs.cols();
@@ -525,9 +521,9 @@
   }
 }
 
-template <typename MatrixType>
+template <typename MatrixType, typename StorageIndex>
 template <typename Rhs>
-void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace(
+void CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::applyZAdjointOnTheLeftInPlace(
     Rhs& rhs) const {
   const Index cols = this->cols();
   const Index nrhs = rhs.cols();
@@ -548,9 +544,9 @@
 }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-template <typename MatrixType_>
+template <typename MatrixType_, typename StorageIndex_>
 template <typename RhsType, typename DstType>
-void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl(
+void CompleteOrthogonalDecomposition<MatrixType_, StorageIndex_>::_solve_impl(
     const RhsType& rhs, DstType& dst) const {
   const Index rank = this->rank();
   if (rank == 0) {
@@ -580,9 +576,9 @@
   dst = colsPermutation() * dst;
 }
 
-template<typename MatrixType_>
+template<typename MatrixType_, typename StorageIndex_>
 template<bool Conjugate, typename RhsType, typename DstType>
-void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
+void CompleteOrthogonalDecomposition<MatrixType_, StorageIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
 {
   const Index rank = this->rank();
 
@@ -611,17 +607,17 @@
 
 namespace internal {
 
-template<typename MatrixType>
-struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
+template<typename MatrixType, typename StorageIndex>
+struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType, StorageIndex> > >
   : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
 {
   enum { Flags = 0 };
 };
 
-template<typename DstXprType, typename MatrixType>
-struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
+template<typename DstXprType, typename MatrixType, typename StorageIndex>
+struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType, StorageIndex> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::Scalar>, Dense2Dense>
 {
-  typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
+  typedef CompleteOrthogonalDecomposition<MatrixType, StorageIndex> CodType;
   typedef Inverse<CodType> SrcXprType;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
   {
@@ -633,9 +629,9 @@
 } // end namespace internal
 
 /** \returns the matrix Q as a sequence of householder transformations */
-template <typename MatrixType>
-typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
-CompleteOrthogonalDecomposition<MatrixType>::householderQ() const {
+template <typename MatrixType, typename StorageIndex>
+typename CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::HouseholderSequenceType
+CompleteOrthogonalDecomposition<MatrixType, StorageIndex>::householderQ() const {
   return m_cpqr.householderQ();
 }
 
@@ -644,7 +640,8 @@
   * \sa class CompleteOrthogonalDecomposition
   */
 template <typename Derived>
-const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject>
+template <typename StorageIndex>
+const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject, StorageIndex>
 MatrixBase<Derived>::completeOrthogonalDecomposition() const {
   return CompleteOrthogonalDecomposition<PlainObject>(eval());
 }
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index ec7e19b..b00fd7b 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -17,19 +17,19 @@
 
 namespace internal {
 
-template<typename MatrixType_> struct traits<FullPivHouseholderQR<MatrixType_> >
+template<typename MatrixType_, typename StorageIndex_> struct traits<FullPivHouseholderQR<MatrixType_, StorageIndex_> >
  : traits<MatrixType_>
 {
   typedef MatrixXpr XprKind;
   typedef SolverStorage StorageKind;
-  typedef int StorageIndex;
+  typedef StorageIndex_ StorageIndex;
   enum { Flags = 0 };
 };
 
-template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
+template<typename MatrixType, typename StorageIndex> struct FullPivHouseholderQRMatrixQReturnType;
 
-template<typename MatrixType>
-struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
+template<typename MatrixType, typename StorageIndex>
+struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType, StorageIndex> >
 {
   typedef typename MatrixType::PlainObject ReturnType;
 };
@@ -59,8 +59,8 @@
   * 
   * \sa MatrixBase::fullPivHouseholderQr()
   */
-template<typename MatrixType_> class FullPivHouseholderQR
-        : public SolverBase<FullPivHouseholderQR<MatrixType_> >
+template<typename MatrixType_, typename StorageIndex_> class FullPivHouseholderQR
+        : public SolverBase<FullPivHouseholderQR<MatrixType_, StorageIndex_> >
 {
   public:
 
@@ -73,12 +73,12 @@
       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
-    typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
+    typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType, StorageIndex_> MatrixQReturnType;
     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
     typedef Matrix<StorageIndex, 1,
                    internal::min_size_prefer_dynamic(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
                    internal::min_size_prefer_fixed(MaxColsAtCompileTime, MaxRowsAtCompileTime)> IntDiagSizeVectorType;
-    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
+    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime, StorageIndex> PermutationType;
     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
     typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
     typedef typename MatrixType::PlainObject PlainObject;
@@ -437,8 +437,8 @@
     Index m_det_p;
 };
 
-template<typename MatrixType>
-typename MatrixType::Scalar FullPivHouseholderQR<MatrixType>::determinant() const
+template<typename MatrixType, typename StorageIndex>
+typename MatrixType::Scalar FullPivHouseholderQR<MatrixType, StorageIndex>::determinant() const
 {
   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@@ -447,8 +447,8 @@
   return m_qr.diagonal().prod() * detQ * Scalar(m_det_p);
 }
 
-template<typename MatrixType>
-typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
+template<typename MatrixType, typename StorageIndex>
+typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, StorageIndex>::absDeterminant() const
 {
   using std::abs;
   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
@@ -456,8 +456,8 @@
   return abs(m_qr.diagonal().prod());
 }
 
-template<typename MatrixType>
-typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
+template<typename MatrixType, typename StorageIndex>
+typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType, StorageIndex>::logAbsDeterminant() const
 {
   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
@@ -470,18 +470,19 @@
   *
   * \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&)
   */
-template<typename MatrixType>
+template<typename MatrixType, typename StorageIndex>
 template<typename InputType>
-FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
+FullPivHouseholderQR<MatrixType, StorageIndex>& FullPivHouseholderQR<MatrixType, StorageIndex>::compute(const EigenBase<InputType>& matrix)
 {
   m_qr = matrix.derived();
   computeInPlace();
   return *this;
 }
 
-template<typename MatrixType>
-void FullPivHouseholderQR<MatrixType>::computeInPlace()
+template<typename MatrixType, typename StorageIndex>
+void FullPivHouseholderQR<MatrixType, StorageIndex>::computeInPlace()
 {
+  eigen_assert(m_qr.cols() <= NumTraits<StorageIndex>::highest());
   using std::abs;
   Index rows = m_qr.rows();
   Index cols = m_qr.cols();
@@ -561,9 +562,9 @@
 }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-template<typename MatrixType_>
+template<typename MatrixType_, typename StorageIndex_>
 template<typename RhsType, typename DstType>
-void FullPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
+void FullPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl(const RhsType &rhs, DstType &dst) const
 {
   const Index l_rank = rank();
 
@@ -595,9 +596,9 @@
   for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
 }
 
-template<typename MatrixType_>
+template<typename MatrixType_, typename StorageIndex_>
 template<bool Conjugate, typename RhsType, typename DstType>
-void FullPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
+void FullPivHouseholderQR<MatrixType_, StorageIndex_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
 {
   const Index l_rank = rank();
 
@@ -634,10 +635,10 @@
 
 namespace internal {
   
-template<typename DstXprType, typename MatrixType>
-struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
+template<typename DstXprType, typename MatrixType, typename StorageIndex>
+struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType, StorageIndex> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType, StorageIndex>::Scalar>, Dense2Dense>
 {
-  typedef FullPivHouseholderQR<MatrixType> QrType;
+  typedef FullPivHouseholderQR<MatrixType, StorageIndex> QrType;
   typedef Inverse<QrType> SrcXprType;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
   {    
@@ -651,11 +652,11 @@
   *
   * \tparam MatrixType type of underlying dense matrix
   */
-template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
-  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
+template<typename MatrixType, typename StorageIndex> struct FullPivHouseholderQRMatrixQReturnType
+  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType, StorageIndex> >
 {
 public:
-  typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
+  typedef typename FullPivHouseholderQR<MatrixType, StorageIndex>::IntDiagSizeVectorType IntDiagSizeVectorType;
   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
   typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
                  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
@@ -712,8 +713,8 @@
 
 } // end namespace internal
 
-template<typename MatrixType>
-inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
+template<typename MatrixType, typename StorageIndex>
+inline typename FullPivHouseholderQR<MatrixType, StorageIndex>::MatrixQReturnType FullPivHouseholderQR<MatrixType, StorageIndex>::matrixQ() const
 {
   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
   return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
@@ -724,10 +725,11 @@
   * \sa class FullPivHouseholderQR
   */
 template<typename Derived>
-const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
+template<typename StorageIndex>
+const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject, StorageIndex>
 MatrixBase<Derived>::fullPivHouseholderQr() const
 {
-  return FullPivHouseholderQR<PlainObject>(eval());
+  return FullPivHouseholderQR<PlainObject, StorageIndex>(eval());
 }
 
 } // end namespace Eigen
diff --git a/Eigen/src/misc/lapacke_helpers.h b/Eigen/src/misc/lapacke_helpers.h
index b6ad6e8..8c190d0 100644
--- a/Eigen/src/misc/lapacke_helpers.h
+++ b/Eigen/src/misc/lapacke_helpers.h
@@ -77,9 +77,9 @@
 }
 
 /// translate UpLo type to the corresponding letter code
-template<UpLoType mode> char translate_mode;
-template<> constexpr char translate_mode<Lower> = 'L';
-template<> constexpr char translate_mode<Upper> = 'U';
+//template<UpLoType mode> char translate_mode;
+//template<> constexpr char translate_mode<Lower> = 'L';
+//template<> constexpr char translate_mode<Upper> = 'U';
 
 
 // ---------------------------------------------------------------------------------------------------------------------
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index 4185f51..9f5e3dc 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -13,9 +13,8 @@
 #include <Eigen/SVD>
 #include "solverbase.h"
 
-template <typename MatrixType>
+template <typename MatrixType, typename StorageIndex>
 void cod() {
-  STATIC_CHECK(( internal::is_same<typename CompleteOrthogonalDecomposition<MatrixType>::StorageIndex,int>::value ));
 
   Index rows = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
   Index cols = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
@@ -28,7 +27,7 @@
       MatrixQType;
   MatrixType matrix;
   createRandomPIMatrixOfRank(rank, rows, cols, matrix);
-  CompleteOrthogonalDecomposition<MatrixType> cod(matrix);
+  CompleteOrthogonalDecomposition<MatrixType, StorageIndex> cod(matrix);
   VERIFY(rank == cod.rank());
   VERIFY(cols - cod.rank() == cod.dimensionOfKernel());
   VERIFY(!cod.isInjective());
@@ -63,14 +62,14 @@
   VERIFY_IS_APPROX(cod_solution, pinv * rhs);
 }
 
-template <typename MatrixType, int Cols2>
+template <typename MatrixType, typename StorageIndex, int Cols2>
 void cod_fixedsize() {
   enum {
     Rows = MatrixType::RowsAtCompileTime,
     Cols = MatrixType::ColsAtCompileTime
   };
   typedef typename MatrixType::Scalar Scalar;
-  typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > COD;
+  typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols>, StorageIndex> COD;
   int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1);
   Matrix<Scalar, Rows, Cols> matrix;
   createRandomPIMatrixOfRank(rank, Rows, Cols, matrix);
@@ -96,12 +95,10 @@
   VERIFY_IS_APPROX(cod_solution, pinv * rhs);
 }
 
-template<typename MatrixType> void qr()
+template<typename MatrixType, typename StorageIndex> void qr()
 {
   using std::sqrt;
 
-  STATIC_CHECK(( internal::is_same<typename ColPivHouseholderQR<MatrixType>::StorageIndex,int>::value ));
-
   Index rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols2 = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE);
   Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1);
 
@@ -110,7 +107,7 @@
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
   MatrixType m1;
   createRandomPIMatrixOfRank(rank,rows,cols,m1);
-  ColPivHouseholderQR<MatrixType> qr(m1);
+  ColPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
   VERIFY_IS_EQUAL(rank, qr.rank());
   VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel());
   VERIFY(!qr.isInjective());
@@ -158,7 +155,7 @@
   }
 }
 
-template<typename MatrixType, int Cols2> void qr_fixedsize()
+template<typename MatrixType, typename StorageIndex, int Cols2> void qr_fixedsize()
 {
   using std::sqrt;
   using std::abs;
@@ -168,7 +165,7 @@
   int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols))-1);
   Matrix<Scalar,Rows,Cols> m1;
   createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
-  ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
+  ColPivHouseholderQR<Matrix<Scalar,Rows,Cols>, StorageIndex> qr(m1);
   VERIFY_IS_EQUAL(rank, qr.rank());
   VERIFY_IS_EQUAL(Cols - qr.rank(), qr.dimensionOfKernel());
   VERIFY_IS_EQUAL(qr.isInjective(), (rank == Rows));
@@ -207,7 +204,7 @@
 // for rank-revealing QR. See
 // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
 // page 3 for more detail.
-template<typename MatrixType> void qr_kahan_matrix()
+template<typename MatrixType, typename StorageIndex> void qr_kahan_matrix()
 {
   using std::sqrt;
   using std::abs;
@@ -227,7 +224,7 @@
     pow_s_i *= s;
   }
   m1 = (m1 + m1.transpose()).eval();
-  ColPivHouseholderQR<MatrixType> qr(m1);
+  ColPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
   MatrixType r = qr.matrixQR().template triangularView<Upper>();
 
   RealScalar threshold =
@@ -247,7 +244,7 @@
   }
 }
 
-template<typename MatrixType> void qr_invertible()
+template<typename MatrixType, typename StorageIndex> void qr_invertible()
 {
   using std::log;
   using std::abs;
@@ -266,7 +263,7 @@
     m1 += a * a.adjoint();
   }
 
-  ColPivHouseholderQR<MatrixType> qr(m1);
+  ColPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
 
   check_solverbase<MatrixType, MatrixType>(m1, qr, size, size, size);
 
@@ -283,11 +280,11 @@
   VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
 }
 
-template<typename MatrixType> void qr_verify_assert()
+template<typename MatrixType, typename StorageIndex> void qr_verify_assert()
 {
   MatrixType tmp;
 
-  ColPivHouseholderQR<MatrixType> qr;
+  ColPivHouseholderQR<MatrixType, StorageIndex> qr;
   VERIFY_RAISES_ASSERT(qr.matrixQR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp))
   VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp))
@@ -303,11 +300,11 @@
   VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
 }
 
-template<typename MatrixType> void cod_verify_assert()
+template<typename MatrixType, typename StorageIndex> void cod_verify_assert()
 {
   MatrixType tmp;
 
-  CompleteOrthogonalDecomposition<MatrixType> cod;
+  CompleteOrthogonalDecomposition<MatrixType, StorageIndex> cod;
   VERIFY_RAISES_ASSERT(cod.matrixQTZ())
   VERIFY_RAISES_ASSERT(cod.solve(tmp))
   VERIFY_RAISES_ASSERT(cod.transpose().solve(tmp))
@@ -323,50 +320,58 @@
   VERIFY_RAISES_ASSERT(cod.logAbsDeterminant())
 }
 
+
+
 EIGEN_DECLARE_TEST(qr_colpivoting)
 {
+  #if defined(EIGEN_USE_LAPACKE)
+  typedef lapack_int StorageIndex;
+  #else
+  typedef int StorageIndex;
+  #endif
+
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST_1( qr<MatrixXf>() );
-    CALL_SUBTEST_2( qr<MatrixXd>() );
-    CALL_SUBTEST_3( qr<MatrixXcd>() );
-    CALL_SUBTEST_4(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
-    CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
-    CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,1,1>, 1 >() ));
+    CALL_SUBTEST_1( (qr<MatrixXf, StorageIndex>)() );
+    CALL_SUBTEST_2( (qr<MatrixXf, StorageIndex>)() );
+    CALL_SUBTEST_3( (qr<MatrixXf, StorageIndex>)() );
+    CALL_SUBTEST_4(( (qr_fixedsize<Matrix<float,3,5>, StorageIndex, 4>)() ));
+    CALL_SUBTEST_5(( (qr_fixedsize<Matrix<double,6,2>, StorageIndex, 3>)() ));
+    CALL_SUBTEST_5(( (qr_fixedsize<Matrix<double,1,1>, StorageIndex, 1>)() ));
   }
 
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST_1( cod<MatrixXf>() );
-    CALL_SUBTEST_2( cod<MatrixXd>() );
-    CALL_SUBTEST_3( cod<MatrixXcd>() );
-    CALL_SUBTEST_4(( cod_fixedsize<Matrix<float,3,5>, 4 >() ));
-    CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,6,2>, 3 >() ));
-    CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,1,1>, 1 >() ));
+    CALL_SUBTEST_1( (cod<MatrixXf, StorageIndex>)() );
+    CALL_SUBTEST_2( (cod<MatrixXf, StorageIndex>)() );
+    CALL_SUBTEST_3( (cod<MatrixXf, StorageIndex>)() );
+    CALL_SUBTEST_4(( (cod_fixedsize<Matrix<float,3,5>, StorageIndex, 4>)() ));
+    CALL_SUBTEST_5(( (cod_fixedsize<Matrix<double,6,2>, StorageIndex, 3>)() ));
+    CALL_SUBTEST_5(( (cod_fixedsize<Matrix<double,1,1>, StorageIndex, 1>)() ));
   }
 
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST_1( qr_invertible<MatrixXf>() );
-    CALL_SUBTEST_2( qr_invertible<MatrixXd>() );
-    CALL_SUBTEST_6( qr_invertible<MatrixXcf>() );
-    CALL_SUBTEST_3( qr_invertible<MatrixXcd>() );
+    CALL_SUBTEST_1( (qr_invertible<MatrixXf, StorageIndex>)() );
+    CALL_SUBTEST_2( (qr_invertible<MatrixXd, StorageIndex>)() );
+    CALL_SUBTEST_6( (qr_invertible<MatrixXcf, StorageIndex>)() );
+    CALL_SUBTEST_3( (qr_invertible<MatrixXcd, StorageIndex>)() );
   }
 
-  CALL_SUBTEST_7(qr_verify_assert<Matrix3f>());
-  CALL_SUBTEST_8(qr_verify_assert<Matrix3d>());
-  CALL_SUBTEST_1(qr_verify_assert<MatrixXf>());
-  CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
-  CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>());
-  CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
+  CALL_SUBTEST_7( (qr_verify_assert<Matrix3f, StorageIndex>) ());
+  CALL_SUBTEST_8( (qr_verify_assert<Matrix3d, StorageIndex>)());
+  CALL_SUBTEST_1( (qr_verify_assert<MatrixXf, StorageIndex>)());
+  CALL_SUBTEST_2( (qr_verify_assert<MatrixXd, StorageIndex>)());
+  CALL_SUBTEST_6( (qr_verify_assert<MatrixXcf, StorageIndex>)());
+  CALL_SUBTEST_3( (qr_verify_assert<MatrixXcd, StorageIndex>)());
 
-  CALL_SUBTEST_7(cod_verify_assert<Matrix3f>());
-  CALL_SUBTEST_8(cod_verify_assert<Matrix3d>());
-  CALL_SUBTEST_1(cod_verify_assert<MatrixXf>());
-  CALL_SUBTEST_2(cod_verify_assert<MatrixXd>());
-  CALL_SUBTEST_6(cod_verify_assert<MatrixXcf>());
-  CALL_SUBTEST_3(cod_verify_assert<MatrixXcd>());
+  CALL_SUBTEST_7( (cod_verify_assert<Matrix3f, StorageIndex>)());
+  CALL_SUBTEST_8( (cod_verify_assert<Matrix3d, StorageIndex>)());
+  CALL_SUBTEST_1( (cod_verify_assert<MatrixXf, StorageIndex>)());
+  CALL_SUBTEST_2( (cod_verify_assert<MatrixXd, StorageIndex>)());
+  CALL_SUBTEST_6( (cod_verify_assert<MatrixXcf, StorageIndex>)());
+  CALL_SUBTEST_3( (cod_verify_assert<MatrixXcd, StorageIndex>)());
 
   // Test problem size constructors
-  CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20));
+  CALL_SUBTEST_9((ColPivHouseholderQR<MatrixXf, StorageIndex>(10, 20)));
 
-  CALL_SUBTEST_1( qr_kahan_matrix<MatrixXf>() );
-  CALL_SUBTEST_2( qr_kahan_matrix<MatrixXd>() );
+  CALL_SUBTEST_1( (qr_kahan_matrix<MatrixXf, StorageIndex>)() );
+  CALL_SUBTEST_2( (qr_kahan_matrix<MatrixXd, StorageIndex>)() );
 }
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp
index cca9a8c..971ec73 100644
--- a/test/qr_fullpivoting.cpp
+++ b/test/qr_fullpivoting.cpp
@@ -12,9 +12,8 @@
 #include <Eigen/QR>
 #include "solverbase.h"
 
-template<typename MatrixType> void qr()
+template<typename MatrixType, typename StorageIndex> void qr()
 {
-  STATIC_CHECK(( internal::is_same<typename FullPivHouseholderQR<MatrixType>::StorageIndex,int>::value ));
 
   static const int Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime;
   Index max_size = EIGEN_TEST_MAX_SIZE;
@@ -28,7 +27,7 @@
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
   MatrixType m1;
   createRandomPIMatrixOfRank(rank,rows,cols,m1);
-  FullPivHouseholderQR<MatrixType> qr(m1);
+  FullPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
   VERIFY_IS_EQUAL(rank, qr.rank());
   VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel());
   VERIFY(!qr.isInjective());
@@ -67,7 +66,7 @@
   }
 }
 
-template<typename MatrixType> void qr_invertible()
+template<typename MatrixType, typename StorageIndex> void qr_invertible()
 {
   using std::log;
   using std::abs;
@@ -88,7 +87,7 @@
     m1 += a * a.adjoint();
   }
 
-  FullPivHouseholderQR<MatrixType> qr(m1);
+  FullPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
   VERIFY(qr.isInjective());
   VERIFY(qr.isInvertible());
   VERIFY(qr.isSurjective());
@@ -108,11 +107,11 @@
   VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
 }
 
-template<typename MatrixType> void qr_verify_assert()
+template<typename MatrixType, typename StorageIndex> void qr_verify_assert()
 {
   MatrixType tmp;
 
-  FullPivHouseholderQR<MatrixType> qr;
+  FullPivHouseholderQR<MatrixType, StorageIndex> qr;
   VERIFY_RAISES_ASSERT(qr.matrixQR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp))
   VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp))
@@ -130,33 +129,35 @@
 
 EIGEN_DECLARE_TEST(qr_fullpivoting)
 {
+  typedef int StorageIndex;
+
   for(int i = 0; i < 1; i++) {
-    CALL_SUBTEST_5( qr<Matrix3f>() );
-    CALL_SUBTEST_6( qr<Matrix3d>() );
-    CALL_SUBTEST_8( qr<Matrix2f>() );
-    CALL_SUBTEST_1( qr<MatrixXf>() );
-    CALL_SUBTEST_2( qr<MatrixXd>() );
-    CALL_SUBTEST_3( qr<MatrixXcd>() );
+    CALL_SUBTEST_5( (qr<Matrix3f, StorageIndex>()) );
+    CALL_SUBTEST_6( (qr<Matrix3d, StorageIndex>()) );
+    CALL_SUBTEST_8( (qr<Matrix2f, StorageIndex>()) );
+    CALL_SUBTEST_1( (qr<MatrixXf, StorageIndex>()) );
+    CALL_SUBTEST_2( (qr<MatrixXd, StorageIndex>()) );
+    CALL_SUBTEST_3( (qr<MatrixXcd, StorageIndex>()) );
   }
 
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST_1( qr_invertible<MatrixXf>() );
-    CALL_SUBTEST_2( qr_invertible<MatrixXd>() );
-    CALL_SUBTEST_4( qr_invertible<MatrixXcf>() );
-    CALL_SUBTEST_3( qr_invertible<MatrixXcd>() );
+    CALL_SUBTEST_1( (qr_invertible<MatrixXf, StorageIndex>()) );
+    CALL_SUBTEST_2( (qr_invertible<MatrixXd, StorageIndex>()) );
+    CALL_SUBTEST_4( (qr_invertible<MatrixXcf, StorageIndex>()) );
+    CALL_SUBTEST_3( (qr_invertible<MatrixXcd, StorageIndex>()) );
   }
 
-  CALL_SUBTEST_5(qr_verify_assert<Matrix3f>());
-  CALL_SUBTEST_6(qr_verify_assert<Matrix3d>());
-  CALL_SUBTEST_1(qr_verify_assert<MatrixXf>());
-  CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
-  CALL_SUBTEST_4(qr_verify_assert<MatrixXcf>());
-  CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
+  CALL_SUBTEST_5( (qr_verify_assert<Matrix3f, StorageIndex>()) );
+  CALL_SUBTEST_6( (qr_verify_assert<Matrix3d, StorageIndex>()) );
+  CALL_SUBTEST_1( (qr_verify_assert<MatrixXf, StorageIndex>()) );
+  CALL_SUBTEST_2( (qr_verify_assert<MatrixXd, StorageIndex>()) );
+  CALL_SUBTEST_4( (qr_verify_assert<MatrixXcf, StorageIndex>()) );
+  CALL_SUBTEST_3( (qr_verify_assert<MatrixXcd, StorageIndex>()) );
 
   // Test problem size constructors
-  CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20));
-  CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(10,20)));
-  CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(Matrix<float,10,20>::Random())));
-  CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(20,10)));
-  CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(Matrix<float,20,10>::Random())));
+  CALL_SUBTEST_7( (FullPivHouseholderQR<MatrixXf,StorageIndex>(10, 20)));
+  CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 10, 20>, StorageIndex>(10, 20)));
+  CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 10, 20>, StorageIndex>(Matrix<float,10,20>::Random())));
+  CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 20, 10>, StorageIndex>(20, 10)));
+  CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 20, 10>, StorageIndex>(Matrix<float,20,10>::Random())));
 }