diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h
index af9a1fe..7e402cc 100644
--- a/Eigen/src/SparseCore/SparsePermutation.h
+++ b/Eigen/src/SparseCore/SparsePermutation.h
@@ -18,69 +18,142 @@
 
 namespace internal {
 
+template<typename ExpressionType, typename PlainObjectType, bool NeedEval = !is_same<ExpressionType, PlainObjectType>::value>
+struct XprHelper
+{
+    XprHelper(const ExpressionType& xpr) : m_xpr(xpr) {}
+    inline const PlainObjectType& xpr() const { return m_xpr; }
+    // this is a new PlainObjectType initialized by xpr
+    const PlainObjectType m_xpr;
+};
+template<typename ExpressionType, typename PlainObjectType>
+struct XprHelper<ExpressionType, PlainObjectType, false>
+{
+    XprHelper(const ExpressionType& xpr) : m_xpr(xpr) {}
+    inline const PlainObjectType& xpr() const { return m_xpr; }
+    // this is a reference to xpr
+    const PlainObjectType& m_xpr;
+};
+
+template<typename PermDerived, bool NeedInverseEval>
+struct PermHelper
+{
+    using IndicesType = typename PermDerived::IndicesType;
+    using PermutationIndex = typename IndicesType::Scalar;
+    using type = PermutationMatrix<IndicesType::SizeAtCompileTime, IndicesType::MaxSizeAtCompileTime, PermutationIndex>;
+    PermHelper(const PermDerived& perm) : m_perm(perm.inverse()) {}
+    inline const type& perm() const { return m_perm; }
+    // this is a new PermutationMatrix initialized by perm.inverse()
+    const type m_perm;
+};
+template<typename PermDerived>
+struct PermHelper<PermDerived, false>
+{
+    using type = PermDerived;
+    PermHelper(const PermDerived& perm) : m_perm(perm) {}
+    inline const type& perm() const { return m_perm; }
+    // this is a reference to perm
+    const type& m_perm;
+};
+
 template<typename ExpressionType, int Side, bool Transposed>
 struct permutation_matrix_product<ExpressionType, Side, Transposed, SparseShape>
 {
-    typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
-    typedef remove_all_t<MatrixType> MatrixTypeCleaned;
+    using MatrixType = typename nested_eval<ExpressionType, 1>::type;
+    using MatrixTypeCleaned = remove_all_t<MatrixType>;
 
-    typedef typename MatrixTypeCleaned::Scalar Scalar;
-    typedef typename MatrixTypeCleaned::StorageIndex StorageIndex;
+    using Scalar = typename MatrixTypeCleaned::Scalar;
+    using StorageIndex = typename MatrixTypeCleaned::StorageIndex;
 
-    enum {
-      SrcStorageOrder = MatrixTypeCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
-      MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
-    };
-    
-    typedef std::conditional_t<MoveOuter,
-        SparseMatrix<Scalar,SrcStorageOrder,StorageIndex>,
-        SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> > ReturnType;
+    // the actual "return type" is `Dest`. this is a temporary type
+    using ReturnType = SparseMatrix<Scalar, MatrixTypeCleaned::IsRowMajor ? RowMajor : ColMajor, StorageIndex>;
+    using TmpHelper = XprHelper<ExpressionType, ReturnType>;
 
-    template<typename Dest,typename PermutationType>
-    static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr)
-    {
-      MatrixType mat(xpr);
-      if(MoveOuter)
-      {
-        SparseMatrix<Scalar,SrcStorageOrder,StorageIndex> tmp(mat.rows(), mat.cols());
-        Matrix<StorageIndex,Dynamic,1> sizes(mat.outerSize());
-        for(Index j=0; j<mat.outerSize(); ++j)
-        {
+    static constexpr bool NeedOuterPermutation = ExpressionType::IsRowMajor ? Side == OnTheLeft : Side == OnTheRight;
+    static constexpr bool NeedInversePermutation = Transposed ? Side == OnTheLeft : Side == OnTheRight;
+
+    template <typename Dest, typename PermutationType>
+    static inline void permute_outer(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) {
+
+        // if ExpressionType is not ReturnType, evaluate `xpr` (allocation)
+        // otherwise, just reference `xpr`
+        // TODO: handle trivial expressions such as CwiseBinaryOp without temporary
+        const TmpHelper tmpHelper(xpr);
+        const ReturnType& tmp = tmpHelper.xpr();
+
+        ReturnType result(tmp.rows(), tmp.cols());
+
+        for (Index j = 0; j < tmp.outerSize(); j++) {
           Index jp = perm.indices().coeff(j);
-          sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = StorageIndex(mat.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros());
+          Index jsrc = NeedInversePermutation ? jp : j;
+          Index jdst = NeedInversePermutation ? j : jp;
+          Index begin = tmp.outerIndexPtr()[jsrc];
+          Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[jsrc + 1] : begin + tmp.innerNonZeroPtr()[jsrc];
+          result.outerIndexPtr()[jdst + 1] += end - begin;
         }
-        tmp.reserve(sizes);
-        for(Index j=0; j<mat.outerSize(); ++j)
-        {
-          Index jp = perm.indices().coeff(j);
-          Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j;
-          Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j;
-          for(typename MatrixTypeCleaned::InnerIterator it(mat,jsrc); it; ++it)
-            tmp.insertByOuterInner(jdst,it.index()) = it.value();
-        }
-        dst = tmp;
-      }
-      else
-      {
-        SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> tmp(mat.rows(), mat.cols());
-        Matrix<StorageIndex,Dynamic,1> sizes(tmp.outerSize());
-        sizes.setZero();
-        PermutationMatrix<Dynamic,Dynamic,StorageIndex> perm_cpy;
-        if((Side==OnTheLeft) ^ Transposed)
-          perm_cpy = perm;
-        else
-          perm_cpy = perm.transpose();
 
-        for(Index j=0; j<mat.outerSize(); ++j)
-          for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it)
-            sizes[perm_cpy.indices().coeff(it.index())]++;
-        tmp.reserve(sizes);
-        for(Index j=0; j<mat.outerSize(); ++j)
-          for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it)
-            tmp.insertByOuterInner(perm_cpy.indices().coeff(it.index()),j) = it.value();
-        dst = tmp;
-      }
+        std::partial_sum(result.outerIndexPtr(), result.outerIndexPtr() + result.outerSize() + 1,
+                         result.outerIndexPtr());
+        result.resizeNonZeros(result.nonZeros());
+
+        for (Index j = 0; j < tmp.outerSize(); j++) {
+          Index jp = perm.indices().coeff(j);
+          Index jsrc = NeedInversePermutation ? jp : j;
+          Index jdst = NeedInversePermutation ? j : jp;
+          Index begin = tmp.outerIndexPtr()[jsrc];
+          Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[jsrc + 1] : begin + tmp.innerNonZeroPtr()[jsrc];
+          Index target = result.outerIndexPtr()[jdst];
+          smart_copy(tmp.innerIndexPtr() + begin, tmp.innerIndexPtr() + end, result.innerIndexPtr() + target);
+          smart_copy(tmp.valuePtr() + begin, tmp.valuePtr() + end, result.valuePtr() + target);
+        }
+        dst = std::move(result);
     }
+
+    template <typename Dest, typename PermutationType>
+    static inline void permute_inner(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) {
+        using InnerPermHelper = PermHelper<PermutationType, NeedInversePermutation>;
+        using InnerPermType = typename InnerPermHelper::type;
+
+        // if ExpressionType is not ReturnType, evaluate `xpr` (allocation)
+        // otherwise, just reference `xpr`
+        // TODO: handle trivial expressions such as CwiseBinaryOp without temporary
+        const TmpHelper tmpHelper(xpr);
+        const ReturnType& tmp = tmpHelper.xpr();
+
+        // if inverse permutation of inner indices is requested, calculate perm.inverse() (allocation)
+        // otherwise, just reference `perm`
+        const InnerPermHelper permHelper(perm);
+        const InnerPermType& innerPerm = permHelper.perm();
+
+        ReturnType result(tmp.rows(), tmp.cols());
+
+        for (Index j = 0; j < tmp.outerSize(); j++) {
+            Index begin = tmp.outerIndexPtr()[j];
+            Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[j + 1] : begin + tmp.innerNonZeroPtr()[j];
+            result.outerIndexPtr()[j + 1] += end - begin;
+        }
+
+        std::partial_sum(result.outerIndexPtr(), result.outerIndexPtr() + result.outerSize() + 1, result.outerIndexPtr());
+        result.resizeNonZeros(result.nonZeros());
+
+        for (Index j = 0; j < tmp.outerSize(); j++) {
+            Index begin = tmp.outerIndexPtr()[j];
+            Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[j + 1] : begin + tmp.innerNonZeroPtr()[j];
+            Index target = result.outerIndexPtr()[j];
+            std::transform(tmp.innerIndexPtr() + begin, tmp.innerIndexPtr() + end, result.innerIndexPtr() + target,
+                           [&innerPerm](StorageIndex i) { return innerPerm.indices().coeff(i); });
+            smart_copy(tmp.valuePtr() + begin, tmp.valuePtr() + end, result.valuePtr() + target);
+        }
+        // the inner indices were permuted, and must be sorted
+        result.sortInnerIndices();
+        dst = std::move(result);
+    }
+
+    template <typename Dest, typename PermutationType, bool DoOuter = NeedOuterPermutation, std::enable_if_t<DoOuter, int> = 0>
+    static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_outer(dst, perm, xpr); }
+
+    template <typename Dest, typename PermutationType, bool DoOuter = NeedOuterPermutation, std::enable_if_t<!DoOuter, int> = 0>
+    static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_inner(dst, perm, xpr); }
 };
 
 }
