dyn.hpp

Go to the documentation of this file.

//=======================================================================
// Copyright (c) 2014-2023 Baptiste Wicht
// Distributed under the terms of the MIT License.
// (See accompanying file LICENSE or copy at
//  http://opensource.org/licenses/MIT)
//=======================================================================

#pragma once

#include "etl/dyn_base.hpp"    //The base class and utilities
#include "etl/direct_fill.hpp" //direct_fill with GPU support

namespace etl {

template <typename T, order SO, size_t D>
struct dyn_matrix_impl final : dense_dyn_base<dyn_matrix_impl<T, SO, D>, T, SO, D>,
                               inplace_assignable<dyn_matrix_impl<T, SO, D>>,
                               expression_able<dyn_matrix_impl<T, SO, D>>,
                               value_testable<dyn_matrix_impl<T, SO, D>>,
                               iterable<dyn_matrix_impl<T, SO, D>, SO == order::RowMajor>,
                               dim_testable<dyn_matrix_impl<T, SO, D>> {
    static constexpr size_t n_dimensions = D;                                      
    static constexpr order storage_order = SO;                                     
    static constexpr size_t alignment    = default_intrinsic_traits<T>::alignment; 

    using this_type              = dyn_matrix_impl<T, SO, D>;                  
    using base_type              = dense_dyn_base<this_type, T, SO, D>;        
    using iterable_base_type     = iterable<this_type, SO == order::RowMajor>; 
    using value_type             = T;                                          
    using dimension_storage_impl = std::array<size_t, n_dimensions>;           
    using memory_type            = value_type*;                                
    using const_memory_type      = const value_type*;                          

    using iterator       = std::conditional_t<SO == order::RowMajor, value_type*, etl::iterator<this_type>>;             
    using const_iterator = std::conditional_t<SO == order::RowMajor, const value_type*, etl::iterator<const this_type>>; 

    template <typename V = default_vec>
    using vec_type = typename V::template vec_type<T>;

private:
    using base_type::_dimensions;
    using base_type::_memory;
    using base_type::_size;

    using base_type::allocate;
    using base_type::check_invariants;
    using base_type::release;

public:
    using base_type::dim;
    using base_type::memory_end;
    using base_type::memory_start;
    using iterable_base_type::begin;
    using iterable_base_type::end;

    // Construction

    dyn_matrix_impl() noexcept : base_type() {
        //Nothing else to init
    }

    dyn_matrix_impl(const dyn_matrix_impl& rhs) noexcept(assert_nothrow) : base_type(rhs) {
        if (_size) {
            _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));
        } else {
            _memory = nullptr;
        }

        if (rhs.is_cpu_up_to_date()) {
            direct_copy(rhs.memory_start(), rhs.memory_end(), memory_start());
        }

