Anyone can explain a bit about permutation_matrix class in lu.hpp, please. I assume that it stores value i at position i.

Current version has an inconvenience: I would need a resize() method to reuse (when possible) the underlying memory of its vector base class.

Maybe a a simpler (minimalistic) approach can be used:

template <typename T = std::size_t>
class permutation_matrix
{
public:
    permutation_matrix(T size) : m_size(size) {}
    T operator() (T i) const {return i;} // check size in debug mode
    T size() const {return m_size();}
private:
    T m_size;
};

As only const operator() is used, there is no need to actually store data i at position i, right?