determinant based on LUP-decomposition

#include <valarray>
#include <vector>
#include <utility>

#include <cmath>
#include <cstddef>
#include <cassert>

template< typename value_type >
struct math
{

    using size_type = std::size_t;

    size_type const dimension_;
    value_type const & eps;

    value_type const zero = value_type(0);
    value_type const one = value_type(1);

private :

    using vector = std::valarray< value_type >;
    using matrix = std::vector< vector >;

    matrix matrix_;

    value_type
    det(matrix & _matrix,
        size_type const _dimension) // hottest function
    { // calculates lower unit triangular matrix and upper triangular
        assert(0 < _dimension);
        value_type det_ = one;
        for (size_type i = 0; i < _dimension; ++i) {
            vector & mi_ = _matrix[i];
            size_type pivot = i;
            {
                using std::abs;
                value_type max_ = abs(mi_[i]);
                size_type j = i;
                while (++j < _dimension) {
                    value_type y_ = abs(_matrix[j][i]);
                    if (max_ < y_) {
                        max_ = std::move(y_);
                        pivot = j;
                    }
                }
                if (!(eps < max_)) { // regular?
                    return zero; // singular
                }
            }
            if (pivot != i) {
                det_ = -det_; // each permutation flips sign of det
                mi_.swap(_matrix[pivot]);
            }
            value_type const & dia_ = mi_[i];
            det_ *= dia_; // det is multiple of diagonal elements
            size_type j = i;
            while (++j < _dimension) {
                vector & mj_ = _matrix[j];
                value_type & mji_ = mj_[i];
                mji_ /= dia_;
                size_type k = i;
                while (++k < _dimension) {
                    mj_[k] -= mji_ * mi_[k];
                }
            }
        }
        return det_;
    }

public :

    math(size_type const _dimension,
         value_type const & _eps)
        : dimension_(_dimension)
        , eps(_eps)
        , matrix_(dimension_)
    {
        assert(1 < dimension_);
        assert(!(eps < zero));
        for (size_type r = 0; r < dimension_; ++r) {
            matrix_[r].resize(dimension_);
        }
    }

    template< typename rhs = matrix >
    void
    operator = (rhs const & _matrix)
    {
        auto irow = std::begin(matrix_);
        for (auto const & row_ : _matrix) {
            auto icol = std::begin(*irow);
            for (auto const & v : row_) {
                *icol = v;
                ++icol;
            }
            ++irow;
        }
    }

    value_type
    det()
    {
        return det(matrix_, dimension_);
    }

};

// main.cpp
#include <iostream>

#include <cstdlib>

int
main()
{
    using value_type = double;
    value_type const eps = std::numeric_limits< value_type >::epsilon();
    std::size_t const dimension_ = 3;
    math< value_type > m(dimension_, eps);
    m = { // example from https://en.wikipedia.org/wiki/Determinant#Laplace.27s_formula_and_the_adjugate_matrix
            {-2.0, 2.0, -3.0},
            {-1.0, 1.0,  3.0},
            { 2.0, 0.0, -1.0}
        };
    std::cout << m.det() << std::endl; // 18
    return EXIT_SUCCESS;
}