diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 607fe261f..28dc0cd47 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -594,7 +594,8 @@ struct ei_image_retval<FullPivLU<_MatrixType> >
         pivots.coeffRef(p++) = i;
     ei_internal_assert(p == rank());
 
-    dst.block(0,0,dst.rows(),rank()) = (originalMatrix()*dec().permutationQ()).block(0,0,dst.rows(),rank());
+    for(int i = 0; i < rank(); ++i)
+      dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
   }
 };
 
@@ -630,8 +631,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
     typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols());
 
     // Step 1
-    for(int i = 0; i < rows; ++i)
-      c.row(dec().permutationP().indices().coeff(i)) = rhs().row(i);
+    c = dec().permutationP() * rhs();
 
     // Step 2
     dec().matrixLU()
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index d3c6b4c7c..e764af529 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -57,10 +57,9 @@ template<typename _MatrixType> class ColPivHouseholderQR
     typedef typename MatrixType::RealScalar RealScalar;
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
     typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
+    typedef PermutationMatrix<ColsAtCompileTime> PermutationType;
     typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
-    typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
     typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
     typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
 
@@ -117,7 +116,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
 
     ColPivHouseholderQR& compute(const MatrixType& matrix);
 
-    const IntRowVectorType& colsPermutation() const
+    const PermutationType& colsPermutation() const
     {
       ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_cols_permutation;
@@ -230,7 +229,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
-    IntRowVectorType m_cols_permutation;
+    PermutationType m_cols_permutation;
     bool m_isInitialized;
     RealScalar m_precision;
     int m_rank;
@@ -271,7 +270,6 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   m_precision = epsilon<Scalar>() * size;
 
   IntRowVectorType cols_transpositions(matrix.cols());
-  m_cols_permutation.resize(matrix.cols());
   int number_of_transpositions = 0;
 
   RealRowVectorType colSqNorms(cols);
@@ -314,9 +312,9 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
     colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2();
   }
 
-  for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k;
+  m_cols_permutation.setIdentity(cols);
   for(int k = 0; k < size; ++k)
-    std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k)));
+    m_cols_permutation.applyTranspositionOnTheLeft(k, cols_transpositions.coeff(k));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;
@@ -368,8 +366,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
        .template triangularView<UpperTriangular>()
        .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols()));
 
-    for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i);
-    for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero();
+    dst = dec().colsPermutation() * c;
   }
 };
 
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index 07df5ed6b..ac0108046 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -58,6 +58,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
     typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
     typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
+    typedef PermutationMatrix<ColsAtCompileTime> PermutationType;
     typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
@@ -112,7 +113,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
 
     FullPivHouseholderQR& compute(const MatrixType& matrix);
 
-    const IntRowVectorType& colsPermutation() const
+    const PermutationType& colsPermutation() const
     {
       ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
       return m_cols_permutation;
@@ -231,7 +232,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
     IntColVectorType m_rows_transpositions;
-    IntRowVectorType m_cols_permutation;
+    PermutationType m_cols_permutation;
     bool m_isInitialized;
     RealScalar m_precision;
     int m_rank;
@@ -273,7 +274,6 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
 
   m_rows_transpositions.resize(matrix.rows());
   IntRowVectorType cols_transpositions(matrix.cols());
-  m_cols_permutation.resize(matrix.cols());
   int number_of_transpositions = 0;
 
   RealScalar biggest(0);
@@ -322,9 +322,9 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
         .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
   }
 
-  for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k;
+  m_cols_permutation.setIdentity(cols);
   for(int k = 0; k < size; ++k)
-    std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k)));
+    m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;
@@ -379,8 +379,8 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
        .template triangularView<UpperTriangular>()
        .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols()));
 
-    for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i);
-    for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero();
+    for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
+    for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
   }
 };
 
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp
index 65d9a071f..c82ba4c7e 100644
--- a/test/qr_fullpivoting.cpp
+++ b/test/qr_fullpivoting.cpp
@@ -51,11 +51,8 @@ template<typename MatrixType> void qr()
   // FIXME need better way to construct trapezoid
   for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
 
-  MatrixType b = qr.matrixQ() * r;
+  MatrixType c = qr.matrixQ() * r * qr.colsPermutation().inverse();
 
-  MatrixType c = MatrixType::Zero(rows,cols);
-
-  for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i);
   VERIFY_IS_APPROX(m1, c);
 
   MatrixType m2 = MatrixType::Random(cols,cols2);