hh check

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

#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 & _rhs)
    {
        auto irow = std::begin(matrix_);
        for (auto const & row_ : _rhs) {
            auto icol = std::begin(*irow);
            for (auto const & v : row_) {
                *icol = v;
                ++icol;
            }
            ++irow;
        }
    }

    template< typename output, typename input >
    void
    mul(output & _output, input const & _input) const
    {
        for (size_type i = 0; i < dimension_; ++i) {
            for (size_type j = 0; j < dimension_; ++j) {
                value_type & result_ = _output[i][j];
                result_ = zero;
                for (size_type k = 0; k < dimension_; ++k) {
                    result_ += _input[i][k] * matrix_[k][j];
                }
            }
        }
    }

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

    template< typename output >
    bool
    orthonormalize(output & _output,
                   size_type const _rank)
    {
        assert(!(dimension_ < _rank));
        for (size_type i = 0; i < _rank; ++i) { // Householder transformation
            value_type sum_ = zero;
            vector & qri_ = matrix_[i];
            for (size_type j = i; j < dimension_; ++j) {
                value_type const & qrij_ = qri_[j];
                sum_ += qrij_ * qrij_;
            }
            using std::sqrt;
            value_type norm_ = sqrt(sum_);
            if (!(eps < norm_)) {
                return false;
            }
            value_type & qrii_ = qri_[i];
            if (qrii_ < zero) {
                qrii_ = -qrii_;
            }
            value_type const factor_ = sqrt(std::move(sum_) + norm_ * qrii_);
            if (!(eps < factor_)) {
                return false;
            }
            qrii_ += std::move(norm_);
            for (size_type k = i; k < dimension_; ++k) {
                qri_[k] /= factor_;
            }
            size_type j = i;
            while (++j < _rank) {
                vector & qrj_ = matrix_[j];
                value_type s_ = zero;
                for (size_type k = i; k < dimension_; ++k) {
                    s_ += qri_[k] * qrj_[k];
                }
                for (size_type k = i; k < dimension_; ++k) {
                    qrj_[k] -= qri_[k] * s_;
                }
            }
        } // matrix_ is packed QR
        for (size_type i = 0; i < _rank; ++i) {
            vector & qi_ = _output[i];
            qi_ = zero;
            qi_[i] = one;
            size_type j = _rank;
            while (0 < j) {
                --j;
                vector const & qrj_ = matrix_[j]; // containing packed QR
                value_type s_ = zero;
                for (size_type k = j; k < dimension_; ++k) {
                    s_ += qrj_[k] * qi_[k];
                }
                for (size_type k = j; k < dimension_; ++k) {
                    qi_[k] -= qrj_[k] * s_;
                }
            }
        } // _output is Q
        return true;
    }

    template< typename output >
    bool
    orthonormalize(output & _output)
    {
        return orthonormalize(_output, dimension_);
    }

    template< typename input >
    void
    matrix_transpose(input & _input)
    {
        for (size_type r = 0; r < dimension_; ++r) {
            for (size_type c = 0; c < dimension_; ++c) {
                std::swap(_input[c][r], _input[r][c]);
            }
        }
    }

    void
    matrix_transpose()
    {
        return matrix_transpose(matrix_);
    }

    template< typename output >
    void
    matrix_sqr(output & _output) const
    {
        for (size_type r = 0; r < dimension_; ++r) {
            vector & lhs_ = _output[r];
            vector const & row_ = matrix_[r];
            for (size_type c = 0; c < dimension_; ++c) {
                lhs_[c] = std::inner_product(std::cbegin(row_), std::cend(row_), std::cbegin(matrix_[c]), zero);
            }
        }
    }

    operator matrix const & () const
    {
        return matrix_;
    }

};

// main.cpp
#include <iostream>

#include <cstdlib>

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

inline
std::ostream &
operator << (std::ostream & _out, matrix const & _matrix)
{
    for (vector const & vector_ : _matrix) {
        for (value_type const & value_ : vector_) {
            if (value_ < 10 * std::numeric_limits< value_type >::epsilon()) {
                _out << 0 << ' ';
            } else {
                _out << value_ << ' ';
            }
        }
        _out << std::endl;
    }
    return _out;
}

int
main()
{
    value_type const eps = std::numeric_limits< value_type >::epsilon();
    std::size_t const dimension_ = 5;
    math< value_type > m(dimension_, eps);
    std::random_device rd;
    std::default_random_engine e(rd());
    std::normal_distribution< value_type > d(m.one, m.one);
    matrix input_;
    input_.resize(dimension_);
    for (vector & vector_ : input_) {
        vector_.resize(dimension_);
        for (value_type & value_ : vector_) {
            value_ = d(e);
        }
    }
    std::cout << "input = \n" << input_ << std::endl;
    matrix output_;
    output_.resize(dimension_);
    for (vector & vector_ : output_) {
        vector_.resize(dimension_, m.zero);
    }
    m = input_;
    assert(m.orthonormalize(output_));
    std::cout << "output = \n" << output_ << std::endl;
    m = output_;
    m.matrix_sqr(output_);
    std::cout << "output = \n" << output_ << std::endl;
    return EXIT_SUCCESS;
}