@@ -143,35 +216,34 @@
 } // end namespace internal
 
 /** \returns the matrix with the permutation applied to the columns
-  */
-template<typename SparseDerived, typename PermDerived>
-inline const Product<SparseDerived, PermDerived, AliasFreeProduct>
-operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm)
-{ return Product<SparseDerived, PermDerived, AliasFreeProduct>(matrix.derived(), perm.derived()); }
+ */
+template <typename SparseDerived, typename PermDerived>
+inline const Product<SparseDerived, PermDerived, AliasFreeProduct> operator*(
+    const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm) {
+  return Product<SparseDerived, PermDerived, AliasFreeProduct>(matrix.derived(), perm.derived());
+}
 
 /** \returns the matrix with the permutation applied to the rows
-  */
-template<typename SparseDerived, typename PermDerived>
-inline const Product<PermDerived, SparseDerived, AliasFreeProduct>
-operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix)
-{ return  Product<PermDerived, SparseDerived, AliasFreeProduct>(perm.derived(), matrix.derived()); }
-
+ */
+template <typename SparseDerived, typename PermDerived>
+inline const Product<PermDerived, SparseDerived, AliasFreeProduct> operator*(
+    const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix) {
+  return Product<PermDerived, SparseDerived, AliasFreeProduct>(perm.derived(), matrix.derived());
+}
 
 /** \returns the matrix with the inverse permutation applied to the columns.
-  */
-template<typename SparseDerived, typename PermutationType>
-inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct>
-operator*(const SparseMatrixBase<SparseDerived>& matrix, const InverseImpl<PermutationType, PermutationStorage>& tperm)
-{
+ */
+template <typename SparseDerived, typename PermutationType>
+inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct> operator*(
+    const SparseMatrixBase<SparseDerived>& matrix, const InverseImpl<PermutationType, PermutationStorage>& tperm) {
   return Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct>(matrix.derived(), tperm.derived());
 }
 
 /** \returns the matrix with the inverse permutation applied to the rows.
-  */
-template<typename SparseDerived, typename PermutationType>
-inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>
-operator*(const InverseImpl<PermutationType,PermutationStorage>& tperm, const SparseMatrixBase<SparseDerived>& matrix)
-{
+ */
+template <typename SparseDerived, typename PermutationType>
+inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct> operator*(
+    const InverseImpl<PermutationType, PermutationStorage>& tperm, const SparseMatrixBase<SparseDerived>& matrix) {
   return Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>(tperm.derived(), matrix.derived());
 }
 
diff --git a/test/sparse_permutations.cpp b/test/sparse_permutations.cpp
index 5974c74..775c560 100644
--- a/test/sparse_permutations.cpp
+++ b/test/sparse_permutations.cpp
@@ -17,6 +17,15 @@
     VERIFY( (#XPR) && nb_transposed_copies==N ); \
   }
 
+static long int nb_temporaries;
+#define EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN {nb_temporaries++;}
+#define VERIFY_TEMPORARY_COUNT(XPR,N) {\
+    nb_temporaries = 0; \
+    XPR; \
+    if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \
+    VERIFY( (#XPR) && nb_temporaries==N ); \
+  }
+
 #include "sparse.h"
 
 template<typename T>
@@ -79,28 +88,52 @@
   VERIFY( is_sorted( ::eval(mat*p) ));
   VERIFY( is_sorted( res = mat*p ));
   VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p), 0);
-  //VERIFY_TRANSPOSITION_COUNT( res = mat*p, IsRowMajor ? 1 : 0 );
+  VERIFY_TEMPORARY_COUNT(     ::eval(mat*p), 1)
   res_d = mat_d*p;
   VERIFY(res.isApprox(res_d) && "mat*p");
 
   VERIFY( is_sorted( ::eval(p*mat) ));
   VERIFY( is_sorted( res = p*mat ));
   VERIFY_TRANSPOSITION_COUNT( ::eval(p*mat), 0);
+  VERIFY_TEMPORARY_COUNT(     ::eval(p*mat), 1);
   res_d = p*mat_d;
   VERIFY(res.isApprox(res_d) && "p*mat");
 
   VERIFY( is_sorted( (mat*p).eval() ));
   VERIFY( is_sorted( res = mat*p.inverse() ));
   VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p.inverse()), 0);