        cpp_assert(rhs.is_cpu_up_to_date() == this->is_cpu_up_to_date(), "dyn_matrix_impl(&) must preserve CPU status");
        cpp_assert(rhs.is_gpu_up_to_date() == this->is_gpu_up_to_date(), "dyn_matrix_impl(&) must preserve GPU status");
    }

    dyn_matrix_impl(dyn_matrix_impl&& rhs) noexcept : base_type(std::move(rhs)) {
        _memory     = rhs._memory;

        rhs._size   = 0;
        rhs._memory = nullptr;
    }

    dyn_matrix_impl(std::initializer_list<value_type> list) noexcept requires(n_dimensions == 1) : base_type(list.size(), {{list.size()}}) {
        _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));

        std::copy(list.begin(), list.end(), begin());
    }

    template <size_c... S>
    explicit dyn_matrix_impl(S... sizes) noexcept requires(sizeof...(S) == D) : base_type(util::size(sizes...), {{static_cast<size_t>(sizes)...}}) {
        _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));
    }

    template <typename... S>
    explicit dyn_matrix_impl(S... sizes) noexcept requires(dyn_detail::is_initializer_list_constructor<S...>::value && sizeof...(S) == D + 1)
            : base_type(util::size(std::make_index_sequence<(sizeof...(S) - 1)>(), sizes...),
                        dyn_detail::sizes(std::make_index_sequence<(sizeof...(S) - 1)>(), sizes...)) {
        _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));

        static_assert(sizeof...(S) == D + 1, "Invalid number of dimensions");

        auto list = cpp::last_value(sizes...);
        std::copy(list.begin(), list.end(), begin());
    }

    template <typename... S>
    explicit dyn_matrix_impl(size_t s1, S... sizes) noexcept requires((sizeof...(S) == D) && cpp::specialization_of<values_t, typename cpp::last_type<size_t, S...>::type>)
            : base_type(util::size(std::make_index_sequence<(sizeof...(S))>(), s1, sizes...),
                        dyn_detail::sizes(std::make_index_sequence<(sizeof...(S))>(), s1, sizes...)) {
        _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));

        auto list = cpp::last_value(sizes...).template list<value_type>();
        std::copy(list.begin(), list.end(), begin());
    }

    template <typename... S>
    explicit dyn_matrix_impl(size_t s1, S... sizes) noexcept requires((sizeof...(S) == D) && !cpp::specialization_of<values_t, typename cpp::last_type<size_t, S...>::type>)
            : base_type(util::size(std::make_index_sequence<(sizeof...(S))>(), s1, sizes...),
                        dyn_detail::sizes(std::make_index_sequence<(sizeof...(S))>(), s1, sizes...)) {
        _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));

        decltype(auto) value = cpp::last_value(s1, sizes...);
        std::fill(begin(), end(), value);
    }

    template <std_container Container>
    explicit dyn_matrix_impl(const Container& container) requires(D == 1 && std::convertible_to<typename Container::value_type, value_type>) :
            base_type(container.size(), {{container.size()}}) {
        _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));

        // Copy the container directly inside the allocated memory
        std::copy_n(container.begin(), _size, _memory);
    }

    dyn_matrix_impl& operator=(const dyn_matrix_impl& rhs) noexcept(assert_nothrow) {
        if (this != &rhs) {
            if (!_size) {
                cpp_assert(!_memory, "_size and _memory are not properly synced");

                _size       = rhs._size;
                _dimensions = rhs._dimensions;

                if (_size) {
                    _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));
                } else {
                    _memory = nullptr;
                }
            } else {
                validate_assign(*this, rhs);
            }

            // TODO Find a better solution
            const_cast<dyn_matrix_impl&>(rhs).assign_to(*this);
        }

        check_invariants();

        cpp_assert(rhs.is_cpu_up_to_date() == this->is_cpu_up_to_date(), "dyn_matrix_impl::operator= must preserve CPU status");
        cpp_assert(rhs.is_gpu_up_to_date() == this->is_gpu_up_to_date(), "dyn_matrix_impl::operator= must preserve GPU status");

        return *this;
    }

    dyn_matrix_impl& operator=(dyn_matrix_impl&& rhs) noexcept {
        if (this != &rhs) {
            if (_memory) {
                release(_memory, alloc_size_mat<T>(_size, dim(n_dimensions - 1)));
            }

            _size       = rhs._size;
            _dimensions = std::move(rhs._dimensions);
            _memory     = rhs._memory;

            rhs._size   = 0;
            rhs._memory = nullptr;
        }

        check_invariants();

        return *this;
    }

    void resize_arr(const dimension_storage_impl& dimensions) {
        auto new_size = std::accumulate(dimensions.begin(), dimensions.end(), size_t(1), std::multiplies<size_t>());

        if (_memory) {
            auto new_memory = allocate(alloc_size_mat<T>(new_size, (dimensions.back())));

            for (size_t i = 0; i < std::min(_size, new_size); ++i) {
                new_memory[i] = _memory[i];
            }

            release(_memory, alloc_size_mat<T>(_size, dim(n_dimensions - 1)));

            _memory = new_memory;
        } else {
            _memory = allocate(alloc_size_mat<T>(new_size, (dimensions.back())));
        }

        _size       = new_size;
        _dimensions = dimensions;
    }

    template <typename... Sizes>
    void resize(Sizes... sizes) requires(sizeof...(Sizes) == n_dimensions) {
        static_assert(sizeof...(Sizes), "Cannot change number of dimensions");

        auto new_size = util::size(sizes...);

        if (_memory) {
            auto new_memory = allocate(alloc_size_mat<T>(new_size, cpp::last_value(sizes...)));

            for (size_t i = 0; i < std::min(_size, new_size); ++i) {
                new_memory[i] = _memory[i];
            }

            release(_memory, alloc_size_mat<T>(_size, dim(n_dimensions - 1)));

            _memory = new_memory;
        } else {
            _memory = allocate(alloc_size_mat<T>(new_size, cpp::last_value(sizes...)));
        }

        _size       = new_size;
        _dimensions = dyn_detail::sizes(std::make_index_sequence<D>(), sizes...);
    }

    template <etl_expr E>
    dyn_matrix_impl& operator=(E&& e) noexcept requires(!std::is_same_v<std::decay_t<E>, dyn_matrix_impl<T, SO, D>> && std::convertible_to<value_t<E>, value_type>) {
        // It is possible that the matrix was not initialized before
        // In the case, get the the dimensions from the expression and
        // initialize the matrix
        if (!_memory) {
            inherit(e);
        } else {
            validate_assign(*this, e);
        }

        // Avoid aliasing issues
        if constexpr (!decay_traits<E>::is_linear) {
            if (e.alias(*this)) {
                // Create a temporary to hold the result
                this_type tmp;

                // Get the proper size
                tmp.inherit(e);

                // Assign the expression to the temporary
                e.assign_to(tmp);

                // Assign the temporary to this matrix
                *this = tmp;
            } else {
                e.assign_to(*this);
            }
        } else {
            // Direct assignment of the expression into this matrix
            e.assign_to(*this);
        }

        check_invariants();

        return *this;
    }

    template <std_container Container>
    dyn_matrix_impl& operator=(const Container& vec) requires(std::convertible_to<typename Container::value_type, value_type>) {
        // Inherit from the dimensions if possible
        if (!_memory && D == 1) {
            // Compute the size and new dimensions
            _size          = vec.size();
            _dimensions[0] = vec.size();

            // Allocate the new memory
            _memory = allocate(alloc_size_mat<T>(_size, dim(0)));
        } else {
            validate_assign(*this, vec);
        }

        std::copy(vec.begin(), vec.end(), begin());

        this->validate_cpu();
        this->invalidate_gpu();

        check_invariants();

        return *this;
    }

    dyn_matrix_impl& operator=(const value_type& value) noexcept {
        direct_fill(*this, value);

        check_invariants();

        return *this;
    }

    ~dyn_matrix_impl() noexcept {
        if (_memory) {
            release(_memory, alloc_size_mat<T>(_size, dim(n_dimensions - 1)));
        }
    }

    void clear() {
        if (_memory) {
            release(_memory, alloc_size_mat<T>(_size, dim(n_dimensions - 1)));
        }

        _memory = nullptr;
        _size   = 0;
    }

    template <typename Y>
    auto& gpu_compute_hint([[maybe_unused]] Y& y) {
        this->ensure_gpu_up_to_date();
        return *this;
    }

    template <typename Y>
    const auto& gpu_compute_hint([[maybe_unused]] Y& y) const {
        this->ensure_gpu_up_to_date();
        return *this;
    }

    void swap(dyn_matrix_impl& other) {
        using std::swap;
        swap(_size, other._size);
        swap(_dimensions, other._dimensions);
        swap(_memory, other._memory);

        //TODO swap is likely screwing up GPU memory!

        check_invariants();
    }

    template <typename V = default_vec>
    ETL_STRONG_INLINE(void)
    store(const vec_type<V> in, size_t i) noexcept {
        V::store(_memory + i, in);
    }

    template <typename V = default_vec>
    ETL_STRONG_INLINE(void)
    storeu(const vec_type<V> in, size_t i) noexcept {
        V::storeu(_memory + i, in);
    }

    template <typename V = default_vec>
    ETL_STRONG_INLINE(void)
    stream(const vec_type<V> in, size_t i) noexcept {
        V::stream(_memory + i, in);
    }

    template <typename V = default_vec>
    ETL_STRONG_INLINE(vec_type<V>)
    load(size_t i) const noexcept {
        return V::load(_memory + i);
    }

    template <typename V = default_vec>
    ETL_STRONG_INLINE(vec_type<V>)
    loadu(size_t i) const noexcept {
        return V::loadu(_memory + i);
    }

    size_t& unsafe_dimension_access(size_t i) {
        cpp_assert(i < n_dimensions, "Out of bounds");
        return _dimensions[i];
    }

    template <exact_dimensions<D> E>
    void inherit_if_null(const E& e) {
        if (!_memory) {
            inherit(e);
        }
    }

    // Assignment functions

    template <typename L>
    void assign_to(L&& lhs) const {
        std_assign_evaluate(*this, std::forward<L>(lhs));
    }

    template <typename L>
    void assign_add_to(L&& lhs) const {
        std_add_evaluate(*this, std::forward<L>(lhs));
    }

    template <typename L>
    void assign_sub_to(L&& lhs) const {
        std_sub_evaluate(*this, std::forward<L>(lhs));
    }

    template <typename L>
    void assign_mul_to(L&& lhs) const {
        std_mul_evaluate(*this, std::forward<L>(lhs));
    }

    template <typename L>
    void assign_div_to(L&& lhs) const {
        std_div_evaluate(*this, std::forward<L>(lhs));
    }

    template <typename L>
    void assign_mod_to(L&& lhs) const {
        std_mod_evaluate(*this, std::forward<L>(lhs));
    }

    // Internals

    void visit([[maybe_unused]] const detail::evaluator_visitor& visitor) const {}