+  VERIFY_TEMPORARY_COUNT(     ::eval(mat*p.inverse()), 1);
   res_d = mat*p.inverse();
   VERIFY(res.isApprox(res_d) && "mat*inv(p)");
 
   VERIFY( is_sorted( (p*mat+p*mat).eval() ));
   VERIFY( is_sorted( res = p.inverse()*mat ));
   VERIFY_TRANSPOSITION_COUNT( ::eval(p.inverse()*mat), 0);
+  VERIFY_TEMPORARY_COUNT(     ::eval(p.inverse()*mat), 1);
   res_d = p.inverse()*mat_d;
   VERIFY(res.isApprox(res_d) && "inv(p)*mat");
 
+  // test non-plaintype expressions that require additional temporary
+  const Scalar alpha(2.34);
+
+  res_d = p * (alpha * mat_d);
+  VERIFY_TEMPORARY_COUNT( res = p * (alpha * mat), 2);
+  VERIFY( res.isApprox(res_d) && "p*(alpha*mat)" );
+
+  res_d = (alpha * mat_d) * p;
+  VERIFY_TEMPORARY_COUNT( res = (alpha * mat) * p, 2);
+  VERIFY( res.isApprox(res_d) && "(alpha*mat)*p" );
+
+  res_d = p.inverse() * (alpha * mat_d);
+  VERIFY_TEMPORARY_COUNT( res = p.inverse() * (alpha * mat), 2);
+  VERIFY( res.isApprox(res_d) && "inv(p)*(alpha*mat)" );
+
+  res_d = (alpha * mat_d) * p.inverse();
+  VERIFY_TEMPORARY_COUNT( res = (alpha * mat) * p.inverse(), 2);
+  VERIFY( res.isApprox(res_d) && "(alpha*mat)*inv(p)" );
+
+  //
+
   VERIFY( is_sorted( (p * mat * p.inverse()).eval() ));
   VERIFY( is_sorted( res = mat.twistedBy(p) ));
   VERIFY_TRANSPOSITION_COUNT( ::eval(p * mat * p.inverse()), 0);
diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h
index cf5828e..5c26133 100644
--- a/unsupported/Eigen/src/SparseExtra/MarketIO.h
+++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h
@@ -189,6 +189,7 @@
         readsizes = true;
         mat.resize(M,N);
         mat.reserve(NNZ);
+        elements.reserve(NNZ);
       }
     }
     else