private:
    template <typename E>
    void inherit([[maybe_unused]] const E& e) {
        if constexpr (etl::decay_traits<E>::is_generator) {
            cpp_unreachable("Impossible to inherit dimensions from generators");
        } else {
            cpp_assert(n_dimensions == etl::dimensions(e), "Invalid number of dimensions");

            // Compute the size and new dimensions
            _size = 1;
            for (size_t d = 0; d < n_dimensions; ++d) {
                _dimensions[d] = etl::dim(e, d);
                _size *= _dimensions[d];
            }

            // Allocate the new memory
            _memory = allocate(alloc_size_mat<T>(_size, dim(n_dimensions - 1)));
        }
    }

    friend std::ostream& operator<<(std::ostream& os, const dyn_matrix_impl& mat) {
        if (D == 1) {
            return os << "V[" << mat.size() << "]";
        }

        os << "M[" << mat.dim(0);

        for (size_t i = 1; i < D; ++i) {
            os << "," << mat.dim(i);
        }

        return os << "]";
    }
};

#ifndef CPP_UTILS_ASSERT_EXCEPTION
static_assert(std::is_nothrow_default_constructible_v<dyn_vector<double>>, "dyn_vector should be nothrow default constructible");
static_assert(std::is_nothrow_copy_constructible_v<dyn_vector<double>>, "dyn_vector should be nothrow copy constructible");
static_assert(std::is_nothrow_move_constructible_v<dyn_vector<double>>, "dyn_vector should be nothrow move constructible");
static_assert(std::is_nothrow_copy_assignable_v<dyn_vector<double>>, "dyn_vector should be nothrow copy assignable");
static_assert(std::is_nothrow_move_assignable_v<dyn_vector<double>>, "dyn_vector should be nothrow move assignable");
static_assert(std::is_nothrow_destructible_v<dyn_vector<double>>, "dyn_vector should be nothrow destructible");
#endif

template <typename T, typename... Sizes>
etl::dyn_matrix<T, sizeof...(Sizes)> make_dyn_matrix(Sizes... sizes) {
    return etl::dyn_matrix<T, sizeof...(Sizes)>(sizes...);
}

template <typename T, order SO, size_t D>
void swap(dyn_matrix_impl<T, SO, D>& lhs, dyn_matrix_impl<T, SO, D>& rhs) {
    lhs.swap(rhs);
}

template <typename Stream, typename T, order SO, size_t D>
void serialize(serializer<Stream>& os, const dyn_matrix_impl<T, SO, D>& matrix) {
    for (size_t i = 0; i < etl::dimensions(matrix); ++i) {
        os << matrix.dim(i);
    }

    for (const auto& value : matrix) {
        os << value;
    }
}

template <typename Stream, typename T, order SO, size_t D>
void deserialize(deserializer<Stream>& is, dyn_matrix_impl<T, SO, D>& matrix) {
    typename std::decay_t<decltype(matrix)>::dimension_storage_impl new_dimensions;

    for (auto& value : new_dimensions) {
        is >> value;
    }

    matrix.resize_arr(new_dimensions);

    for (auto& value : matrix) {
        is >> value;
    }
}

} //end of namespace etl